Cases of creation of NonlinearSystems have been broken

[TorkelE]

This used to work, but doesn't anymore:

using ModelingToolkit, OrdinaryDiffEq, Test
using ModelingToolkit: t_nounits as t, D_nounits as D

# Defines the model.
@parameters k1 k2 V0
@variables X1(t) X2(t) Y1(t) Y2(t) V(t) = V0 W(t)
@parameters v = V w = W Γ[1:2] = missing [guess = [1.0, 1.0]]
eqs = [
    0 ~ k1 * (Γ[1] - X1) - k2 * X1,
    0 ~ k1 * (Γ[2] - Y1) - k2 * Y1,
    0 ~ 1 - V,
    0 ~ 1 - W,
]
observed = [X2 ~ Γ[1] - X1, Y2 ~ Γ[2] - Y1]
@mtkbuild osys = ODESystem(eqs, t, [X1, X2, Y1, Y2, V, W], [k1, k2, V0, v, w, Γ]; observed)


@mtkbuild nlsys = NonlinearSystem(eqs, [X1, X2, Y1, Y2, V, W], [k1, k2, V0, v, w, Γ]; observed)


# Create the`ODEProblem`s.
u0 = [X1 => 1.0, X2 => 2.0, Y1 => 3.0, Y2 => 4.0, V => 5.0, W => 6.0]
ps = [k1 => 0.1, k2 => 0.2, V0 => 3.0]
nlprob = NonlinearProblem(nlsys, u0, ps);

' This works
oprob = ODEProblem(osys, u0, 10.0, ps);

[AayushSabharwal]

Supporting NonlinearProblem(::ODESystem) is fairly trivial and I have it working locally. However, the MWE expects Gamma to be solved for and doesn't include equations to solve for it. NonlinearProblem only runs parameter initialization, and thus observed equations are not included in initialization.

[TorkelE]

Actually just the creation of the NonlinearProblem from the NonlinearSystem fails here:

nlprob = NonlinearProblem(nlsys, u0, ps);

ERROR: MethodError: Cannot `convert` an object of type Missing to an object of type Real
The function `convert` exists, but no method is defined for this combination of argument types.

[AayushSabharwal]

That's what I mentioned

NonlinearProblem only runs parameter initialization, and thus observed equations are not included in initialization.

It can't solve for Gamma, because there is no parameter equation solving for Gamma. I do see a way to do it by turning X1 into Initial(X1) but will have to draw up a proper plan for making this work.

[TorkelE]

I thought the idea was that the observed equations could be used to solve for Gamma, given that we knew the other parts? I.e. previously we would (and still can) encode default values for Gamma which are used to compute these from the initial conditions. Then, the observed are used to compute X1 and Y2. However, with the recent update the idea is that we no longer provide the default values?

This all works for ODEs, right, so there is some part of the update where we lost NonlinearSystem support?

[AayushSabharwal]

I thought the idea was that the observed equations could be used to solve for Gamma, given that we knew the other parts?

For time-dependent systems, yes. However, running initialization to solve for x doesn't make sense in a NonlinearSystem whose sole purpose is to solve for x. Thus, it runs only parameter initialization. Since it excludes variables, it consequently excludes observed equations from initialization. We didn't "lose" nonlinearsystem support, we fixed a pretty big bug which was causing solves to fail. Currently, we only support a limited form of parameter initialization in nonlinear systems. This will be expanded to handle the case in the issue.

[TorkelE]

Right, so practically, how should we handle this in Catalyst? I.e. what kind of information should we/should not provide? Do correct me if I am wrong, but at least the cases of this that we had in our tests suite in Catalyst used to work without problems.

[AayushSabharwal]

Yeah this would have worked for some time, but it was incredibly broken and thus we hotfixed it out. For time-independent systems, initialization only solves for parameters using parameter dependencies and initialization_eqs. The latter are assumed to be parameter-only equations. We don't validate this but really should.

[TorkelE]

Well, I am glad stuff got fixed at least!

Do you have any idea of when this might be available again?

[AayushSabharwal]

I can't promise a timeline, unfortunately. I'm aware of the issue, but am currently focused on fixing a plethora of CI issues across SciML.

