Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Closed
TorkelE opened this issue Sep 24, 2023 · 8 comments

Comments

@TorkelE
Copy link

TorkelE commented Sep 24, 2023

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
Copy link
Member

rveltz commented Sep 25, 2023

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
Copy link
Author

TorkelE commented Sep 25, 2023

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

@rveltz
Copy link
Member

rveltz commented Sep 25, 2023

it is tagged, so v> 0.3

@TorkelE
Copy link
Author

TorkelE commented Sep 25, 2023

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
Copy link
Author

TorkelE commented Sep 25, 2023

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
Copy link
Member

rveltz commented Sep 25, 2023

excellent!

@rveltz
Copy link
Member

rveltz commented Feb 3, 2024

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

@TorkelE
Copy link
Author

TorkelE commented Feb 3, 2024

Yes, let's do

@TorkelE TorkelE closed this as completed Feb 3, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants