Can no longer simulate coupled `ReactionSystem` `ODESystem`

[TorkelE]

Probably something that has been changed in how stuff should be declared, but this now errors.

@parameters p d a b
@species X(t)
@variables A(t)
eqs = [
    Reaction(p, nothing, [X]),
    Reaction(d, [X], nothing),
    a*A^2 ~ X + b
]
@named coupled_rs = ReactionSystem(eqs, t)
coupled_rs = complete(coupled_rs)
osys = structural_simplify(convert(ODESystem, coupled_rs))

u0 = [X => 0.1, A => 10.0]
tspan = (0.0, 1000.0)
ps = [p => 1.0, d => 0.5, a => 2.0, b => 16.0]

oprob = ODEProblem(osys, u0, tspan, ps) # ┌ Warning: Initialization system is overdetermined.
osol = solve(oprob, Rosenbrock23()) # retcode: InitialFailure
osol.u[end] ≈ [2.0, 3.0] # false


oprob = ODEProblem(osys, [X => 0.1], tspan, ps) # ERROR: Initial condition underdefined. Some are missing from the variable map.

[AayushSabharwal]

osys has 2 unknowns and an algebraic equation 0 ~ b + X(t) - a*(A(t)^2). So it only needs one initial constraint, but u0 has two. Thus the initialization is overdetermined. The initial condition X => 0.1, A => 10.0 also doesn't satisfy this algebraic equation, thus leading to the InitialFailure.

We need a better error for the last case. It's complaining that A doesn't have a guess. If you do guesses(osys)[A] = 1.0 it works and solves successfully.

[TorkelE]

You mean like this?

using Catalyst
t = default_t()
@parameters p d a b
@species X(t)
@variables A(t) [guess = 1.0]
eqs = [
    Reaction(p, nothing, [X]),
    Reaction(d, [X], nothing),
    a*A^2 ~ X + b
]
@named coupled_rs = ReactionSystem(eqs, t)
coupled_rs = complete(coupled_rs)
osys = structural_simplify(convert(ODESystem, coupled_rs))

oprob = ODEProblem(osys, [X => 0.1], tspan, ps; guesses = [A => 1]) # ERROR: UndefVarError: `equations` not defined in `SciMLBase`

(note the new error)

[AayushSabharwal]

Yeah. That error is due to a bad change in SciMLBase. Bumping it down to 2.74 for the time being should work. Also worth noting is that error happens when displaying oprob, so the problem should still be there and solvable.

[TorkelE]

Sounds good, but on a different note, do we even do tests in MTK anymore? The latest updates literally have me getting test failures across 10 different Catalyst test sets. Accidentally introducing bugs happens and is fine, but currently I spend most of my weekends untangling the latest bunch of SciML bugs.

I realise that MTK is a complex system and it is a bit of a learning experience figuring out how it works, but if any even somewhat noticeable update generates multiple bugs then there is a systematic error. Introducing bugs is parts of programming and the way we avoid them is the CI tests. Given that even fairly trivial things slip through, it seems clear to me that the MTK tests (while extensive) just are not at all exhaustive enough (not a direct criticism to you @AayushSabharwal, from what I can tell you are probably better at adding tests than many other).

Like going through the MTK tests and trying to figure out what is missing is a momentous task. But still, it seems to me that a reliable test set is probably one of the most urgent things MTK actually needs. Actually relying on downstream ones like Catalyst's ones would also take us part of the way. However, for some reason it was decided that merging with downstream test failures was the way to go (and now Catalyst is too many bugs behind to be useful). I still think this was a mistake.

[ChrisRackauckas]

Okay, I don't know what is triggering the toxicity here but there is no error left so closing. The only error was in printing, and it was already addressed. Catalyst also doesn't test printing and maybe it should too.

[isaacsas]

Should the error in the non-overdetermined case now be fixed? I'm still getting

julia> oprob = ODEProblem(osys, u0, tspan, ps)
ERROR: Initial condition underdefined. Some are missing from the variable map.
Please provide a default (`u0`), initialization equation, or guess
for the following variables:

Set(Any[A(t)])

when I only pass an initial value for X(t).

[isaacsas]

Example:

using Catalyst
t = default_t()
@parameters p d a b
@species X(t)
@variables A(t)
eqs = [
    Reaction(p, nothing, [X]),
    Reaction(d, [X], nothing),
    a*A^2 ~ X + b
]
@named coupled_rs = ReactionSystem(eqs, t)
coupled_rs = complete(coupled_rs)
osys = structural_simplify(convert(ODESystem, coupled_rs))
u0 = [X => 0.1]
tspan = (0.0, 1000.0)
ps = [p => 1.0, d => 0.5, a => 2.0, b => 16.0]
oprob = ODEProblem(osys, u0, tspan, ps)  # errors about A missing from variable map
oprob = ODEProblem(osys, u0, tspan, ps; guesses = [A => 1]) # works

u0 = [X => 0.1, A => sqrt(16.1/2)]
oprob = ODEProblem(osys, u0, tspan, ps)  # overdetermined warning

How is a user supposed to know that passing a guess is needed? This means they now have to examine the final system and compare what variables it has vs. what they initially constructed to avoid warnings or errors?

[isaacsas]

Is it not possible to automatically put in a guess for variables structural_simplify eliminates that are missing initial values?

[AayushSabharwal]

How is a user supposed to know that passing a guess is needed?

Typically a good idea is to provide a guess for all variables that you don't expect the user to give initial values for. Algebraic variables should basically always have guesses, and in general I'm a proponent for always giving guesses. This is because while typically the user is expected to provide initial values for differential variables and MTK will solve for the rest, which variables are differential may not always be apparent. It's also sometimes possible to calculate initial values for differential variables if an initial value is given for an observed variable whose equation contains the differential variable, so guesses for differential variables can still come in handy. In terms of the modelica spec, defaults are start values with fixed = true, and guesses are start values with fixed = false.

Is it not possible to automatically put in a guess for variables structural_simplify eliminates that are missing initial values?

A is not eliminated from the system - it is an algebraic variable. Observed variables automatically have guesses computed from their equations: e.g. if y ~ 2x + 3 is an observed variable, then the guess for y is taken as 2x + 3 if not explicitly provided. We have a plan to add an enum option to the problem constructor to choose between the following options for handling missing guesses:

• error
• zero
• one
• random value

[ChrisRackauckas]

But even then, the enum is a fallback and not generally good to use, many of the options lead to numerical instability

[ChrisRackauckas]

Should the error in the non-overdetermined case now be fixed? I'm still getting

I don't know what you mean by "fixed" here. This is very clearly and unambiguously a user error. No initial guess was provided, it errors asking for that guess. If you for example were trying to solve an ODE or nonlinear system of 3 variables and only had a u0 of length 2, we would agree that should error right? This is exactly the situation. You can't make a Newton method work with that.

How is a user supposed to know that passing a guess is needed? This means they now have to examine the final system and compare what variables it has vs. what they initially constructed to avoid warnings or errors?

No, that's not true at all. You can always place a guess on anything, and a guess is ignored if it's not necessary. So if you really want to be "safe" you could just put a guess on every algebraic variable (some standard library components do this for that reason). But since many will simplify away, it's not necessary, you only need to do it on the ones that cannot be explicitly solved for. If it's simpler in your head to just always put one for any algebraic variable then that's fine but there's no reason to impose it because in many cases like 90% of variables don't need one and never need one.