Using the SciML Symbolic Indexing Interface

This tutorial will cover ways to use the symbolic indexing interface for types that implement it. SciML's core types (problems, solutions, and iterator (integrator) types) all support this symbolic indexing interface which allows for domain-specific interfaces (such as ModelingToolkit, Catalyst, etc.) to seamlessly blend their symbolic languages with the types obtained from SciML. Other tutorials will focus on how users can make use of the interface for their own DSL, this tutorial will simply focus on what the user experience looks like for DSL which have set it up.

We recommend any DSL implementing the symbolic indexing interface to link to this tutorial as a full description of the functionality.

Note

While this tutorial focuses on demonstrating the symbolic indexing interface for ODEs, note that the same functionality works across all of the other problem types, such as optimization problems, nonlinear problems, nonlinear solutions, etc.

Symbolic Indexing of Differential Equation Solutions

Consider the following example:

using ModelingToolkit, OrdinaryDiffEq, SymbolicIndexingInterface, Plots
using ModelingToolkit: t_nounits as t, D_nounits as D

@parameters σ ρ β
@variables x(t) y(t) z(t) w(t)

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

@mtkbuild sys = ODESystem(eqs, t)

\[ \begin{align} \frac{\mathrm{d} \mathtt{xˍt}\left( t \right)}{\mathrm{d}t} &= \left( - x\left( t \right) + y\left( t \right) \right) \sigma \\ \frac{\mathrm{d} y\left( t \right)}{\mathrm{d}t} &= - y\left( t \right) + x\left( t \right) \left( - z\left( t \right) + \rho \right) \\ \frac{\mathrm{d} z\left( t \right)}{\mathrm{d}t} &= x\left( t \right) y\left( t \right) - z\left( t \right) \beta \\ \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} &= \mathtt{xˍt}\left( t \right) \end{align} \]

The system has 4 state variables, 3 parameters, and one observed variable:

ModelingToolkit.observed(sys)

\[ \begin{align} w\left( t \right) &= x\left( t \right) + y\left( t \right) + z\left( t \right) \end{align} \]

Solving the system,

u0 = [D(x) => 2.0,
    x => 1.0,
    y => 0.0,
    z => 0.0]

p = [σ => 28.0,
    ρ => 10.0,
    β => 8 / 3]

tspan = (0.0, 100.0)
prob = ODEProblem(sys, u0, tspan, p, jac = true)
sol = solve(prob, Tsit5())

retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 1495-element Vector{Float64}:
   0.0
   0.00014131524176735154
   0.0015544676594408668
   0.012956381586203063
   0.038266699884064254
   0.0785229194408219
   0.12210864510623465
   0.1715537595773225
   0.23112906724447269
   0.30252384441787916
   ⋮
  99.50500400949036
  99.57665671063391
  99.64283644778946
  99.72166970896653
  99.79967565070119
  99.86801146250268
  99.93304886418163
  99.99013920811814
 100.0
u: 1495-element Vector{Vector{Float64}}:
 [2.0, 0.0, 0.0, 1.0]
 [1.9960454103704428, 0.0014132521265988688, 9.986094350108851e-8, 1.0002823510089436]
 [1.956746201469907, 0.01555657123298723, 1.2096178746316906e-5, 1.0030752466310615]
 [1.6563911009012329, 0.13029732325873117, 0.0008464888055104521, 1.0236450170514986]
 [1.1011942450107057, 0.3874130541723886, 0.007445311072875279, 1.0582098434939375]
 [0.5557298342165394, 0.795954115942288, 0.031202922260444813, 1.0901348112728222]
 [0.44950103893773047, 1.2290079071325886, 0.07390401064741846, 1.1102092737279297]
 [0.9266727130791688, 1.706238215969287, 0.14159578142525983, 1.1416871413493421]
 [2.278100719563122, 2.2773793680509047, 0.2511409516498877, 1.2331771178073494]
 [4.867420100033599, 3.0146936570333036, 0.44047996356952857, 1.482457085693269]
 ⋮
 [80.19626954598822, 0.2231933697492472, 6.377997673391619, -10.921009114194087]
 [93.91126317449019, -1.943511875150559, 5.714426336544729, -4.580903836595237]
 [92.47052726207596, -2.2128840651669464, 4.966588950256005, 1.6531217356898207]
 [78.34377161396489, 0.27369559171564467, 3.6720084014137733, 8.448762638974289]
 [59.61988063036174, 5.2078889703092095, 5.364597977976095, 13.8343282979256]
 [42.06608076700321, 6.848720789431849, 10.923031490807823, 17.327761617672973]
 [17.746354561947708, 2.5614609605308467, 14.71014279775863, 19.336066112359077]
 [-13.425210358521637, -1.9987158362085793, 12.635789047318434, 19.496276813392893]
 [-19.399440458388117, -2.411052621022261, 11.888639580483389, 19.334495167624354]

