varmap_to_vars output and MTKParameters in ODEProblem are Different

[timknab]

Using varmap_to_vars to get the lowered array does not seem to return the same ordering as MTKParameters within an ODEProblem.

using ModelingToolkit, Plots, DifferentialEquations, Test
using ModelingToolkit: t_nounits as t, D_nounits as D
@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]

@named sys = ODESystem(eqs, t)
simpsys = structural_simplify(sys)
tspan = (0.0, 10.0)
prob = ODEProblem(simpsys, [], tspan)

println(ModelingToolkit.varmap_to_vars([α => 1.5, β=>1.0, γ => 3.0, δ => 1.0, x0 => 2.0],parameters(simpsys)))
println(prob.p[1])

@assert ModelingToolkit.varmap_to_vars([α => 1.5, β=>1.0, γ => 3.0, δ => 1.0, x0 => 2.0],parameters(simpsys)) == prob.p[1]

[2.0, 1.5, 1.0, 1.0, 3.0]
[1.5, 1.0, 1.0, 3.0, 2.0]
ERROR: AssertionError: ModelingToolkit.varmap_to_vars([α => 1.5, β => 1.0, γ => 3.0, δ => 1.0, x0 => 2.0], parameters(simpsys)) == prob.p[1]

Here it just seems to have moved x0 to the head of the parameter vector, but I have seen it much more "jumbled" in a more complicated model.

I found this when I tried to use ModelingToolkit.get_u0_p to generate new u0 and p vectors to use in remake because remake does not currently seem to update u0 parametrically based on p.

u0new, pnew = ModelingToolkit.get_u0_p(simpsys, [], [α => 1.5, β=>1.0, γ => 3.0, δ => 1.0, x0 => 7.0])
println(pnew)

[7.0, 1.5, 1.0, 1.0, 3.0]

Again, the value for x0 is not in the same position as in the MTKParameters vector.

Environment info

Pkg.status()
  [0c46a032] DifferentialEquations v7.13.0
  [06fc5a27] DynamicQuantities v0.12.3
  [961ee093] ModelingToolkit v9.4.0
  [91a5bcdd] Plots v1.40.1

[timknab]

A fix seems to be to just use MTKParameters directly.

ModelingToolkit.MTKParameters(simpsys,[α => 1.5, β=>1.0, γ => 3.0, δ => 1.0, x0 => 2.0])

and for get_u0_p

u0new, pnew, defs = ModelingToolkit.get_u0_p(simpsys, [], [α => 1.5, β=>1.0, γ => 3.0, δ => 1.0, x0 => 7.0])
ModelingToolkit.MTKParameters(simpsys,defs)

Maybe the documentation should be updated to reflect this? Would also be great to have remake handle parametric u0 (#2103) which was marked closed by #2403 but doesn't seem to have solved the issue.

[AayushSabharwal]

println(prob.p[1])

MTKParameters is not meant to be indexed. The fact that this works is not public API, and it only does so because it makes codegen easier. Parameters do not have an ordering, since they aren't stored in a single array.

I found this when I tried to use ModelingToolkit.get_u0_p to generate new u0 and p vectors to use in remake because remake does not currently seem to update u0 parametrically based on p.

If you need to generate new u0 and p, I recommend get_u0 for u0 and MTKParameters for p (two steps instead of one). Although it would be better to avoid this entirely as much as possible. remake has been fixed as of SciML/SciMLBase.jl#644. It also supports initial conditions dependent on other parameters/initial conditions. If you find cases that do not work, please open an issue.

[ChrisRackauckas]

We should update the FAQ