try to bump BifurcationKit

[isaacsas]

No description provided.

[isaacsas]

@vyudu could you try to get this extension working with BifurcationKit 0.4? You can see the other open PR about splitting, and also what was recently merged into MTK to get it working there for help.

[ChrisRackauckas]

ModelingToolkit.jl now has a BirfurcationKit extension. Is there anything the Catalyst one does that is special that makes it so that it's better than targeting the standard MTK one, or is it legacy?

[isaacsas]

Didn’t @TorkelE write the original MTK one? Ours has always been a thin wrapper around the MTK one that handles a few things needed when converting from a ReactionSystem, like elimination of conserved equations.

[TorkelE]

Yeah, I wrote the original MTK extension as well. The Catalyst one basically just handles the elimination of conservation laws to make non-singular Jacobians.

[ChrisRackauckas]

The Catalyst one basically just handles the elimination of conservation laws to make non-singular Jacobians.

Why would it complete instead of SS?

[TorkelE]

sorry, I don't think I understand what you mean

[ChrisRackauckas]

https://github.com/SciML/Catalyst.jl/blob/master/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl#L25-L27

If you just structurally simplified it would eliminate the conservation laws and it wouldn't need something special?

[TorkelE]

So conservation law elimination never generates any DAEs, it generates a pure ODE with the conservation laws removed, hence structural simplify is never required.

[ChrisRackauckas]

But you're building a nonlinear system, you can eliminate a lot more?

[ChrisRackauckas]

Also that wouldn't SCC decompose?

[TorkelE]

sorry, I don't think I follow? Generally it is not possible to generate more than the number of conservation laws? I have tried to generate these using only structural simplify, but it doesn't seem to be able to find them even in the simplest of cases right now.

[ChrisRackauckas]

It's a nonlinear system so you have a larger set of algebraic relations: every ODE relation turns into an algebraic relation. Because of this, there's many more equations to simplify than the conservation laws.

[TorkelE]

One could look at it again, but we still need the conservation laws as SS cannot handle those. We could then test adding SS on top, but that would just be another layer?

[ChrisRackauckas]

but we still need the conservation laws as SS cannot handle those

Why can't it handle those?

[TorkelE]

I don't know. E.g. this is basically the simplest case:

using Catalyst
rs = @reaction_network begin
    (k1,k2), X1 <--> X2
end
osys = convert(ODESystem, rs)
osys = structural_simplify(osys)
equations(osys)

however, structural_simplify generates a 2-equation ODE, while Catalyst's conservation law elimination can create a single-equation ODE + and observable.

[ChrisRackauckas]

But for BK you convert to a nonlinear system.

[TorkelE]

Hmm, true, it seems that SS works differently for the different system types.

[ChrisRackauckas]

For example with

using Catalyst
rs = @reaction_network begin
    (1,1), X1 <--> X2
end
osys = convert(NonlinearSystem, rs)
osys = structural_simplify(osys)
equations(osys)

it would eliminate more than the conservation law, it would reduce it all the way down to an analytical solution.

@AayushSabharwal do you know why this case with the k1 k2 wouldn't eliminate even with allow_parameters =true ?

1-element Vector{Equation}:
 0 ~ k1*X1(t) - k2*X2(t)

a linear parameter equation should eliminate.

[ChrisRackauckas]

Hmm, true, it seems that SS works differently for the different system types.

It always does elimination on the algebraic system set. For DAEs that's just the subset of equations that are algebraic, so it would eliminate the conservation laws. On a nonlinear system, like solving for a steady state, that would be all of the equations, so you'd get a much larger benefit. Especially with the coming SCC codegen, you'd get even more benefits from how it would split the equation.

[TorkelE]

What is the SCC codegen? Is it something new?

[ChrisRackauckas]

SciML/NonlinearSolve.jl#388

it's our big new thing that's going to accelerate most of our solvers by an order of magnitude. It will take a few months to roll it out, but we're getting close to the end for it.

