ModelingToolkit: Indexing a parameter by another parameter in an ODEProblem

General Usage
,
1

I have an ODEProblem where I need to index into a parameter using another separate parameter. Consider for example

using ModelingToolkit, DifferentialEquations
function rober(du, u, p, t)
    y₁, y₂, y₃ = u
    k₁, k₂, k₃, d, v = p
    kpa = d[v]
    du[1] = -k₁ * y₁ + k₃ * y₂ * y₃
    du[2] = k₁ * y₁ - k₂ * y₂^2 - k₃ * y₂ * y₃
    du[3] = k₂ * y₂^2
    nothing
end
just_a_dict = Dict([1,2,3,4,5,10,29] .=> [[1,2,3], [1,2,3], [4,5,1], [2,3,1], [4,1,9],[5,1,2],[2,3,1]])
prob = ODEProblem(rober, [1.0, 0.0, 0.0], (0.0, 1e5), (0.04, 3e7, 1e4, just_a_dict, 1))
modelingtoolkitize(prob)

This leads to the following

ERROR: MethodError: no method matching getindex(::Num, ::Num)
Closest candidates are:
  getindex(::Number) at C:\Users\licer\AppData\Local\Programs\Julia-1.7.3\share\julia\base\number.jl:95
  getindex(::Number, ::Integer) at C:\Users\licer\AppData\Local\Programs\Julia-1.7.3\share\julia\base\number.jl:96
  getindex(::Number, ::Integer...) at C:\Users\licer\AppData\Local\Programs\Julia-1.7.3\share\julia\base\number.jl:101
  ...
Stacktrace:
 [1] rober(du::Vector{Num}, u::Vector{Num}, p::NTuple{5, Num}, t::Num)
   @ Main c:\Users\licer\.julia\dev\FVM\dev\sparsity.jl:25
 [2] (::ODEFunction{true, typeof(rober), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing})(::Vector{Num}, ::Vararg{Any})
   @ SciMLBase C:\Users\licer\.julia\packages\SciMLBase\IJbT7\src\scimlfunctions.jl:1624
 [3] modelingtoolkitize(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Tuple{Float64, Float64, Float64, Dict{Int64, Vector{Int64}}, Int64}, ODEFunction{true, typeof(rober), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})       
   @ ModelingToolkit C:\Users\licer\.julia\packages\ModelingToolkit\tMgaW\src\systems\diffeqs\modelingtoolkitize.jl:48
 [4] modelingtoolkitize(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Tuple{Float64, Float64, Float64, Dict{Int64, Vector{Int64}}, Int64}, ODEFunction{true, typeof(rober), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem})
   @ ModelingToolkit C:\Users\licer\.julia\packages\ModelingToolkit\tMgaW\src\systems\diffeqs\modelingtoolkitize.jl:7
 [5] top-level scope
   @ c:\Users\licer\.julia\dev\FVM\dev\sparsity.jl:33

Is there any way to allow this behaviour with kpa = d[v]?

I also have indexing like u[q[1]] or s[k][1], with q a provided variable and k::Num, which gives similar errors, like

using ModelingToolkit, DifferentialEquations
function rober(du, u, p, t)
    y₁, y₂, y₃ = u
    k₁, k₂, k₃, d, v, q = p
    u[q[1]]
    kpa = d[v]
    du[1] = -k₁ * y₁ + k₃ * y₂ * y₃
    du[2] = k₁ * y₁ - k₂ * y₂^2 - k₃ * y₂ * y₃
    du[3] = k₂ * y₂^2
    nothing
end
just_a_dict = Dict([1,2,3,4,5,10,29] .=> [[1,2,3], [1,2,3], [4,5,1], [2,3,1], [4,1,9],[5,1,2],[2,3,1]])
prob = ODEProblem(rober, [1.0, 0.0, 0.0], (0.0, 1e5), (0.04, 3e7, 1e4, just_a_dict, 1, [1, 2]))
modelingtoolkitize(prob)

julia> modelingtoolkitize(prob)
ERROR: ArgumentError: invalid index: α₅ of type Num

though this error is much less clear to me for how it might be resolved / whether it’s the same issue as above.

2

I think the solution is to just tell MTK not to trace through the lookup. Create a function that does the dictionary indexing:

lookup_kpa(dict, index) = dict[index]

Register that function so it is not traced:

@register_symbolic lookup_kpa(dict, index)

And then use that instead of directly indexing.

function ...
    ...
    kpa = lookup_kpa(d, v)
    ...
