BifKit extension, conservation laws, and changing initial conditions

[isaacsas]

We need to check if it is actually working when one is using an initial condition as the bifurcation parameter. I suspect that it has no way to know it needs to update the conserved constant in this case (i.e. it is just changing that parameter in the underlying parameter vector), and hence will give incorrect results.

[isaacsas]

This doesn't work, but at least it errors:

rn = @reaction_network begin
              @parameters M
              @species A(t) = M
              (1.0,1.0), A <--> B
              end

osys = complete(convert(ODESystem, rn; remove_conserved = true))
u_guess = [osys.B => 2.0]
p_start = [osys.M => 3.0]
bprob = BifurcationProblem(osys, u_guess, p_start, osys.M; plot_var = osys.A)

giving

ERROR: ArgumentError: SymbolicUtils.BasicSymbolic{Real}[M] are either missing from the variable map or missing from the system's unknowns/parameters list.
Stacktrace:
 [1] throw_missingvars_in_sys(vars::Vector{SymbolicUtils.BasicSymbolic{Real}})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/XHRvt/src/utils.jl:765
 [2] promote_to_concrete(vs::Vector{SymbolicUtils.BasicSymbolic{Real}}; tofloat::Bool, use_union::Bool)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/XHRvt/src/utils.jl:784
 [3] varmap_to_vars(varmap::Vector{…}, varlist::Vector{…}; defaults::Dict{…}, check::Bool, toterm::Function, promotetoconcrete::Nothing, tofloat::Bool, use_union::Bool)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/XHRvt/src/variables.jl:171
 [4] BifurcationProblem(::NonlinearSystem, ::Vector{…}, ::Vector{…}, ::Num; plot_var::Num, record_from_solution::Function, jac::Bool, kwargs::@Kwargs{})
   @ MTKBifurcationKitExt ~/.julia/packages/ModelingToolkit/XHRvt/ext/MTKBifurcationKitExt.jl:104
 [5] BifurcationProblem(::ODESystem, ::Vector{Pair{Num, Float64}}, ::Vararg{Any}; kwargs::@Kwargs{plot_var::Num})
   @ MTKBifurcationKitExt ~/.julia/packages/ModelingToolkit/XHRvt/ext/MTKBifurcationKitExt.jl:152
 [6] top-level scope
   @ REPL[35]:1
Some type information was truncated. Use `show(err)` to see complete types.

[isaacsas]

Continuing with respect to Γ[1], the conservation constant, seems to (silently) give wrong answers:

using Catalyst, BifurcationKit, Plots
rn = @reaction_network begin
              (1.0,1.0), A <--> B
              end
osys = structural_simplify(convert(ODESystem, rn; remove_conserved = true))
u_guess = [osys.A => 1.0, osys.B => 2.0]
p_start = [osys.Γ[1] => 3.0]  # A = 1.0 then
bprob = BifurcationProblem(osys, u_guess, p_start, osys.Γ[1]; plot_var = osys.A)
p_span = (3.0, 20.) # A goes from 1.0 to 18.0
opts_br = ContinuationPar(p_min = p_span[1], p_max = p_span[2], max_steps = 1000)
bif_dia = bifurcationdiagram(bprob, PALC(), 2, opts_br; bothside = true)
plot(bif_dia; xlabel = "A(0) + B(0)", label = "A(t)")

which gives:

Actually this is right and I'm just being dumb. So that is good.

[isaacsas]

Seems like we need to disallow changing initial conditions and/or the conservation constant in the current extension setup.

[isaacsas]

@rveltz is there any way you could use remake for updating parameters and initial conditions in the bifurcation solvers (at least when using problems coming via MTK)? I guess right now you just modify the bifurcation parameter, but in cases like this one needs to hook into the SciML remake functionality for everything to update correctly. If you have an initial condition given by a symbolic parameter you want to make the bifurcation parameter, or a conservation law over which you want to vary the conserved constant, you need to have any other dependent parameters/ICs also update accordingly.