[ChrisRackauckas]

Note that stiff ODEs will benefit from it too, but it's a very detailed discussion as to how.

[isaacsas]

it would eliminate more than the conservation law, it would reduce it all the way down to an analytical solution.

But without specifying a constant for the conservation law, or fixing the value of one of the unknowns appearing in it, there are an infinite set of solutions. Does the SCC process introduce constants to represent the analytical solutions in terms of, and provide an interface for specifying their value? i.e for the A <--> B reaction, without choosing a value for one of A, B, or the conserved constant, there is a infinite solution space at SS.

The process here is setup (when working) to introduce such a constant into the reduced set of equations, thereby making it a parameter of the generated system (with a corresponding default in terms of the unknowns it can be defined from).

[isaacsas]

I think for now we should keep doing conservation law elimination on our end when generating steady state systems. But we should add calling structural_simplify or making it an optional kwarg whereever we don't currently do that.

More generally, in problem conversion (i.e. in ODEProblem and NonlinearProblem), I don't recall why we set structural_simplify = false as the default. I'm pretty sure we discussed this and @TorkelE wanted to turn it on by default, so I must have had some reason I said we shouldn't. But if no one can remind me of what it is I'd be fine with changing the default to structural_simplify = true in problem generation.

[isaacsas]

And on this specific example, if we eliminate the conservation law in Catalyst, structural_simplify does then fully solve the system.

[isaacsas]

I'd also be fine with making conservation law removal a kwarg (it is for SciML problems, just not for extension problems I think).

[AayushSabharwal]

For example with

using Catalyst
rs = @reaction_network begin
    (1,1), X1 <--> X2
end
osys = convert(NonlinearSystem, rs)
osys = structural_simplify(osys)
equations(osys)

it would eliminate more than the conservation law, it would reduce it all the way down to an analytical solution.

@AayushSabharwal do you know why this case with the k1 k2 wouldn't eliminate even with allow_parameters =true ?

The MWE above doesn't have any parameters, and reduces to zero equations. However, it is weird in that both equations are identical:

julia> rs = @reaction_network begin
           (1,1), X1 <--> X2
       end
julia> osys = convert(NonlinearSystem, rs)
Model ##ReactionSystem#360 with 2 equations
# ...

julia> osys2 = structural_simplify(osys)
Model ##ReactionSystem#360 with 0 equations
# ...

julia> equations(osys)
2-element Vector{Equation}:
 0 ~ -X1(t) + X2(t)
 0 ~ X1(t) - X2(t)

julia> observed(osys2)
2-element Vector{Equation}:
 X2(t) ~ -0.0
 X1(t) ~ X2(t)

The modified version does in fact not simplify fully:

julia> rs = @reaction_network begin
           (k1,k2), X1 <--> X2
       end
julia> osys = convert(NonlinearSystem, rs)
Model ##ReactionSystem#367 with 2 equations
# ...

julia> osys2 = structural_simplify(osys)
Model ##ReactionSystem#367 with 1 equations
# ...

julia> equations(osys)
2-element Vector{Equation}:
 0 ~ -k1*X1(t) + k2*X2(t)
 0 ~ k1*X1(t) - k2*X2(t)

julia> observed(osys2)
1-element Vector{Equation}:
 X1(t) ~ (k2*X2(t)) / k1

julia> equations(osys2)
1-element Vector{Equation}:
 0 ~ k1*X1(t) - k2*X2(t)

The identical equations seem to be due to the fact that convert(ODESystem, rs) has the following equations:

2-element Vector{Equation}:
 Differential(t)(X1(t)) ~ -k1*X1(t) + k2*X2(t)
 Differential(t)(X2(t)) ~ k1*X1(t) - k2*X2(t)

And at steady state both rates are zero. As such, the first case without parameters the solution seems incorrect since the system is underdetermined. X2 can take any value and X1 ~ X2 whereas the simplification would suggest X1 = X2 = 0 is the only solution.

[isaacsas]

