(potentially long-term) Complete integration with ModelingTookit's nonlinear problem

[TorkelE]

It is currently possible to define a set of non-linear equations in ModelingTookit, and use these as input to solve a bifurcation problem (https://docs.sciml.ai/ModelingToolkit/stable/tutorials/bifurcation_diagram_computation/). This depends on an extension which converts MTK's NonlinearProblem to a f and a J function that BifurcationKit can interpret internally. Would it make sense to actually store the NonlinearProblem all the way through BifurcationKit's internal process?

I could see some advantages with this. Some really mundane, like parameter's "names" being stored, so these can easier for labeling plot x-axises with parameters' names. Furthermore, if the full bifurcation path is stored, SymbolicIndexingInterface could be used to get the output path's values for a specific variable, e.g.

br = continuation(prob, ...)
br[:X] #gives the full path of `X` values.
br[:Y] #gives the full path of `Y` values.

Also, NonlinearSolve.jl has recently seen a lot of work (https://arxiv.org/abs/2403.16341). Here, BifurcationKit could directly leverage the Newton (and other nonlinear solvers) implemented in NonlinearSolve.jl (it would be even cooler if NonlinearSolve.jl could utilise features implemented in BifurcationKit only, like deflated Newton).

It might require some work, but I think it is worth a thought, I think it could turn out really nice.

[rveltz]

You can already have this feature br.X when X is in record_from_solution

[rveltz]

As for NonlinearSolve, I am thinking about this with @Datseris
My main blocking point for now is that I want to allow a vector inferface which is not an AbstractArray like KrylovKit do.

[TorkelE]

You mean defining values in plain vector form like

u = [1.0, 2.0, 3.0]

instead of a map

u = [X => 1.0, Y => 2.0, Z => 3.0]

?

[rveltz]

No, define

prob = BifurcationProblem(..., record_from_solution = (x,p) -> (X = x[1], Y = x[2]))

[TorkelE]

Right, that makes sense. Yeah, that would require some modification (or supporting alternative approaches simultaneously, but that also seems messy).

[rveltz]

The biggest advantage would be to tackle parameter derivatives not with finite differences

[TorkelE]

Actually, I checked again and vector interface should be supported for NonlinearSystems. Presumable, in your example, x would end up a SciML integrator object. For these, both integrator[X] and integrator.u[1] should be supported (Furthermore, integrator[1] would be easy to support, if you'd need that).

[rveltz]

Or really?? you mean this https://jutho.github.io/KrylovKit.jl/stable/ and VectorInterface.jl ?

[TorkelE]

Sorry, I think I am missunderstanding something. Not familiar with those packages, so there are probably stuff going on that I don't really understand.

[rveltz]

Should we close this? It is seems that BK is a now an external dep of MTK

[Datseris]

huh? where? how?

[rveltz]

Sorry I misread.

Would it make sense to actually store the NonlinearProblem all the way through BifurcationKit's internal process?

You could save the NonlinearProblem in your own BifurcationProblem, it has to be a subtype of BK.AbstractBifurcationProblem. You could then dispatch _newton in your problem type to use NonlinearSolve. For newton_palc, it is also possible but a tiny bit more technical. It becomes difficult if you want to target large dim and small dim in one interface and the best is then to rely on Bordered Linear solvers which are not in LinearSolve.jl but present in BifurcationKit

You could also avoid using an Accessors.@optic for the parameters and just provide getp