[isaacsas]

So should we make a separate issue for the changes discussed in the last few comments here that mean we can't use conservation law elimination with NonlinearSystems? That is blocker for making a new Catalyst release.

[ChrisRackauckas]

What's the issue?

[TorkelE]

This issue (basically observables and defaults in NonlinearSystems) disables the elimination of conservation laws for NonlinearSystems in Catalyst. While you can still solve NonlinearProblems without conservation laws (although undesired), this elimination is required to e.g. have a non-singular Jacobian. This is essential for quite a few cases, e.g. computing bifurcation diagrams or steady state stabilities. Basically, until this is solved, a decent chunk of Catalyst features relating to steady state and stability analysis is broken.

[isaacsas]

From above:

For time-independent systems, initialization only solves for parameters using parameter dependencies and initialization_eqs. The latter are assumed to be parameter-only equations. We don't validate this but really should.

To solve for the values of conserved constants we need to use initial conditions. Previously, MTK used u0maps passed to NonlinearProblems to solve for them, i.e. if we provided u0maps for all species, even those turned into observables, the value of the parameters representing conserved constants were calculated. This is apparently no longer possible if I understand @AayushSabharwal's posts above.

[ChrisRackauckas]

No the point that is being made is that equation was omitted. I assume that what was meant was:

eqs = [
    0 ~ k1 * (Γ[1] - X1) - k2 * X1,
    0 ~ k1 * (Γ[2] - Y1) - k2 * Y1,
    0 ~ 1 - V,
    0 ~ 1 - W,
]
initialization_eqs = [
  Γ[1] ~ X1 + X2
  Γ[2] ~ Y1 + Y2
]
@mtkbuild osys = ODESystem(eqs, t, [X1, X2, Y1, Y2, V, W], [k1, k2, V0, v, w, Γ]; initialization_eqs)

?

See, this is exactly why the initial values are parameters.

However, the MWE expects Gamma to be solved for and doesn't include equations to solve for it. NonlinearProblem only runs parameter initialization, and thus observed equations are not included in initialization.

Since the initial values are parameters, Γ[1] ~ Initial(X1) + Initial(Y1) and Γ[2] ~ Initial(X2) + Initial(Y2) will be generated as the equations to solve to get the Gammas, and that will be resolved before the ODE/nonlinear system is solved.

[TorkelE]

Just checking, have observed been deprecated for NonlinearSystems? Or is there a reason we cannot use these anymore?

