Updating parameter/variable values does not update values of dependant defaults

[TorkelE]

I thought this one was reported, but didn't actually find it anywhere. This is both a silent error (not just a crash), and something that we have Catalyst users that report errors with, so would be good to have a fix asap so that we can finally transition Catalyst to the MTK v9.

using ModelingToolkit
@variables t X(t)
@parameters p d = X
D = Differential(t)

eqs = [D(X) ~ p - d*X]
@mtkbuild osys = ODESystem(eqs, t)

u0 = [X => 1.0,]
ps = [p => 1.0]
oprob = ODEProblem(osys, u0, 1.0, ps)

oprob.ps[:d] # 1.0
oprob[:X] # 1.0
oprob[:X] = 2.0
oprob.ps[:d] # 1.0, should be 2.0

[ChrisRackauckas]

It wasn't listed as a dependent parameter, just has a default value

[TorkelE]

I am not sure I follow? Here, d's value is set to, by default, be whatever value Xhas. If the value of X is updated so should the value ofd`. I think this used to work before, but seems to be broken in the MTK v9 transition

[SebastianM-C]

The thing is that the default values for parameters encode that they should only initially be that and not that it's a relationship that needs to be enforced. To enforce the relationship one needs to use the parameter dependencies keyword argument as @ChrisRackauckas pointed out.

This might be a bit confusing since for variables the defaults refer to an actual constraint and remake seems to treat defaults as enforcing equations, but setp does the correct thing (IMO) and does not update parameters based on defaults.

[TorkelE]

Up until now the default value are the value of species/parameters at the beginning of the simulation (after which they can change). The the above example one would expect the initial value of d at the beginning to be the same as X, whatever X is. However, I think we are both in agreement that d should not change throughout a simulation as X changes, I agree that that is not what default are. Generally, there should be no difference in the behaviours of remake and setp/setu.

Parameter dependencies is, if I understand it right, a different thing. However, there doesn't seem to be any documentation of it?

[hersle]

I would expect

prob = remake(oprob; u0 = [X => 2.0], use_defaults = true)
@assert prob.ps[d] == 2.0

to do what the original example tries to do. I'm not on my computer to check if it works, though. IIRC, the behavior of remake here could depend on whether p is passed to it.

[AayushSabharwal]

A clarification on terminology: parameter dependencies and dependent defaults are two different things.

Dependent defaults refers to the default value of a variable depending on another variable (this can be an unknown or a parameter, or even an expression involving a mixture of the two). Note that giving X a default value of p, and updating p won't update X in the ODEProblem.

Parameter dependencies can be thought of as observed for parameters. The relationships defined here are maintained throughout the simulation. For example, if [p1 => 2p2] is a parameter dependency and p2 is modified in a callback, p1 will also be updated accordingly. Right now, dependent parameters (p1) are stored in the MTKParameters object. However, there is a future plan to treat these the same way we do observed: specify p1 ~ 2p2 as an equation, and compute p1 on the fly as necessary.

However, I can see the inconsistency here, where defaults for unknowns are considered when creating an initialization system for the ODEProblem, but ones for parameters aren't. This ties directly in to #2665, which requires allowing parameters in the ODESystem to be unknowns in the initialization system. Note that while the above issue is something I'm working on, this will mean that oprob[X] = 2.0 won't update oprob.ps[d] directly, this update would instead be performed during initialization of the integrator (init(oprob, Tsit5())).

[AayushSabharwal]

to do what the original example tries to do. I'm not on my computer to check if it works, though. IIRC, the behavior of remake here could depend on whether p is passed to it.

Correct, you need to pass p = Dict() here to update the parameters.

[TorkelE]

I still think there is an issue here that can and should be solved before the creation of the integrator and simulation system initialization. I.e.

@variables t X(t)
@parameters p d = X

...

oprob.ps[:d] # 1.0
oprob[:X] # 1.0
oprob[:X] = 2.0
oprob.ps[:d] # This should yield `2.0`

Everything else is going to be greatly confusing to users, and seems like asking for more problems down the line (and also cause problems for other packages, like DynamicalSystems, which have some functionality to use ODEProblems). I.e. whatever happens to set the values when ODEProblem is called should probably happen when values are updated. Then it might be useful to have some kind of flag/metadata in ODEProblem to inform setp that there are some default dependencies going on.

[AayushSabharwal]

The thing is in general oprob can't represent the true initialization. It already doesn't do it for unknowns, where the variables that are solved for in the initialization problem are just set to 0 in the ODEProblem. Parameters would probably proceed similarly. d = X is a simple default initialization, but in general this may not be invertible. What if the user changes d? The initialization problem would solve for X, but this wouldn't be reflected at the ODEProblem level.

[TorkelE]

If the user changes d then that is its new value (since it now has a value, the default is no longer used). Initialization problems never really were a thing until the latest update, so I am not really sure why their introduction should pose a problem to default values?

[AayushSabharwal]

Consider an example where f(x, y) = 0 is an algebraic equation and both x and y are unknowns. If a user specifies a guess for x, and passes y => 42.0 to the ODEProblem constructor, x will be initialized to 0 in the ODEProblem. This does not mean the initial value of x is 0. It just means that when the integrator is initialized, it will solve a nonlinear system to calculate x, using the algebraic equation. Modifying x in the ODEProblem won't be reflected in the integrator.

If I now specify a default for x ~ 2y + z, this is now an additional equation in the nonlinear problem which must be satisfied. In short, defaults affect initialization.

#2665 raises a valid concern where parameters can also be unknowns in the initialization problem. Thus, they will effectively get the same treatment. If @parameters d = f(X) is a default value, it will be an equation in the initialization system, and will always be satisfied. Thus while updating X will change d (and this may be trivial) but updating d will also update X when the integrator is initialized (and this is not trivial to do at the ODEProblem level).

[TorkelE]

Right, so I agree that defaults will have effects down the line. However, this was a thing before we started working with initialization and focus on DAEs, I don't see what the introduction of this means that we should throw away old features that people have been utilising.

Generally, it is not that defaults permits you to suddenly do something entirely different, e.g.

using ModelingToolkit
@variables t X(t)
@parameters p d = X
D = Differential(t)

eqs = [D(X) ~ p - d*X]
@mtkbuild osys = ODESystem(eqs, t)

u0 = [X => 1.0,]
ps = [p => 1.0]
oprob = ODEProblem(osys, u0, 1.0, ps)
oprob[:X] = 2.0

should just be the same as setting X => 2.0 in the first place:

using ModelingToolkit
@variables t X(t)
@parameters p d = X
D = Differential(t)

eqs = [D(X) ~ p - d*X]
@mtkbuild osys = ODESystem(eqs, t)

u0 = [X => 2.0,]
ps = [p => 1.0]
oprob = ODEProblem(osys, u0, 1.0, ps)

which is hardly controversial.

[AayushSabharwal]

I don't see what the introduction of this means that we should throw away old features that people have been utilising.

Did this work before? Updating X and d automatically changing?