[rveltz]

Do you want the same name because there is BK.re_make

I was worried about type piracy so that s why the new name. I can name it as in SciMLBase if needed

[isaacsas]

Does that hook into the symbolic SciMLBase.remake?

[rveltz]

no, it does not extend it

[isaacsas]

Yeah, so I suspect that if one has parameters that are functions of initial conditions, initial conditions that are functions of parameters, or parameters that are functions of other parameters, one may be likely to get incorrect results (because these may need to get updated as the bifurcation parameter changes).

[rveltz]

If this https://github.com/SciML/SciMLBase.jl/blob/0f8ec1f501e916e5f8c39d94ddc02757d574cc9c/src/remake.jl#L27 is enough, it wont take me long.

But if you want this one

https://github.com/SciML/SciMLBase.jl/blob/0f8ec1f501e916e5f8c39d94ddc02757d574cc9c/src/remake.jl#L80

then I have to think

[isaacsas]

Maybe this would be easier by having a reeval type function that a user can pass in (in this case Catalyst/MTK)? That could handle updating all other parameters when one changes during a solve? Then we could perhaps pass a wrapped version of SciMLBase.remake? @AayushSabharwal is the expert on this functionality in SciMLBase, but it would be nice to get this working and updating correctly.

[isaacsas]

I guess trying to worry about changing a parameter in an initial condition is a bit much, but I'll try to work up a simple example where there are dependent parameters but things aren't updating correctly when the bifurcation parameter changes.

[rveltz]

OK!

[vyudu]

    rn = @reaction_network begin
        @parameters k ksq = k^2
        (k, ksq), A <--> B
    end

    rn = complete(rn)
    u0_guess = [:A => 1., :B => 1.] 
    p_start = [:k => 2.]

    bprob = BifurcationProblem(rn, u0_guess, p_start, :k; plot_var = :A, u0 = [:A => 5., :B => 3.])
    p_span = (0.1, 6.0)
    opts_br = ContinuationPar(dsmin = 0.0001, dsmax = 0.001, ds = 0.0001, max_steps = 10000, p_min = p_span[1], p_max = p_span[2], n_inversion = 4)
    bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside = true)
    plot(bif_dia, xlabel = "k", ylabel = "A", xlims = (0, 6), ylims=(0,8))

will give the following plot. It looks like ksq is fixed at the original value of k^2 = 2^2 = 4, and does not change when k is varied. (Notice that the equilibrium where [A] = [B] occurs at k = 4).

[isaacsas]

Yeah, this is the kind of issue I was worried about. Thanks for the simple example.

[rveltz]

Hi,

Can you tell me what you would like to see implemented? Did you settle on this?

[isaacsas]

We still need to get a function, rebuild into Catalyst. But somehow we need to be able to recalculate a parameter dependency like the example above when a parameter is varied during continuation. Is there a way to pass BifurcationKit a function that handles updating dependent parameter values when you have a new value arising for the continuation parameter?

[isaacsas]

(I think this is likely an issue in MTK too, but haven't checked.)

[ChrisRackauckas]

No, BifurcationKit just needs to properly account for all equations. This is now in SciMLBse and we'll fix the extension.

[rveltz]

Sorry for being stiff!

I don't understand your problem. Say you want to solve F(x,p) you need to specify the parameters and the initial guess x0. BK function re_make allows you to do that by modifying the NifurcationProbelm so what is the problem?

[ChrisRackauckas]

It's not applying the initializationprobdata which specifies an implicit equation to solve for the parameters and initial conditions in the case where those are implicitly defined through a set of equations. Don't worry though as we are rapidly refactoring this system to be a scimlbase standard so we can more easily apply it to all problems. We are right now generalizing it to all sciml problems and then we will shortly after handle the extensions

[rveltz]

OK

Dont hesitate to ping me for BK