julia> observed(osys2)
2-element Vector{Equation}:
 X2(t) ~ -0.0
 X1(t) ~ X2(t)

I agree with @AayushSabharwal , this seems like a silent correctness issue to set X2 = 0 instead of giving an error that the system is not uniquely solvable or introducing a constant as Catalyst does.

[AayushSabharwal]

So does catalyst turn this into X1 ~ c, X2 ~ X1?

[isaacsas]

Along those lines I would have thought a more general solution would use a user's u0 guess in building a NonlinearProblem to choose the solution manifold:

julia> nprob = NonlinearProblem(osys2, [osys2.X1 => 1.0, osys2.X2 => 1.0])
NonlinearProblem with uType Nothing. In-place: true
u0: nothing

julia> sol = solve(nprob)
retcode: Success
u: Float64[]

julia> sol[osys2.X1]
-0.0

But it seems to ignore it here.

Catalyst, gives:

julia> osys3 = convert(ODESystem, rs; remove_conserved = true)
Model ##ReactionSystem#235 with 1 equations
Unknowns (1):
  X1(t)
Parameters (1):
  Γ[1] [defaults to X1(t) + X2(t)]

julia> equations(osys3)
1-element Vector{Equation}:
 Differential(t)(X1(t)) ~ -2X1(t) + Γ[1]

julia> observed(osys3)
1-element Vector{Equation}:
 X2(t) ~ -X1(t) + Γ[1]

which eliminates X2 but introduces a new constant defined via the initial conditions of X1 and X2.

[isaacsas]

Ideally then a user could either pass X1(0) and X2(0) to determine gamma (via u0) or could directly set gamma if they happen to know its value.

[isaacsas]

Ideally then a user could either pass X1(0) and X2(0) to determine gamma (via u0) or could directly set gamma if they happen to know its value.

(Which does work currently outside of remake, I was just commenting on the intended workflow we had in mind when providing this functionality.)

[AayushSabharwal]

You're comparing two different things?

The broken osys2 will not solve correctly. Infact, structural_simplify should have thrown an error before it even got to the solve phase, since it requires users to explicitly state that their systems are underdetermined by passing fully_determined = false.

julia> osys3 = convert(ODESystem, rs; remove_conserved = true)

Here, you're explicitly adding an equation X1 + X2 ~ Gamma[1]. MTK has no reason to do so, since it doesn't have the domain knowledge that Catalyst has.

[AayushSabharwal]

Ah, I see, you're saying that since D(X1) = -D(X2) then X1 + X2 ~ Gamma[1]. How does catalyst detect this case?

[isaacsas]

We exploit the very specific structure of ODEs arising from reactions (i.e. all first order) to find linear combinations of the derivatives that are zero. This is the code that does it, but it actually uses the integer nullspace code in MTK (Yingbo was nice enough to modify it so we could use it here too):

Catalyst.jl/src/network_analysis.jl

Line 700 in 6de14a0

function conservationlaws(rs::ReactionSystem)

We give a warning about mixing non-reaction ODEs into the model when doing such detection, since we only check the subset of ODEs arising from reactions.

[isaacsas]

If you look at https://docs.sciml.ai/Catalyst/stable/model_creation/network_analysis/ you can see the standard form of reaction models, $\frac{dx}{dt} = N v(x(t))$. Looking for conservation laws then can be done by examining the left null space of the constant matrix $N$.

[AayushSabharwal]

I see, thanks for the clarification. So yeah, I guess MTK can't easily do this in general. Given the new equation, though, it will simplify it out.

[ChrisRackauckas]

Yes, Catalyst needs to derive the conservation laws, but the simplification part should be automatic.

[isaacsas]

Yes, Catalyst needs to derive the conservation laws, but the simplification part should be automatic.

Hopefully the simplification would also be correct, and not introduce new assumptions like X2 = 0. (Or at least not silently make such an assumption.)

[ChrisRackauckas]

It can't make any assumptions that aren't in the equations or that's a bug.