We can obtain the timeseries of any time-dependent variable using getindex

sol[x]

1495-element Vector{Float64}:
   1.0
   1.0002823510089436
   1.0030752466310615
   1.0236450170514986
   1.0582098434939375
   1.0901348112728222
   1.1102092737279297
   1.1416871413493421
   1.2331771178073494
   1.482457085693269
   ⋮
 -10.921009114194087
  -4.580903836595237
   1.6531217356898207
   8.448762638974289
  13.8343282979256
  17.327761617672973
  19.336066112359077
  19.496276813392893
  19.334495167624354

This also works for arrays or tuples of variables, including observed quantities and independent variables, for interpolating solutions, and plotting:

sol[[x, y]]

1495-element Vector{Vector{Float64}}:
 [1.0, 0.0]
 [1.0002823510089436, 0.0014132521265988688]
 [1.0030752466310615, 0.01555657123298723]
 [1.0236450170514986, 0.13029732325873117]
 [1.0582098434939375, 0.3874130541723886]
 [1.0901348112728222, 0.795954115942288]
 [1.1102092737279297, 1.2290079071325886]
 [1.1416871413493421, 1.706238215969287]
 [1.2331771178073494, 2.2773793680509047]
 [1.482457085693269, 3.0146936570333036]
 ⋮
 [-10.921009114194087, 0.2231933697492472]
 [-4.580903836595237, -1.943511875150559]
 [1.6531217356898207, -2.2128840651669464]
 [8.448762638974289, 0.27369559171564467]
 [13.8343282979256, 5.2078889703092095]
 [17.327761617672973, 6.848720789431849]
 [19.336066112359077, 2.5614609605308467]
 [19.496276813392893, -1.9987158362085793]
 [19.334495167624354, -2.411052621022261]

sol[(t, w)]

1495-element Vector{Tuple{Float64, Float64}}:
 (0.0, 1.0)
 (0.00014131524176735154, 1.001695702996486)
 (0.0015544676594408668, 1.0186439140427952)
 (0.012956381586203063, 1.1547888291157402)
 (0.038266699884064254, 1.4530682087392015)
 (0.0785229194408219, 1.917291849475555)
 (0.12210864510623465, 2.413121191507937)
 (0.1715537595773225, 2.989521138743889)
 (0.23112906724447269, 3.7616974375081416)
 (0.30252384441787916, 4.937630706296101)
 ⋮
 (99.50500400949036, -4.319818071053221)
 (99.57665671063391, -0.8099893752010663)
 (99.64283644778946, 4.4068266207788795)
 (99.72166970896653, 12.394466632103708)
 (99.79967565070119, 24.406815246210904)
 (99.86801146250268, 35.099513897912644)
 (99.93304886418163, 36.60766987064855)
 (99.99013920811814, 30.13335002450275)
 (100.0, 28.812082127085482)

sol(1.3, idxs=x)

-7.739577070673492

sol(1.3, idxs=[x, w])

2-element Vector{Float64}:
 -7.739577070673492
 -5.727345621422577

sol(1.3, idxs=[:y, :z])

2-element Vector{Float64}:
 -1.688874681282943
  3.7011061305338573

plot(sol, idxs=x)

If necessary, Symbols can be used to refer to variables. This is only valid for symbolic variables for which hasname returns true. The Symbol used must match the one returned by getname for the variable.

hasname(x)

true

getname(x)

:x

sol[(:x, :w)]

