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 https://github.com/SciML/ModelingToolkit.jl/issues/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.

https://github.com/SciML/ModelingToolkit.jl/issues/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?

[TorkelE]

Well, it was the intended behaviour, there were bugs going around before v9 as well. In this case, we only had the bug being reported this spring.

[isaacsas]

@ChrisRackauckas can we get a decision on whether this should be supported? This is currently a blocker on making a Catalyst 14 release with MTK9 support and there seems to be a lot of back and forth in this thread.

Right now, and in MTK8 and earlier, it was supported to declare a parameter as a function of unknowns/variables via defaults and have this implicitly mean the parameter was defined in terms of the initial values of the unknowns. It is true that it appears remake never worked with this, silently being bugged and not updating such parameters in the past, but otherwise this feature worked and was used by Catalyst. Currently, such parameter definitions are still supported and work except in remake type uses as described above. But it is inconsistent to support defining parameters in such a way except for updating. If this is not meant to be allowed then it should really explicitly error during system construction when such a parameter is detected.

If defining, and consistently updating, parameters that are given by functions of the initial values of unknowns is no longer supported we need to deprecate conservation law elimination in Catalyst 14 and let Catalyst users know it will no longer be supported. Otherwise we are stuck on MTK8 until this is resolved.

[isaacsas]

Just to add, if this is something that is supposed to work with parameter_dependencies we are happy to switch to using that to define conservation constants instead of using defaults, but it isn't clear to me if this is considered supported in that context either from the preceding discussion.

[isaacsas]

And modifying @TorkelE's example to use parameter_dependencies seems to error suggesting there is no support for using initial values of unknowns in defining a dependent parameter. Or is there some other way this should be done now?

@mtkbuild osys = ODESystem(eqs, t; parameter_dependencies=Dict(d => 2*X))
oprob = ODEProblem(osys, u0, 1.0, ps)

gives