end

I don’t know of a way to tell modelingtoolkitize that a parameter has a type other than Num, so it is always going to look for the wrong getindex.

3

Thank you. This does work, but for the latter problem in my question where I’d need to index as e.g. u[v[1]]], this doesn’t seem to work anymore? For example,

using ModelingToolkit, DifferentialEquations, Symbolics
@inline lookup_val(dict, idx) = dict[idx]
@register_symbolic lookup_val(dict, idx)
function rober(du, u, p, t)
    y₁, y₂, y₃ = u
    k₁, k₂, k₃, d, v = p
    kpa = lookup_val(d, v)
    #kpa2 = u[kpa] # <-- here
    kpa2 = lookup_val(u, kpa)
    kpa2 = lookup_val(u, kpa[1])
    du[1] = -k₁ * y₁ + k₃ * y₂ * y₃
    du[2] = k₁ * y₁ - k₂ * y₂^2 - k₃ * y₂ * y₃
    du[3] = k₂ * y₂^2
    nothing
end
just_a_dict = Dict([1,2,3,4,5,10,29] .=> [[1,2,3], [1,2,3], [4,5,1], [2,3,1], [4,1,9],[5,1,2],[2,3,1]])
prob = ODEProblem(rober, [1.0, 0.0, 0.0], (0.0, 1e5), (0.04, 3e7, 1e4, just_a_dict, 1))
modelingtoolkitize(prob)

Neither of the two kpa2 lines seem to work, giving the error

julia> modelingtoolkitize(prob)
ERROR: ArgumentError: invalid index: lookup_val(α₄, α₅) of type Num
Stacktrace:
  [1] to_index(i::Num)
    @ Base .\indices.jl:300
  [2] to_index(A::Vector{Num}, i::Num)
    @ Base .\indices.jl:277
  [3] to_indices
    @ .\indices.jl:333 [inlined]
  [4] to_indices
    @ .\indices.jl:325 [inlined]
  [5] getindex(A::Vector{Num}, I::Num)
    @ Base .\abstractarray.jl:1218
  [6] lookup_val
    @ c:\Users\licer\.julia\dev\FVM\dev\sparsity.jl:54 [inlined]
  [7] rober(du::Vector{Num}, u::Vector{Num}, p::NTuple{5, Num}, t::Num)
    @ Main c:\Users\licer\.julia\dev\FVM\dev\sparsity.jl:61
  [8] (::ODEFunction{true, typeof(rober), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, 
Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing})(::Vector{Num}, ::Vararg{Any})
    @ SciMLBase C:\Users\licer\.julia\packages\SciMLBase\IJbT7\src\scimlfunctions.jl:1624
  [9] modelingtoolkitize(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Tuple{Float64, Float64, Float64, Dict{Int64, Vector{Int64}}, Int64}, ODEFunction{true, 
typeof(rober), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit C:\Users\licer\.julia\packages\ModelingToolkit\tMgaW\src\systems\diffeqs\modelingtoolkitize.jl:48
 [10] modelingtoolkitize(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Tuple{Float64, Float64, Float64, Dict{Int64, Vector{Int64}}, Int64}, ODEFunction{true, 
typeof(rober), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem})
    @ ModelingToolkit C:\Users\licer\.julia\packages\ModelingToolkit\tMgaW\src\systems\diffeqs\modelingtoolkitize.jl:7
 [11] top-level scope
    @ c:\Users\licer\.julia\dev\FVM\dev\sparsity.jl:70

I’m not sure what the difference here is. Do I need to somehow register lookup_val as being classed as an index? I imagine that’s not possible since everything is Num.

4

In this case, you are indexing with kpa[1] while kpa is assumed to be a scalar. modelingtoolkitize works by calling your function with a type Num that redefines operations to record an unevaluated symbolic representation. Since no return type was specified in the registration, it was assumed to be a scalar Num. You can specify a return type to sort this out, but you also have to specify that you want the symbolic scalar lookup to accept a vector:

lookup_arr(d, v) = d[v] 
@register_symbolic lookup_arr(d, v)::Vector{Num}
lookup_val(a, v) = a[v]
@register_symbolic lookup_val(a::Vector{Num}, v)

function ...
    ...
    kpa = lookup_arr(d, v)
    kpa2 = lookup_val(u, kpa[1])
    ...
end

You could also just hide all the indexing behind one function with something like

lookup_val(u, d, v) = u[d[v]]
@register_symbolic lookup_val(u::Vector{Num}, d, v)