1495-element Vector{Tuple{Float64, Float64}}:
 (1.0, 1.0)
 (1.0002823510089436, 1.001695702996486)
 (1.0030752466310615, 1.0186439140427952)
 (1.0236450170514986, 1.1547888291157402)
 (1.0582098434939375, 1.4530682087392015)
 (1.0901348112728222, 1.917291849475555)
 (1.1102092737279297, 2.413121191507937)
 (1.1416871413493421, 2.989521138743889)
 (1.2331771178073494, 3.7616974375081416)
 (1.482457085693269, 4.937630706296101)
 ⋮
 (-10.921009114194087, -4.319818071053221)
 (-4.580903836595237, -0.8099893752010663)
 (1.6531217356898207, 4.4068266207788795)
 (8.448762638974289, 12.394466632103708)
 (13.8343282979256, 24.406815246210904)
 (17.327761617672973, 35.099513897912644)
 (19.336066112359077, 36.60766987064855)
 (19.496276813392893, 30.13335002450275)
 (19.334495167624354, 28.812082127085482)

Note how when indexing with an array, the returned type is a Vector{Array{Float64}}, and when using a Tuple, the returned type is Vector{Tuple{Float64, Float64}}. To obtain the value of all state variables, we can use the shorthand:

sol[solvedvariables] # equivalent to sol[variable_symbols(sol)]

1495-element Vector{Vector{Float64}}:
 [2.0, 0.0, 0.0, 1.0]
 [1.9960454103704428, 0.0014132521265988688, 9.986094350108851e-8, 1.0002823510089436]
 [1.956746201469907, 0.01555657123298723, 1.2096178746316906e-5, 1.0030752466310615]
 [1.6563911009012329, 0.13029732325873117, 0.0008464888055104521, 1.0236450170514986]
 [1.1011942450107057, 0.3874130541723886, 0.007445311072875279, 1.0582098434939375]
 [0.5557298342165394, 0.795954115942288, 0.031202922260444813, 1.0901348112728222]
 [0.44950103893773047, 1.2290079071325886, 0.07390401064741846, 1.1102092737279297]
 [0.9266727130791688, 1.706238215969287, 0.14159578142525983, 1.1416871413493421]
 [2.278100719563122, 2.2773793680509047, 0.2511409516498877, 1.2331771178073494]
 [4.867420100033599, 3.0146936570333036, 0.44047996356952857, 1.482457085693269]
 ⋮
 [80.19626954598822, 0.2231933697492472, 6.377997673391619, -10.921009114194087]
 [93.91126317449019, -1.943511875150559, 5.714426336544729, -4.580903836595237]
 [92.47052726207596, -2.2128840651669464, 4.966588950256005, 1.6531217356898207]
 [78.34377161396489, 0.27369559171564467, 3.6720084014137733, 8.448762638974289]
 [59.61988063036174, 5.2078889703092095, 5.364597977976095, 13.8343282979256]
 [42.06608076700321, 6.848720789431849, 10.923031490807823, 17.327761617672973]
 [17.746354561947708, 2.5614609605308467, 14.71014279775863, 19.336066112359077]
 [-13.425210358521637, -1.9987158362085793, 12.635789047318434, 19.496276813392893]
 [-19.399440458388117, -2.411052621022261, 11.888639580483389, 19.334495167624354]

This does not include the observed variable w. To include observed variables in the output, the following shorthand is used:

sol[allvariables] # equivalent to sol[all_variable_symbols(sol)]

1495-element Vector{Vector{Float64}}:
 [2.0, 0.0, 0.0, 1.0, 1.0]
 [1.9960454103704428, 0.0014132521265988688, 9.986094350108851e-8, 1.0002823510089436, 1.001695702996486]
 [1.956746201469907, 0.01555657123298723, 1.2096178746316906e-5, 1.0030752466310615, 1.0186439140427952]
 [1.6563911009012329, 0.13029732325873117, 0.0008464888055104521, 1.0236450170514986, 1.1547888291157402]
 [1.1011942450107057, 0.3874130541723886, 0.007445311072875279, 1.0582098434939375, 1.4530682087392015]
 [0.5557298342165394, 0.795954115942288, 0.031202922260444813, 1.0901348112728222, 1.917291849475555]
 [0.44950103893773047, 1.2290079071325886, 0.07390401064741846, 1.1102092737279297, 2.413121191507937]
 [0.9266727130791688, 1.706238215969287, 0.14159578142525983, 1.1416871413493421, 2.989521138743889]
 [2.278100719563122, 2.2773793680509047, 0.2511409516498877, 1.2331771178073494, 3.7616974375081416]
 [4.867420100033599, 3.0146936570333036, 0.44047996356952857, 1.482457085693269, 4.937630706296101]
 ⋮
 [80.19626954598822, 0.2231933697492472, 6.377997673391619, -10.921009114194087, -4.319818071053221]
 [93.91126317449019, -1.943511875150559, 5.714426336544729, -4.580903836595237, -0.8099893752010663]
 [92.47052726207596, -2.2128840651669464, 4.966588950256005, 1.6531217356898207, 4.4068266207788795]
 [78.34377161396489, 0.27369559171564467, 3.6720084014137733, 8.448762638974289, 12.394466632103708]
 [59.61988063036174, 5.2078889703092095, 5.364597977976095, 13.8343282979256, 24.406815246210904]
 [42.06608076700321, 6.848720789431849, 10.923031490807823, 17.327761617672973, 35.099513897912644]
 [17.746354561947708, 2.5614609605308467, 14.71014279775863, 19.336066112359077, 36.60766987064855]
 [-13.425210358521637, -1.9987158362085793, 12.635789047318434, 19.496276813392893, 30.13335002450275]
 [-19.399440458388117, -2.411052621022261, 11.888639580483389, 19.334495167624354, 28.812082127085482]