ERROR: UndefVarError: `X` not defined in `ModelingToolkit`
Suggestion: check for spelling errors or missing imports.
Stacktrace:
  [1] macro expansion
    @ ~/.julia/packages/SymbolicUtils/qyMYa/src/code.jl:418 [inlined]
  [2] macro expansion
    @ ~/.julia/packages/Symbolics/Eas9m/src/build_function.jl:546 [inlined]
  [3] macro expansion
    @ ~/.julia/packages/SymbolicUtils/qyMYa/src/code.jl:375 [inlined]
  [4] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:163 [inlined]
  [5] macro expansion
    @ ./none:0 [inlined]
  [6] generated_callfunc
    @ ./none:0 [inlined]
  [7] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{…})(::RecursiveArrayTools.ArrayPartition{…}, ::Vector{…}, ::Vector{…})
    @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:150
  [8] ModelingToolkit.MTKParameters(sys::ODESystem, p::Dict{…}, u0::Vector{…}; tofloat::Bool, use_union::Bool)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/cAhZr/src/systems/parameter_buffer.jl:152
  [9] ModelingToolkit.MTKParameters(sys::ODESystem, p::Dict{SymbolicUtils.BasicSymbolic{…}, Float64}, u0::Vector{Pair{…}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/cAhZr/src/systems/parameter_buffer.jl:13
 [10] process_DEProblem(constructor::Type, sys::ODESystem, u0map::Vector{…}, parammap::Vector{…}; implicit_dae::Bool, du0map::Nothing, version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, use_union::Bool, tofloat::Bool, symbolic_u0::Bool, u0_constructor::typeof(identity), guesses::Dict{…}, t::Float64, warn_initialize_determined::Bool, build_initializeprob::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/cAhZr/src/systems/diffeqs/abstractodesystem.jl:943
 [11] process_DEProblem
    @ ~/.julia/packages/ModelingToolkit/cAhZr/src/systems/diffeqs/abstractodesystem.jl:835 [inlined]
 [12] (ODEProblem{…})(sys::ODESystem, u0map::Vector{…}, tspan::Float64, parammap::Vector{…}; callback::Nothing, check_length::Bool, warn_initialize_determined::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/cAhZr/src/systems/diffeqs/abstractodesystem.jl:1086
 [13] (ODEProblem{…})(sys::ODESystem, u0map::Vector{…}, tspan::Float64, parammap::Vector{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/cAhZr/src/systems/diffeqs/abstractodesystem.jl:1076
 [14] (ODEProblem{true})(::ODESystem, ::Vector{Pair{Num, Float64}}, ::Vararg{Any}; kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/cAhZr/src/systems/diffeqs/abstractodesystem.jl:1063
 [15] (ODEProblem{true})(::ODESystem, ::Vector{Pair{Num, Float64}}, ::Vararg{Any})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/cAhZr/src/systems/diffeqs/abstractodesystem.jl:1062
 [16] ODEProblem(::ODESystem, ::Vector{Pair{Num, Float64}}, ::Vararg{Any}; kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/cAhZr/src/systems/diffeqs/abstractodesystem.jl:1052
 [17] ODEProblem(::ODESystem, ::Vector{Pair{Num, Float64}}, ::Vararg{Any})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/cAhZr/src/systems/diffeqs/abstractodesystem.jl:1051
 [18] top-level scope
    @ REPL[32]:1
Some type information was truncated. Use `show(err)` to see complete types.

[ChrisRackauckas]

It should be parameter_dependencies yes. But it only handles parameters right now.

[AayushSabharwal]

It should be parameter_dependencies yes. But it only handles parameters right now.

No, it shouldn't be parameter dependencies

The way to think about this is, parameter_dependencies are equations that are part of the ODESystem involving only parameters (and the lhs is a single parameter). They are like observed equations. If I have @variables x(t) y(t) and have an equation y ~ 2 * x, y is not part of the state but is calculated on the fly using the observed equation. Similarly, if p1 => 2 * p2 is a parameter dependency, p1 is not part of the parameters but is calculated on the fly using observed equations.

Defaults such as @parameters d = X are equations used only during initialization (not yet merged, but implemented in a few PRs). This is identical behavior to @variables X(t) = Y(t). If an explicit value for d is provided to ODEProblem, it will keep that value. If no value is provided, d is an unknown in the initialization system, d ~ X is an equation, and d will be solved for during the initialization stage, which happens when the integrator is created (init(prob, Alg()) is called).

[TorkelE]

But for this case we should be good, right? Let me explain our application

Conserved quantities in chemical reaction networks

We have a reaction system with a single thing (X) which can exist in two states (X1 and X2). Basically we have these two reactions:

k1, X1 --> X2
k2, X2 --> X1

This creates an ODE with two equations

dX1/dt = -k1*X1 + k2*X2
dX2/dt = k1*X1 - k2*X2

However, we then make the observation that the quantity X1 + X2 is conserved, i.e. does not change throughout the simulation. If we then create a parameter, Γ, to represent this conserved quantity we get an equation Γ = X1 + X2. From this we can eliminate one variable (and thus ODE) from our system, creating:

dX1/dt = -k1*X1 + k2*(Γ - X1)

and representing the species X2 as an observable

X2 = Γ - X1

(this is a bit similar to what structural_simplify does, however, chemical reaction networks have specialised structures which we can utilise to do this more efficiently)

Computing conserved quantities using defaults

Now the only question is, how do we compute Γ? Well, it is a constant Γ = X1 + X2, and can hence be computed from the system's initial condition. I.e. under the hood, we declare the parameter Γ

@parameters Γ = X1 + X2

Since X1 + X2 is changed throughout the simulation, we can _compute it from the initial condition. I.e., we do not want to have this a dynamic equation throughout the system. We just want to make sure that, when the simulation is initialised, the value of Γ is computed from the initial condition.

The problem

So the problem is that, right now, the value of parameters that have initial conditions are set when a ODEProblem is declared, but never updated when setu or remake is called. That means that they have the wrong value when you create an integrator from the ODEProblem. Basically, what is needed, is that when setu or remake is called, the value of unknowns/parameters that depend on the updated thing must also be updated.

[TorkelE]

Basically, the problem we see here is the same as we see when we have a parametric initial condition:

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

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

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

oprob.ps[X0] # 2.0
oprob[X] # 2.0
oprob.ps[X0] = 3.0
oprob[X] # 2.0, but should be 3.0

While the initial example I gave her might be a bit non-intuitive, I think it should be clear why the behaviour in this example is un-desired.

[hersle]

Maybe I am missing something, but is there a common need for parameter_dependencies? Is not its use case already covered entirely by defaults? A system consists of a set of variables and parameters, so there are 2 x 2 = 4 possible "types" of "defaults assignments" that could be made between them during initialization. Let us try to cover all of them, both with creating problems with ODEProblem and updating them with remake:

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

@testset begin "all 2x2 combinations of variable <-> parameter defaults"

vars = @variables x1(t) x2(t)
pars = @parameters p3 p4
eqs = [D(x1) ~ 0, D(x2) ~ 0]

@named sys1 = ODESystem(eqs, t, vars, pars; defaults = [x1 => x2*1/2]) # var(var) default
@named sys2 = ODESystem(eqs, t, vars, pars; defaults = [x2 => p3*2/3]) # var(par) default
@named sys3 = ODESystem(eqs, t, vars, pars; defaults = [p3 => x2*3/2]) # par(var) default
@named sys4 = ODESystem(eqs, t, vars, pars; defaults = [p4 => p3*4/3]) # par(par) default
sys1, sys3, sys2, sys4 = complete.([sys1, sys3, sys2, sys4]) # complete all systems

tspan = (0.0, 1.0)

# first, create and initialize all problems
prob1 = ODEProblem(sys1, [           x2 => 2.0], tspan, [p3 => 3.0, p4 => 4.0]) # calculate missing x1 from defaults, works
prob2 = ODEProblem(sys2, [x1 => 1.0,          ], tspan, [p3 => 3.0, p4 => 4.0]) # calculate missing x2 from defaults, works
prob3 = ODEProblem(sys3, [x1 => 1.0, x2 => 2.0], tspan, [           p4 => 4.0]) # calculate missing p3 from defaults, works
prob4 = ODEProblem(sys4, [x1 => 1.0, x2 => 2.0], tspan, [p3 => 3.0           ]) # calculate missing p4 from defaults, works

for prob in [prob1, prob2, prob3, prob4]
    @test prob[x1] == 1.0 && prob[x2] == 2.0 && prob.ps[p3] == 3.0 && prob.ps[p4] == 4.0
end

# second, update the variables and parameters by switching their signs
prob1 = remake(prob1, u0 = [            x2 => -2.0], p = [p3 => -3.0, p4 => -4.0], use_defaults = true) # works
prob2 = remake(prob2, u0 = [x1 => -1.0,           ], p = [p3 => -3.0, p4 => -4.0], use_defaults = true) # fails
prob3 = remake(prob3, u0 = [x1 => -1.0, x2 => -2.0], p = [            p4 => -4.0], use_defaults = true) # works
prob4 = remake(prob4, u0 = [x1 => -1.0, x2 => -2.0], p = [p3 => -3.0            ], use_defaults = true) # works

for prob in [prob1, prob2, prob3, prob4]
    @test prob[x1] == -1.0 && prob[x2] == -2.0 && prob.ps[p3] == -3.0 && prob.ps[p4] == -4.0
end

end

This all works for me on master, except the line with remake(prob2, ...), as reported in #2715. The use of parameter_dependencies sounds fully equivalent to defaults in case 4 to me, and I only find it confusing to introduce a special term for this.

EDIT: @TorkelE's example just above is also equivalent to the case remake(prob2, ...). Sorry, I happened to post simultaneously and did not mean to "bury" your comment :sweat_smile:

[hersle]

Regarding @TorkelE's example with

oprob.ps[X0] = 3.0 # with default X = X0, should now oprob[X] == 3.0?

vs.

oprob = remake(oprob; p = [X0 => 3.0], u0 = Dict(), use_defaults = true) # with default X = X0, now oprob[X] == 3.0

As a user, I would not expect oprob[X] == 3.0 after the first method. What if the system is huge and has many defaults? For performance reasons, I think it is somewhat unfair to expect the library to reinitialize the system and recalculate all dependent defaults for every single line in

oprob.ps[X0] = 3.0
oprob.ps[Y0] = 4.0
oprob.ps[Z0] = 5.0
# perhaps many more in a bigger system ...

Instead, I think the second method with one call to remake communicates very clearly that it will reinitialize the problem once in a performant way. In other words, I would argue that you are "on your own" when modifying oprob.ps or oprob.u0 of a problem with dependent defaults.

[TorkelE]

If there are performance involved, I think @hersle makes a good point. To keep setu and setp performant that might be required.

Generally, in Catalyst, I think we recommend that remake is used for all modifications to systems, and then describe setp and setu to the user as something for situations where performance is important.

[AayushSabharwal]

@TorkelE

The problem So the problem is that, right now, the value of parameters that have initial conditions are set when a ODEProblem is declared, but never updated when setu or remake is called. That means that they have the wrong value when you create an integrator from the ODEProblem. Basically, what is needed, is that when setu or remake is called, the value of unknowns/parameters that depend on the updated thing must also be updated.

So there is no problem? If you have a default for Gamma = X1 + X2, it will be used as an equation during initialization. If the user specifies X1 and X2 to the ODEProblem during creation, Gamma will default to X1 + X2 and be updated to X1 + X2 during initialization, in case either of X1 or X2 change. If Gamma changes, I'm not entirely sure what will happen.

While the initial example I gave her might be a bit non-intuitive, I think it should be clear why the behaviour in this example is un-desired.

I can see why this is undesirable, but it's undesirable specifically because of the semantics Catalyst wants to enforce. Your concern is that the values in the ODEProblem don't reflect the constraints enforced, but the thing is they're not meant to and never have. These algebraic constraints are enforced during initialization. The ODEProblem just contains what the user knows about the initial state, and how the initialization should update said state to be consistent with the aforementioned algebraic constraints. Updating when a value changes based on defaults is

1. slow
   • Imagine a large system, with many such dependent defaults. Every single time you modify any part of the ODEProblem, all the dependent defaults would have to be recalculated. This is simply too slow.
2. Not technically correct
   • What if I change Gamma in the example above? What if there are more complicated constraints that can't be consistently enforced by simple substitution? What if changing a variable violates algebraic constraints?

As long as you have the default Gamma = X1 + X2, it will be enforced (to the best of the solver's ability, in case of underdetermined or overdetermined initialization) when the problem is initialized.

[AayushSabharwal]

@hersle

but is there a common need for parameter_dependencies?

Yes. Parameters can be changed during callbacks. parameter_dependencies updates the dependent parameters when this happens. It's also less prone to errors, since defaults don't actually enforce a dependency, they just assign a default value during initialization which can be overriden. With a default of p2 => 2p1 both p1 and p2 are true parameters, and can be changed independently. If I modify p2, it may not equal p1. Parameter dependencies specifically disallow this behavior. If p2 => 2p1 is a parameter dependency, if I run prob.ps[p1] = 2.0 then prob.ps[p2] will be 4.0.

This all works for me on master, except the line with remake(prob2, ...), as reported in https://github.com/SciML/ModelingToolkit.jl/issues/2715

That is a bug, and something I intend to fix eventually. remake should respect these defaults if asked to.

[hersle]

Thank you! I did not think of that use, and this clears up all my confusion. For reference, it sounds to me like the TL;DR of this thread is:

• To enforce relationships between variables and/or parameters during initialization, use defaults.
• To enforce relationships between parameters at all times (e.g. when one is modified in a callback), use parameter_dependencies.
• To update a problem with new initial conditions or parameters that respect defaults, use remake(prob; u0 = [x1 => 1.0, ...], p = [p1 => 1.0, ...], use_defaults = true).
• Defaults are not respected when modifying e.g. prob.ps[p1] = 1.0. Doing so would generally be too slow and not flexible enough to support more complicated constraints.
• Currently, remake does not properly handle defaults where variables depend on parameters. This is a bug. See examples in this thread or #2715 .

[AayushSabharwal]

Yeah. Also an additional point - the relationships during initialization won't be fully respected until https://github.com/SciML/ModelingToolkit.jl/issues/2665 is closed. Specifically, relationships that require parameters to be changed won't work.

[TorkelE]

I think under ideal circumstances we could have a dependency graph to ensure that defaults are updated, but I understand that that is lots of work given the limited resources that we have.

In the case a ODEProblem represents not the initial condition of the simulation, but simply the inputs, shouldn't

using ModelingToolkit
using OrdinaryDiffEq
@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)

yield nothing or similar when you type

oprob.ps[:d]

, since no value was given to ODEProblem? I.e. generally I think that this equality should hold:

using ModelingToolkit
@variables t X(t)
@parameters d = X
D = Differential(t)
eqs = [D(X) ~ d*X]
@mtkbuild osys = ODESystem(eqs, t)
oprob = ODEProblem(osys, [X => 1.0], (0.0, 1.0), [])

u0 = [X => 2.0]
tspan = (0.0, 1.0)
p = []
oprob1 = remake(oprob; u0, tspan, p)
oprob2 = ODEProblem(osys, u0, tspan, p)
oprob1 == oprob2

(i.e. using remake or ODEProblem with the same input should give the same results)

[TorkelE]

For now, having intended behaviour in actual simulations should be good, but I am correct that this does not work for defaults using either parameters/variables? I.e. this

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

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

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

oprob.ps[X0] # 2.0
oprob[X] # 2.0
oprob.ps[X0] = 3.0
oprob[X] # 2.0, but should be 3.0

using Plots
sol = solve(oprob)
plot(sol)

gives an initial value of X to 2.0, not 3.0. It does not seem to be exactly the issue linked, but maybe it is under the hood?

[AayushSabharwal]

yield nothing or similar when you type

Why? You gave it X, it calculated d based on the defaults because it wasn't given another value. Actually respecting X ~ d is something that isn't implemented yet, and even then nothing won't be returned, it'll be 0 at worst. Storing and/or returning nothing is also a lot of work.

I.e. generally I think that this equality should hold:

It will if you pass use_defaults = true to remake. It's unfortunate, but it's required since the alternative is breaking for SciMLBase.

For now, having intended behaviour in actual simulations should be good, but I am correct that this does not work for defaults using either parameters/variables?

It's the same thing. The change will (should) happen when the ODEProblem is initialized.

[TorkelE]

But you said that ODEProblem does not represent the initial value, but the value which was provided to ODEProblem? If that is the case, and that the default does not represent the default value, but the default value at initiation, then d has not value, and no value should be returned (rather than its default value).

If instead ODEProblem represent the current value that should be used, then it makes sense to update the value of d when the values its default depend on updates. Then I can accept that implementing a proper dependency graph and check for this is not worth the effort now, however, I am just trying to wrap my head around what ODEProblem actually represents, because it seems it might have changed.

[isaacsas]

If I am understanding this correctly won’t this initialization approach require all integrator type objects across all of SciML to add such code? I.e, this will require updates to OrdinaryDiffEq, StochasticDiffEq, JumpProcesses, etc. Is this planned? It would be unfortunate if only ODE models supported defining parameters and initial conditions in terms of specified initial conditions given this feature makes sense for many mathematical model types.

[AayushSabharwal]

But you said that ODEProblem does not represent the initial value, but the value which was provided to ODEProblem?

Yes, but it also needs to create the parameter object and store the value. Storing nothing is type unstable (at best) and definitely breaking

If I am understanding this correctly won’t this initialization approach require all integrator type objects across all of SciML to add such code?

Unfortunately, yes. ODEs support it, but I'm not aware of the required semantics for the rest of the ecosystem to be able to do so.

[TorkelE]

SDEProblem, DiscreteProblem, NonlinearProblem, and SteadyStateProblems should all work similarily? JumpProblems might be slightly modified since they store this inside a DiscretProblem. OptimizationProblems I am unsure about.

[isaacsas]

I guess the minimal semantics I would want as an MTK user for non-DAE models are pretty basic; I can have a mix of parameters and initial conditions of two types; some are "independent" and I fully specify them numerically via mappings or defaults, and some are "dependent" and defined via explicit formulas involving a mix of independent parameters and independent initial conditions. It seems like such functionality should be able to be provided in a universal way and not require customization at the level of each integrator (though perhaps integrators need to be modified to call some common SciMLBase/DiffEqBase initialization routines). (But I also admit @AayushSabharwal that you have certainly thought about this much more than I, so perhaps I am missing fundamental issues in providing such functionality integrators could easily opt into.)