Nicer parameter indexing

Consider the following example:

import Plots
using ModelingToolkit, DifferentialEquations

@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t

eqs = [D(x) ~ σ * (y - x),
       D(y) ~ x * (ρ - z) - y,
       D(z) ~ x * y - β * z]

lorenz1 = ODESystem(eqs, name=:lorenz1)
lorenz2 = ODESystem(eqs, name=:lorenz2)

@variables a
@parameters γ
connections = [0 ~ lorenz1.x + lorenz2.y + a * γ]
connected = ODESystem(connections, t, [a], [γ], systems=[lorenz1,lorenz2])

u0 = [lorenz1.x => 1.0,
      lorenz1.y => 0.0,
      lorenz1.z => 0.0,
      lorenz2.x => 0.0,
      lorenz2.y => 1.0,
      lorenz2.z => 0.0,
      a => 2.0]

The parameters are than specified as a vector of pairs, corresponding to the parameters of the ODESystem.

p  = [lorenz1.σ => 10.0,
      lorenz1.ρ => 28.0,
      lorenz1.β => 8 / 3,
      lorenz2.σ => 10.0,
      lorenz2.ρ => 28.0,
      lorenz2.β => 8 / 3,
      γ => 2.0]

tspan = (0.0, 100.0)
prob = ODEProblem(connected, u0, tspan, p)
sol = solve(prob, Rodas5())

It feels a bit inconsistent to speficy p as a vector of pairs but inside the integrator it is just a vector.
In a callback one would have to use e.g. int.p[1] which in this case corresponds to gamma.

julia> prob.p
7-element Array{Float64,1}:
  2.0
 10.0
 28.0
  2.6666666666666665
 10.0
 28.0
  2.6666666666666665

It seems to be constructed like this:

[
    connected.ps...,
    connected.systems[1].ps...,
    connected.systems[2].ps...,
]

Is there an easier way of referencing the parameters e.g. for usage in a callback function?

1 Like

We could use ComponentArrays.jl to make this better, and that might be what we do in the future.

2 Likes

Has this been implemented yet (if so, example?) or if not is there an issue that I could subscribe to? A search turned up https://github.com/SciML/ModelingToolkit.jl/pull/702 which I think mentions this in passing but I can’t tell if it’s any more than that.

You can do it like:

https://github.com/SciML/ModelingToolkit.jl/issues/318#issuecomment-796632465

and that should get automated hopefully soon.

1 Like