Remake does not update parametric initial conditions

[timknab]

Remake still does not appear to update initial conditions that are dependent on parameters. Following up from this issue where supposedly this was fixed in #644

using ModelingToolkit, Plots, DifferentialEquations, Test
using ModelingToolkit: t_nounits as t, D_nounits as D
using SciMLBase
@parameters α=1.5 β=1.0 γ=3.0 δ=1.0 x0 = 2.0
@variables x(t) = x0 y(t) = 0.0

eqs = [D(x) ~ α * x - β * x * y
       D(y) ~ -γ * y + δ * x * y]

@mtkbuild sys = ODESystem(eqs, t)
tspan = (0.0, 10.0)
prob = ODEProblem(sys, [], tspan)
prob2 = SciMLBase.remake(prob, p = [x0 => 9.0])
println(prob2.u0)

[2.0, 0.0]

Would expect prob2.u0 = [9.0, 0.0]

Pkg.status()
  [0c46a032] DifferentialEquations v7.13.0
  [961ee093] ModelingToolkit v9.4.0
  [91a5bcdd] Plots v1.40.2
  [0bca4576] SciMLBase v2.30.0 `https://github.com/SciML/SciMLBase.jl.git#master`

versioninfo()
Julia Version 1.10.2
Commit bd47eca2c8a (2024-03-01 10:14 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (x86_64-apple-darwin22.4.0)
  CPU: 12 × Intel(R) Core(TM) i7-9750H CPU @ 2.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, skylake)
Threads: 12 default, 0 interactive, 6 GC (on 12 virtual cores)
Environment:
  JULIA_PKG_USE_CLI_GIT = true
  JULIA_EDITOR = code

[AayushSabharwal]

Thanks for pointing this out. I should have phrased the reply in the linked issue a bit better. remake right now supports parametric initial conditions if they are passed to remake. i.e. passing u0 = [x => 2p1], p => [p1 => 3] will set the value of x to 6. To support doing this for default values in the system, we need SciML/SymbolicIndexingInterface.jl#47 which is something I'll get around to finishing soon

[timknab]

Ah ok, sorry for the misunderstanding

I'll keep an eye on SciML/SymbolicIndexingInterface.jl#47 and can easily work around it using get_u0 in the meantime

Thanks!

[AayushSabharwal]

This now works:

julia> prob2 = remake(prob, u0 = Dict(), p = [x0 => 9.0], use_defaults = true)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 10.0)
u0: 2-element Vector{Float64}:
 9.0
 0.0

u0 = Dict() is necessary since otherwise it doesn't process any updates to u0. use_defaults = true is necessary so that the model defaults take precedence over the existing values of variables in the problem.

[ChrisRackauckas]

use_defaults = true is necessary so that the model defaults take precedence over the existing values of variables in the problem.

That shouldn't be the case, dependencies should always be handled?

[AayushSabharwal]

That would be breaking. Prior to this remake rework, it didn't take defaults into account. If we do that by default, it breaks MTK tests

[ChrisRackauckas]

I wouldn't call that breaking. It's a bug fix. I think it's clear that any user would expect u0[1] = x0 here, and the fact that it doesn't after remake is unintended.

[AayushSabharwal]

So if the defaults also include constant values y => 10.0 should this be used instead of the existing values in the problem passed to remake? Or do we only consider dependent defaults?

[AayushSabharwal]

Also, I'm calling it breaking because MTK tests rely on this behavior of ignoring defaults. I think they'll also pass if we only ignore constant-valued defaults, but am not sure

[ChrisRackauckas]

I think we'd ideally error if any given constraints cannot be satisfied. I'm not sure how easy that is to do though, but that would be the most ideal situation.

[AayushSabharwal]

Should y => 1.0 be treated as a constraint? Wouldn't this severely restrict how problems can be remade?

[ChrisRackauckas]

No, just functions of other parameters.

[AayushSabharwal]

And if the user explicitly provides an alternate value (e.g. x => 1.0 in the above example)? Error?

[ChrisRackauckas]

I'd prefer error yes

[ChrisRackauckas]

Actually, it should act like the initialization stuff. So you always want to make sure the "current" equations are satisfied. Anything put into ODEProblem or remake is an override. If someone does p1 ~ 2*p2, then that should always be satisfied, so p2 => 2.0 should change p1. If someone then does p1 => 2.0, they are overriding this definition, since they used to have @parameters p1 = 2*p2 and now they did a late override to remove that equation.

[AayushSabharwal]

"Removing that equation" is modifying the system, which remake doesn't do. Parameter dependencies are intrinsic to the system since they only exist at the MTK level. Parameter defaults can be overridden, which remake will do.