Also, when I try to construct the system this way, indexing seems broken for e.g. X1`:

eqs = [
    0 ~ k1 * (Γ[1] - X1) - k2 * X1,
    0 ~ k1 * (Γ[2] - Y1) - k2 * Y1,
    0 ~ 1 - V,
    0 ~ 1 - W,
]
initialization_eqs = [
  Γ[1] ~ X1 + X2
  Γ[2] ~ Y1 + Y2
]
@mtkbuild nlsys = NonlinearSystem(eqs, [X1, X2, Y1, Y2, V, W], [k1, k2, V0, v, w, Γ]; initialization_eqs)
u0 = [X1 => 1.0, X2 => 2.0, Y1 => 3.0, Y2 => 4.0, V => 5.0, W => 6.0]
ps = [k1 => 0.1, k2 => 0.2, V0 => 3.0]
nlprob = NonlinearProblem(nlsys, u0, ps);
nlprob[X1] # 0.0, should be 1.0.

[AayushSabharwal]

Since the initial values are parameters, Γ[1] ~ Initial(X1) + Initial(Y1) and Γ[2] ~ Initial(X2) + Initial(Y2) will be generated as the equations to solve to get the Gammas, and that will be resolved before the ODE/nonlinear system is solved.

That doesn't happen right now, which is what @TorkelE is referring to. The reason is that it's a little annoying to do since we can't just include all the observed equations. We need the ones that only depend on variables whose value is given as a constant and even then I'm not sure it's always valid. If the user has an observed equation y ~ 2x in a nonlinear system and they provide [x => 1, y => 1] as defaults, do we really want to InitialFailure just because Initial(y) ~ 2Initial(x) isn't satisfied? If we say we only run parameter initialization, why do the values of unknowns fail initialization? I'm not sure where to draw the line here. The argument "just the observed equations that include the parameters we're solving for" doesn't work because what if we have

y ~ 2x
z ~ p - y

and defaults [x => 1.0, z => 1.0]? We need to include both equations.

[AayushSabharwal]

Just checking, have observed been deprecated for NonlinearSystems? Or is there a reason we cannot use these anymore?

observed is very much not deprecated, that would be incredibly breaking and honestly just require special handling in structural simplification to not use observed. Why do you say they can't be used?

[AayushSabharwal]

initialization_eqs = [
  Γ[1] ~ X1 + X2
  Γ[2] ~ Y1 + Y2
]
@mtkbuild nlsys = NonlinearSystem(eqs, [X1, X2, Y1, Y2, V, W], [k1, k2, V0, v, w, Γ]; initialization_eqs)

initialization_eqs containing unknowns is currently UB for NonlinearSystem. As I mentioned earlier, we really should validate this.

[AayushSabharwal]

nlprob[X1] # 0.0, should be 1.0.

No. The system has no unknowns - it is linearly solvable by MTK. So any inital guesses you provide don't hold because MTK analytically solved the system.

[TorkelE]

observed is very much not deprecated, that would be incredibly breaking and honestly just require special handling in structural simplification to not use observed. Why do you say they can't be used?

Sounds good. But in that case, why can I no longer declare observables here as I have previously done? I.e. previously we have just made X2 and Y2 observables of the system. But now we are recommended not to do this?

No. The system has no unknowns - it is linearly solvable by MTK. So any inital guesses you provide don't hold because MTK analytically solved the system.

This is a bit confusing no? I.e. in this case MTK solves my nonlinear problem already when I create it, meaning that the guesses are not stored as values. But in other cases the problem is solved when I make the solve command? I imagine this has to do when you can algebraically do the solution and not?

Finally some other weird stuff. If you do

nlsol = solve(nlprob, NewtonRaphson())

then you got

nlsol[X1] # -0.0 (not 0.0)

like, it is the same thing, so it doesn't matter, but still feels a bit weird.

Also, using this approach I can no longer evaluate X2. I.e.

nlprob[X2] # ERROR: ArgumentError: Symbol X2(t) is not present in the system.

(same for the solution

[isaacsas]

@ChrisRackauckas Just to add on to what @TorkelE is saying. In limited testing I have found that passing initialization_eqs to NonlinearProblem for the conserved constants doesn't give the correct answer (presumably due to the unknowns within the initialziation equations as @AayushSabharwal mentioned above).

[AayushSabharwal]

Sounds good. But in that case, why can I no longer declare observables here as I have previously done? I.e. previously we have just made X2 and Y2 observables of the system. But now we are recommended not to do this?

observed are ignored and overwritten when you call structural_simplify. So if you're identifying observed variables on your own, don't call structural_simplify or include them as normal equations. This behavior has existed for a long time.

This is a bit confusing no? I.e. in this case MTK solves my nonlinear problem already when I create it, meaning that the guesses are not stored as values. But in other cases the problem is solved when I make the solve command? I imagine this has to do when you can algebraically do the solution and not?

What's confusing? You gave a linear system to MTK, MTK solved it for you. Why would you want to provide guesses for something that's linearly solvable? If you don't want MTK to simplify this (for whatever reason) don't call structural_simplify? This is identical to giving ODESystem a purely algebraic system and asking why it doesn't integrate over the time domain - the solution is constant. There's nothing to solve.

Also, using this approach I can no longer evaluate X2. I.e.

This is because you provide the equation for X2 to observed and then call structural_simplify, which overwrites it. So it's not part of the system. Again, we've had this behavior for a long time.

[AayushSabharwal]

In limited testing I have found that passing initialization_eqs to NonlinearProblem for the conserved constants doesn't give the correct answer (presumably due to the unknowns within the initialziation equations as @AayushSabharwal mentioned above).

Yes. If the initialization_eqs are pure parameter equations it will work. If they contain unknowns, it's UB. I agree this is not great API, and have mentioned that we should validate this.

[ChrisRackauckas]

observed is very much not deprecated, that would be incredibly breaking and honestly just require special handling in structural simplification to not use observed. Why do you say they can't be used?

Well two things though. (1) we of course aren't deprecating it, but (2) we purposefully do not document the observed keyword argument. I'm actually surprised @TorkelE found out it exists 😅. What is documented, like on the very first tutorial https://docs.sciml.ai/ModelingToolkit/stable/tutorials/ode_modeling/, is that you simply state your equations. MTK will move things to observed during structural simplify. You almost certainly don't want to be writing this yourself, unless you're specifically writing a custom pass.

That doesn't happen right now, which is what @TorkelE is referring to. The reason is that it's a little annoying to do since we can't just include all the observed equations. We need the ones that only depend on variables whose value is given as a constant and even then I'm not sure it's always valid. If the user has an observed equation y ~ 2x in a nonlinear system and they provide [x => 1, y => 1] as defaults, do we really want to InitialFailure just because Initial(y) ~ 2Initial(x) isn't satisfied? If we say we only run parameter initialization, why do the values of unknowns fail initialization? I'm not sure where to draw the line here. The argument "just the observed equations that include the parameters we're solving for" doesn't work because what if we have

Good point. So then what about letting the user saying it explicitly? Since that's what's wanted here, and we now have the grammar for it:

eqs = [
    0 ~ k1 * (Γ[1] - X1) - k2 * X1,
    0 ~ k1 * (Γ[2] - Y1) - k2 * Y1,
    0 ~ 1 - V,
    0 ~ 1 - W,
]
initialization_eqs = [
  Γ[1] ~ Initial(X1) + Initial(X2)
  Γ[2] ~ Initial(Y1) + Initial(Y2)
]
@mtkbuild nlsys = NonlinearSystem(eqs, [X1, X2, Y1, Y2, V, W], [k1, k2, V0, v, w, Γ]; initialization_eqs)

For the nonlinear system, the Initial(X1) is the guess. So this is saying, use the guesses to define the constant Γ and then solve the system, which IIUC is what's wanted. I haven't checked, but this seems reasonable to support and would be the right way for Catalyst to specify what it wants?

[AayushSabharwal]

X2 doesn't exist in that system. Assuming the user provides Gamma[1] ~ X1 + X2 in eqs as well, this can be made to work even if it is a little redundant.

[AayushSabharwal]

This is a bit of a contrived case though. We're talking about a DSL being able to identify the equations it needs to use to solve for the parameters it wants to solve for. Is this something easily identifiable in general?

[AayushSabharwal]

The reason I ask is that while the proposed syntax works it might not be feasible longer term. If later on we figure out a way to identify the observed equations to use they might conflict with the user-provided equations and lead to a falsely overdetermined system (because equation equality is a superficial check)

[isaacsas]

1. Regarding observed, if it is now being considered purely internal, can we have a user_observed field? We've been using observed for years, and it was something we specifically discussed using long long ago pre-initialization (i.e. this was something I believe we discussed with Yingbo, and updates to MTK/indexing were specifically made based on comments from these discussions). I agree at first it was treated as purely internal, but the understanding we came to back then is that it would be ok for us to use it to store eliminated variables from conservation laws. If that is now changing so be it, though it would be nice if such a change happened in a breaking release, but it would be great to then replace this loss of functionality with something that works equivalently... For now I guess we can just drop structural_simplify use throughout Catalyst and the docs to avoid this issue.
2. I don't think it makes sense to start creating new equations for parameters within the core set of system equations(sys). The whole point of conservation law elimination is to remove equations to lower the system's dimensionality (and to produce a non-singular Jacobian). Moreover, these constants should be fully solvable / calculable before solving the associated nonlinear system since they are given by linear relations among the initial/guess values.

[AayushSabharwal]

We've been using observed for years, and it was something we specifically discussed using long long ago pre-initialization

The issue is not using observed - it's using observed and calling structural_simplify. I don't know how long ago you're talking about, but as far as I know ever since I've worked on MTK we've had this behavior through v8 and v9 so it's not exactly breaking. If you can MWE the behavior you want working on an older version and show it not working now, I can fix that. If you're calling structural_simplify, why do you feel it necessary to explicitly handle conservation laws like this? MTK is capable of doing it anyway.

[AayushSabharwal]

The whole point of conservation law elimination is to remove equations to lower the system's dimensionality

MTK will eliminate those laws.

Moreover, these constants should be fully solvable / calculable before solving the associated nonlinear system since they are given by linear relations among the initial/guess values.

As I illustrated earlier, it's not whether they're solvable, but how to choose the equations to solve them from.

[isaacsas]

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
@variables A(t) B(t)
eqs = [D(A) ~ B - A, D(B) ~ A - B]
@named osys = ODESystem(eqs, t)
osys2 = structural_simplify(osys)

julia> equations(osys2)
2-element Vector{Equation}:
 Differential(t)(A(t)) ~ B(t) - A(t)
 Differential(t)(B(t)) ~ -B(t) + A(t)

julia> full_equations(osys2)
2-element Vector{Equation}:
 Differential(t)(A(t)) ~ B(t) - A(t)
 Differential(t)(B(t)) ~ -B(t) + A(t)

Seems like no conservation law was detected or removed to me.

[TorkelE]

@AayushSabharwal You are right in your initial comment about structural_simplify. I have never really used it, and carelessly assumed that it would act as a normal complete check in my above example. To me, that does seem like an undesirable effect (i.e. silently disabling observables when structural_simplify is used). However, it seems it idea might be to remove it entirely in favour of initialization_eqs (although I am getting confused by initialization_eqs, observables, guesses, default values, and normal algebraic equations, I trust there is a distinction, but I am a bit unsure what is intended to be used where).

[ChrisRackauckas]

if it is now being considered purely internal, can we have a user_observed field? We've been using observed for years, and it was something we specifically discussed using long long ago pre-initialization (i.e. this was something I believe we discussed with Yingbo, and updates to MTK/indexing were specifically made based on comments from these discussions). I agree at first it was treated as purely internal, but the understanding we came to back then is that it would be ok for us to use it to store eliminated variables from conservation laws. If that is now changing so be it, though it would be nice if such a change happened in a breaking release, but it would be great to then replace this loss of functionality with something that works equivalently... For now I guess we can just drop structural_simplify use throughout Catalyst and the docs to avoid this issue.

There is a whole lot to unpack here:

1. Yingbo specifically made comments saying that no one should ever use the observed keyword argument directly. He repeated this all of the time in chats and that's why the tutorials are the way they are: there are comments about removing the direct use of observed.
2. observed as a keyword argument was never documented, at least since v6 IIRC, again because of the implications with structural simplify. This was noted repeatedly in the code reviews to the repo and any reference to it was specifically removed from examples back in the v6 period to avoid potential issues from using internals.
3. This is not new. Look at the first tutorial in v6.0: https://docs.sciml.ai/ModelingToolkit/v6.0/tutorials/ode_modeling/. https://github.com/SciML/ModelingToolkit.jl/releases/tag/v6.0.0 it has been almost half of a decade. This is not new. How to make an un-solved observed quantity is literally the first thing we have shown to every new user since 2021. That literally is not being used figuratively. It is specifically shown in the first tutorial because it's a major feature of MTK, and it is highlighted and discussed.
4. " can we have a user_observed field?" there is observed(sys). This is public, it's documented, but completely unrelated. We are not discussing accessors, we're discussing the direct use in the constructor.
5. Dropping structural_simplify is a really bad idea. It means you will have numerical errors and equations that cannot be solved in many cases. We cannot recommend people to just be incorrect.

Look, I know Catalyst relied on some internals in the v5 period. But let's not use that as justification that we made a breaking change mid v9. I understand where you're coming from, Catalyst mostly learned a very older version of MTK and incrementally moved, and some things didn't move when it should've, but that doesn't mean that things just randomly changed. You're talking about things from 4-5 years ago and 4-5 breaking updates ago. Maybe it doesn't feel like that long ago but the history is all shown there, it really was that long ago. Catalyst now feels like MTK is unstable because these incremental changes means it's using things weirdly/incorrectly. Some things just happened to work in special cases when using internals. But... let's not rely on that. Let's do things right. Let's get Catalyst up to date so it's using the actual public APIs from v8 forwards and everything will feel much more stable! Isn't that what we all want?

Seems like no conservation law was detected or removed to me.

You didn't add the equation. The whole point is that the system doesn't guess things and solves exactly the mathematical description provided.

I have never really used it, and carelessly assumed that it would act as a normal complete check in my above example. To me, that does seem like an undesirable effect (i.e. silently disabling observables when structural_simplify is used).

No observables are "disabled"

X2 doesn't exist in that system. Assuming the user provides Gamma[1] ~ X1 + X2 in eqs as well, this can be made to work even if it is a little redundant.

Okay yes, it's under-conditioned. The complete mathematical description of the system is:

using ModelingToolkit, OrdinaryDiffEq, Test
using ModelingToolkit: t_nounits as t, D_nounits as D

# Defines the model.
@parameters k1 k2 V0
@variables X1(t) X2(t) Y1(t) Y2(t) V(t) = V0 W(t)
@parameters v = V w = W Γ[1:2] = missing [guess = [1.0, 1.0]]
eqs = [
    D(X1) ~ k1 * (Γ[1] - X1) - k2 * X1,
    D(X2)~ k1 * (Γ[2] - Y1) - k2 * Y1,
    D(Y1) ~ 1 - V,
    D(Y2) ~ 1 - W,
    Γ[1] ~ X1 + X2
    Γ[2] ~ Y1 + Y2
]

initialization_eqs = [
  Γ[1] ~ Initial(X1) + Initial(X2)
  Γ[2] ~ Initial(Y1) + Initial(Y2)
]

@mtkbuild osys = ODESystem(eqs, t, [X1, X2, Y1, Y2, V, W], [k1, k2, V0, v, w, Γ]; initialization_eqs)

where then NonlinearSystem(osys) generates the equivalent to:

eqs = [
    D(X1) ~ k1 * (Γ[1] - X1) - k2 * X1,
    D(X2)~ k1 * (Γ[2] - Y1) - k2 * Y1,
    D(Y1) ~ 1 - V,
    D(Y2) ~ 1 - W,
    Γ[1] ~ X1 + X2
    Γ[2] ~ Y1 + Y2
]

initialization_eqs = [
  Γ[1] ~ Initial(X1) + Initial(X2)
  Γ[2] ~ Initial(Y1) + Initial(Y2)
]

@mtkbuild osys = NonlinearSystem(eqs, t, [X1, X2, Y1, Y2, V, W], [k1, k2, V0, v, w, Γ]; initialization_eqs)

Is there anything ambiguous about that?

[AayushSabharwal]

although I am getting confused by initialization_eqs, observables, guesses, default values, and normal algebraic equations, I trust there is a distinction, but I am a bit unsure what is intended to be used where

Observables is something MTK figures out. It's not something you specify to the system directly. All equations that are true about your model during the entire integration should be provided as the equations of the system. MTK is capable of figuring out which ones to make observed. initialization_eqs are equations that are only true at the initial time (t=0 in most cases). They don't hold throughout the integration, but are used when initializing the system. Default values are initial value constraints on the respective variables that can be overridden by the u0 provided to ODEProblem. Guesses are the u0 provided to the nonlinear solver when running initialization if that variable is an unknown for the initialization system. Algebraic equations are (typically nonlinear) equations that MTK can't solve analytically for a set of variables. If MTK were able to solve them analytically, they would be observed equations. Since MTK can't solve them analytically, it relies on DAE solvers to solve them numerically.

However, it seems it idea might be to remove it entirely in favour of initialization_eqs

No, the equations shouldn't be moved directly to initialization_eqs. See the example Chris wrote for how we might want to express this.

Seems like no conservation law was detected or removed to me.

[D(A) ~ B - A, A + B ~ Gamma] will eliminate B as an observed quantity.

Is there anything ambiguous about that?

Consider the case when we just have this one initialization equation in the nonlinear system

G ~ Initial(X) + Initial(Y)

Where X is an unknown in the system, and Y is an observed variable. If the user provides X => 1.0, G => 1.0 as u0, this may or may not satisfy the initialization equation depending on the observed equation used for Y. Do we not include the equation in this case, since Y wasn't explicitly given a constant initial value? Suppose the user provides initial values for all three. Do we error if it doesn't satisfy the equation? What if there are other initialization equations. How do we know it isn't possible to solve for Initial(Y) from them? Do we want to?

[ChrisRackauckas]

this may or may not satisfy the initialization equation depending on the observed equation used for Y.

This is why the observed equations are in the initialization system's construction. So if it's not possible then the system is unsatisfiable.

Suppose the user provides initial values for all three. Do we error if it doesn't satisfy the equation?

Yes.

How do we know it isn't possible to solve for Initial(Y) from them? Do we want to?

If it's possible then we do, yes. It's just a system of equations. We write out the whole list and solve it if possible, or show it's not satisfiable.

[AayushSabharwal]

This is why the observed equations are in the initialization system's construction. So if it's not possible then the system is unsatisfiable.

Yeah but the entire reason for this syntax is we can't include all observed variables in the initialization system of NonlinearSystem. I think we should just use the initialization_eqs, and of all the Initial(x) used in them, if x => Initial(x) isn't in u0map we turn Initial(x) into an unknown.

[ChrisRackauckas]

we can't include all observed variables in the initialization system of NonlinearSystem

Why not? All of those equations have to hold at steady state as well.

[AayushSabharwal]

But not necessarily for the initial guess, right? If we include all observed of the ODE as initial equations in the Nonlinear system we'll very likely end up with an InitialFailure

[TorkelE]

So, my understanding is that you are requesting that we change from directly inputting the eliminated quantities as observables, we instead add the conservation law equations, creating a DAE, and then letting structural_simplify do the step of the elimination and observables creation? I agree that this is nice in the sense that it makes Catalyst and MTK more uniform with how this stuff is handled. However, what I am a bit worried about is how this is going to change for non-ODE systems. I.e. I really do not know how structural simplification generalises to SDEs, and for Jump systems, this is just plain not possible.

Even if this was something that technically was not meant to be accessible, I still think it would be nice to support normal observables that can be used, especially for how these work really effortlessly with stochastic systems (i.e. while I do not know about non-CRN modelling software, this is something that at least in the CRN field all other tools support). This is something that basically works already (with the exception of some cases like structural_simplify that @AayushSabharwal pointed out). While I realise it might not be required for your vision of MTK and might mean some extra work, it really would be useful for us, and we are a part of the SciML system. We could also help with some of the implementation. Ultimately, I think we could figure out some way to adapt things to make pure observables available within the MTK framework (without going through DAEs), and from the Catalyst side we would be really really grateful if we could have a think about this.

[ChrisRackauckas]

I really do not know how structural simplification generalises to SDEs

SDEs are handled the same way. You just add Brownian variables https://docs.sciml.ai/ModelingToolkit/stable/tutorials/stochastic_diffeq/

and for Jump systems

We still need to talk about what to do on Levy processes in general, but I would think we can make that just work the same with a jump variable. In the current reaction setup, we could have structural simplify just removes the equalities from eqs and moves them to observed. It could also tear a resulting system or something if we have a DAE type jump processes etc., but I think we'll be moving here in the very near future. Jump variables might be a cool way to do this in the future.

[ChrisRackauckas]

I still think it would be nice to support normal observables that can be used, especially for how these work really effortlessly with stochastic systems (i.e. while I do not know about non-CRN modelling software, this is something that at least in the CRN field all other tools support)

Letting structural simplification do this is much easier in the sense of making it support different noise types. That's the reason for the Brownian form.

Ultimately, I think we could figure out some way to adapt things to make pure observables available within the MTK framework (without going through DAEs), and from the Catalyst side we would be really really grateful if we could have a think about this.

What would be the advantage?

[TorkelE]

Is there actually any theory on what happens with the structural simplification algorithms when applied to SDEs? It also doesn't really work with Catalyst, as those SDEs depend on creating a System from the equations (which may or may not contain noise). However, for Catalyst the entire point is that a ReactionSystem can generate all system types (making convert(System, rs::ReactionSystem) ambiguous). Although it should be straightforward to simply allow SDESystems to be created directly from equations containing Brownian motions/Wiener processes.

Also, is then plans to just disable conservation laws for JumpSystems until further notice? I think if that is the case, we would rather just stay in our current MTK version where conservation laws work, and then once there is a tearing implementation for Jumps we will look at transfering Catalyst to the latest MTK version.

[ChrisRackauckas]

Is there actually any theory on what happens with the structural simplification algorithms when applied to SDEs?

Yes. https://www.sciencedirect.com/science/article/pii/S0377042703010239 justifies the mass matrix index 1 DAE solvers in StochasticDiffEq and https://arxiv.org/abs/2102.00605 justifies the index reduction algorithm. Note that we specifically avoid what's required to build the counterexample by not allowing noise terms in the algebraic equations.

Also, is then plans to just disable conservation laws for JumpSystems until further notice? I think if that is the case, we would rather just stay in our current MTK version where conservation laws work, and then once there is a tearing implementation for Jumps we will look at transfering Catalyst to the latest MTK version.

Uhh what? There is at least no regression from previous versions.... where is this even coming from? What is this even referencing?

[isaacsas]

Just to clarify, we don't allow conserved species removal for jump models since it would probably be less computationally efficient (doing it for the chemical master equation does make sense).

[TorkelE]

Sam is right, I got int confussed, I meant observables generally.

Yes. https://www.sciencedirect.com/science/article/pii/S0377042703010239 justifies the mass matrix index 1 DAE solvers in StochasticDiffEq and https://arxiv.org/abs/2102.00605 justifies the index reduction algorithm. Note that we specifically avoid what's required to build the counterexample by not allowing noise terms in the algebraic equations.

It this is basically just reduction of SADEs? I.e. pure SDEs won't be touched?

Uhh what? There is at least no regression from previous versions.... where is this even coming from? What is this even referencing?

Maybe I overreacted, if observables still is going to work like they are currently I agree there is no major problem. It would prefer to handle these uniformly at least across Catalyst. I.e. until structural_simplify is reliably supported for Jumps (and doesn't affect performance), I'd like to support observables like we currently do (i.e. providing them directly, rather than creating algebraic equations which structural_simplify then handles).

[ChrisRackauckas]

It this is basically just reduction of SADEs? I.e. pure SDEs won't be touched?

Yes. x ~ y will be reduced, but if there's a noise equation in there it won't touch or try to reduce it. Likely if you have an algebraic equation with noise it's just unknown what to do, with some of those equations unsolvable like discussed in that paper. So that case just isn't handled by structural simplify at all. That also means you cannot have an observed with extra noise, because if all observed come from structural simplification then it won't move a Brownian, which also fixes the fact that it's not clear how that would be defined either (or how that would work, whether it's always well-defined, etc.)

Maybe I overreacted, if observables still is going to work like they are currently I agree there is no major problem. It would prefer to handle these uniformly at least across Catalyst. I.e. until structural_simplify is reliably supported for Jumps (and doesn't affect performance), I'd like to support observables like we currently do (i.e. providing them directly, rather than creating algebraic equations which structural_simplify then handles).

And just to be clear here, as Sam says, jump models don't use structural reduction because it's not necessarily helpful. With ODEs, enforcing x^2 + y^2 = 1 by defining y = sqrt(1 - x^2) can help, emphasis on can because of dummy derivative singularities etc. see the whole JuliaCon workshop on that. But with jump processes, you have discrete updates, +1, -1, and because of that by design you have conserved quantities in the solver. So you don't get a numerical accuracy improvement by removing them. Though in theory in some cases you could subset your model by eliminating some variables, that's just a minor performance thing.

[TorkelE]

Right, that makes sense, got it. I think that is one reason I think it would be useful to keep the option of general observables open, although I realise it might not been as trivial to integrate with the remaining system as I had hoped.