Evaluating expressions

getsym also generates functions for expressions if the object passed to it supports observed. For example:

getsym(prob, x + y + z)(prob)

1.0

To evaluate this function using values other than the ones contained in prob, we need an object that supports state_values, parameter_values, current_time. SymbolicIndexingInterface provides the ProblemState type, which has trivial implementations of the above functions. We can thus do:

temp_state = ProblemState(; u = [0.1, 0.2, 0.3, 0.4], p = parameter_values(prob))
getsym(prob, x + y + z)(temp_state)

0.8999999999999999

Note that providing all of the state vector, parameter object and time may not be necessary if the function generated by observed does not access them. ModelingToolkit.jl generates functions that access the parameters regardless of whether they are used in the expression, and thus it needs to be provided to the ProblemState.

Parameter Indexing: Getting and Setting Parameter Values

Parameters cannot be obtained using this syntax, and instead require using getp and setp.

Note

The reason why parameters use a separate syntax is to be able to ensure type stability of the sol[x] indexing. Without separating the parameter indexing, the return type of symbolic indexing could be anything a parameter can be, which is general is not the same type as state variables!

σ_getter = getp(sys, σ)
σ_getter(sol) # can also pass `prob`

28.0

Note that this also supports arrays/tuples of parameter symbols:

σ_ρ_getter = getp(sys, (σ, ρ))
σ_ρ_getter(sol)

(28.0, 10.0)

Now suppose the system has to be solved with a different value of the parameter β.

β_setter = setp(sys, β)
β_setter(prob, 3)

The updated parameter values can be checked using parameter_values.

parameter_values(prob)

ModelingToolkit.MTKParameters{Vector{Float64}, Vector{Float64}, Tuple{}, Tuple{}, Tuple{}, Tuple{}}([10.0, 3.0, 28.0], [0.0, 0.0, 1.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0], (), (), (), ())

When solving the new system, note that the parameter getter functions still work on the new solution object.

sol2 = solve(prob, Tsit5())
σ_getter(sol)

28.0

σ_ρ_getter(sol)

(28.0, 10.0)

To set the entire parameter vector at once, setp can be used (note that the order of symbols passed to setp must match the order of values in the array).

setp(prob, parameter_symbols(prob))(prob, [29.0, 11.0, 2.5])
parameter_values(prob)

ModelingToolkit.MTKParameters{Vector{Float64}, Vector{Float64}, Tuple{}, Tuple{}, Tuple{}, Tuple{}}([29.0, 11.0, 2.5], [0.0, 0.0, 1.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0], (), (), (), ())

Note

These getters and setters generate high-performance functions for the specific chosen symbols or collection of symbols. Caching the getter/setter function and reusing it on other problem/solution instances can be the key to achieving good performance. Note that this caching is allowed only when the symbolic system is unchanged (it's fine for the numerical values to have changed, but not the underlying equations).

Re-creating a buffer

To re-create a buffer (of unknowns or parameters) use remake_buffer. This allows changing the type of values in the buffer (for example, changing the value of a parameter from Float64 to Float32).

remake_buffer(sys, prob.p, Dict(σ => 1f0, ρ => 2f0, β => 3f0))

ModelingToolkit.MTKParameters{Vector{Float32}, Vector{Float64}, Tuple{}, Tuple{}, Tuple{}, Tuple{}}(Float32[2.0, 3.0, 1.0], [0.0, 0.0, 1.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0], (), (), (), ())