bifurcationkit / BifurcationKit.jl

A Julia package to perform Bifurcation Analysis
https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable
Other
288 stars 34 forks source link

Add `BifurcationProblem` dispatch for ModelingToolkit problems (possibly in an extension?) #118

Closed TorkelE closed 5 months ago

TorkelE commented 9 months ago

Currently, when I want to e.g. compute a bifurcation diagram for a ModelingToolkit-type structure (as e.g. generated by Catalyst) I have to manually extract the f and jac functions, e.g:

using Catalyst
rn = @reaction_network begin
    (v0 + v*(S * X)^n / ((S*X)^n + (D*A)^n + K^n), d), ∅ ↔ X
    (X/τ, 1/τ), ∅ ↔ A
end

p = Dict(:S => 1., :D => 9., :τ => 1000., :v0 => 0.01,
         :v => 2., :K => 20., :n => 3, :d => 0.05)
bif_par = :S           # bifurcation parameter
p_span = (0.1, 20.)    # interval to vary S over
plot_var = :X          # we will plot X vs S

p_bstart = copy(p)
p_bstart[bif_par] = p_span[1]
u0 = [:X => 1.0, :A => 1.0]

oprob = ODEProblem(rn, u0, (0.0, 0.0), p_bstart; jac = true)
F = (u,p) -> oprob.f(u, p, 0)
J = (u,p) -> oprob.f.jac(u, p, 0)

# get S and X as symbolic variables
@unpack S, X = rn

# find their indices in oprob.p and oprob.u0 respectively
bif_idx  = findfirst(isequal(S), parameters(rn))
plot_idx = findfirst(isequal(X), species(rn))

using BifurcationKit, Plots, LinearAlgebra, Setfield

bprob = BifurcationProblem(F, oprob.u0, oprob.p, (@lens _[bif_idx]);
                           recordFromSolution = (x, p) -> x[plot_idx], J = J)

It seems like it would be easier to simply create a dispatch for BifurcationProblem which automatically finds the F and J functions? If this sounds useful, I am happy to make a PR to BifurcationKit creating an extension (loaded if e.g. ModelingToolkit is loaded) that adds this dispatch.

Possibly, it could also take the bifurcation parameter (and maybe under some circumstances, the plotting variable), and automate the bif_idx and plot_idx computations as well.

rveltz commented 9 months ago

Yes, that would be very nice if you could do that! Note that BK@0.3.0 has a new syntax:

bprob = BifurcationProblem(F, oprob.u0, oprob.p, (@lens _[bif_idx]);
                           record_from_solution = (x, p) -> x[plot_idx], J = J)

Also, no need to load Setfield, it is re-exported from BifurcationKit

TorkelE commented 9 months ago

Sounds good, I presume this is the syntax on master? I will make a PR as soon as it is finished.

rveltz commented 9 months ago

it is tagged, so v> 0.3

TorkelE commented 9 months ago

Yeah, should have double checked that properly. Congratulations on the 0.3 release, looking forward to try around with it. Once this new interface is ready I will rewrite the tutorial over at Catalyst.

TorkelE commented 9 months ago

Also, given that BifurcationKit already has SciMLBase (which contains the ModelingToolkit types) in its dependencies, creating an extension for the new dispatches would be redundant. I will just add it to the base package.

rveltz commented 9 months ago

excellent!

rveltz commented 5 months ago

should we close this? It seems it is solved on MTK side, right?

TorkelE commented 5 months ago

Yes, let's do