# ModelingToolkit: Indexing a parameter by another parameter in an ODEProblem

**URL:** <https://discourse.julialang.org/t/modelingtoolkit-indexing-a-parameter-by-another-parameter-in-an-odeproblem/83569>\
**Category:** General Usage\
**Tags:** modelingtoolkit, differentialequation\
**Created:** [June 30, 2022, 2:54pm UTC](https://discourse.julialang.org/t/modelingtoolkit-indexing-a-parameter-by-another-parameter-in-an-odeproblem/83569 "2022-06-30T14:54:20Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![DanielVandH](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielvandh/32/31134_2.png) [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)\
**Post date:** [June 30, 2022, 2:54pm UTC](https://discourse.julialang.org/t/modelingtoolkit-indexing-a-parameter-by-another-parameter-in-an-odeproblem/83569/1 "2022-06-30T14:54:20Z")

</div>

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

```julia
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

```julia
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

```julia
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
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.

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [June 30, 2022, 8:21pm UTC](https://discourse.julialang.org/t/modelingtoolkit-indexing-a-parameter-by-another-parameter-in-an-odeproblem/83569/2 "2022-06-30T20:21:04Z")

</div>

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

```julia
lookup_kpa(dict, index) = dict[index]

```

Register that function so it is not traced:

```julia
@register_symbolic lookup_kpa(dict, index)

```

And then use that instead of directly indexing.

```julia
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.

---

<div class="post-metadata">

**Author:** ![DanielVandH](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielvandh/32/31134_2.png) [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)\
**Post date:** [July 1, 2022, 3:08am UTC](https://discourse.julialang.org/t/modelingtoolkit-indexing-a-parameter-by-another-parameter-in-an-odeproblem/83569/3 "2022-07-01T03:08:26Z")

</div>

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,

```julia
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
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`.

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [July 1, 2022, 3:54pm UTC](https://discourse.julialang.org/t/modelingtoolkit-indexing-a-parameter-by-another-parameter-in-an-odeproblem/83569/4 "2022-07-01T15:54:25Z")

</div>

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:

```julia
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

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

```