Support more general parameters

[ranocha]

For example, the following code fails using Julia v1.6.0-rc3 and ModelingToolkit v5.13.6.

julia> using ModelingToolkit, OrdinaryDiffEq

julia> struct MyParameter
           a::Float64
       end

julia> f(du, u, p, t) = du .= p.a .* u
f (generic function with 1 method)

julia> ode = ODEProblem(f, ones(2), (0.0, 1.0), MyParameter(2.0))
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 1.0)
u0: [1.0, 1.0]

julia> sol = solve(ode, Tsit5())# the problem is set up correctly

julia> ode_mtk = modelingtoolkitize(ode)
ERROR: MethodError: no method matching keys(::MyParameter)
Closest candidates are:
  keys(::Union{Tables.AbstractColumns, Tables.AbstractRow}) at ~/.julia/packages/Tables/UxLRG/src/Tables.jl:179
  keys(::IndexStyle, ::AbstractArray, ::AbstractArray...) at abstractarray.jl:327
  keys(::IterativeSolvers.ConvergenceHistory) at ~/.julia/packages/IterativeSolvers/TpeDx/src/history.jl:223
  ...
Stacktrace:
 [1] eachindex(itrs::MyParameter)
   @ Base ./abstractarray.jl:253
 [2] modelingtoolkitize(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, MyParameter, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/7KM6m/src/systems/diffeqs/modelingtoolkitize.jl:21
 [3] top-level scope
   @ REPL[9]:1

Is it possible to use ModelingToolkit.jl to deal with general structs as parameter (p in the ODEProblem)?

[ChrisRackauckas]

#538 and #197 are related.

[ChrisRackauckas]

Note that it's extremely easy to workaround this by hand though.

using ModelingToolkit, OrdinaryDiffEq
struct MyParameter{T}
    a::T
end
Base.eachindex(p::MyParameter) = Base.OneTo(1)
Base.size(p::MyParameter) = (1,)
Base.length(p::MyParameter) = 1
Base.iterate(p::MyParameter,x=1) = Base.iterate([p.a],x)

f(du, u, p, t) = du .= p.a .* u
f(du, u, p::Array, t) = f(du, u, MyParameter(p...), t)
ode = ODEProblem(f, ones(2), (0.0, 1.0), MyParameter(2.0))
modelingtoolkitize(ode)

It's much harder to programmatically figure out where all of the parameter values are if than it is for the user to just give it. So it's going to be a difficult one to fix, while being rather easy for anyone to get a practical solution, so this one may sit for awhile.

[ranocha]

Thanks for the help, @ChrisRackauckas! So do I understand it correctly that modelingtoolkitize somehow needs to visit every scalar parameter? How can other types of parameters be handled that do not reduce to (sets of) scalars, e.g. functions? For example

Without changes at the user side (same error as above)

julia> using ModelingToolkit, OrdinaryDiffEq

julia> function f(du, u, p, t)
           @. du = -p(u)
       end
f (generic function with 1 method)

julia> my_parameter(x) = x
my_parameter (generic function with 1 method)

julia> ode = ODEProblem(f, ones(2), (0.0, 1.0), my_parameter)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 1.0)
u0: [1.0, 1.0]

julia> modelingtoolkitize(ode)
ERROR: MethodError: no method matching keys(::typeof(my_parameter))
Closest candidates are:
  keys(::Union{Tables.AbstractColumns, Tables.AbstractRow}) at ~/.julia/packages/Tables/UxLRG/src/Tables.jl:179
  keys(::IndexStyle, ::AbstractArray, ::AbstractArray...) at abstractarray.jl:327
  keys(::IterativeSolvers.ConvergenceHistory) at ~/.julia/packages/IterativeSolvers/TpeDx/src/history.jl:223
  ...
Stacktrace:
 [1] eachindex(itrs::Function)
   @ Base ./abstractarray.jl:253
 [2] modelingtoolkitize(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, typeof(my_parameter), ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/7KM6m/src/systems/diffeqs/modelingtoolkitize.jl:21
 [3] top-level scope
   @ REPL[5]:1

Using the same approach as suggested for MyParameter

julia> Base.eachindex(p::typeof(my_parameter)) = Base.OneTo(1)

julia> Base.size(p::typeof(my_parameter)) = (1,)

julia> Base.length(p::typeof(my_parameter)) = 1

julia> Base.iterate(p::typeof(my_parameter), x=1) = Base.iterate([p], x)

julia> f(du, u, p::Array, t) = f(du, u, p[1], t)
f (generic function with 2 methods)

julia> modelingtoolkitize(ode)
ERROR: Sym α₁ is not callable. Use @syms α₁(var1, var2,...) to create it as a callable. See ?@fun for more options
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:33
  [2] (::Sym{Real, Nothing})(args::Term{Real, Nothing})
    @ SymbolicUtils ~/.julia/packages/SymbolicUtils/p3pNS/src/types.jl:199
  [3] (::Num)(args::Num)
    @ Symbolics ~/.julia/packages/Symbolics/vetEm/src/num.jl:13
  [4] _broadcast_getindex_evalf
    @ ./broadcast.jl:648 [inlined]
  [5] _broadcast_getindex
    @ ./broadcast.jl:621 [inlined]
  [6] _getindex
    @ ./broadcast.jl:645 [inlined]
  [7] _broadcast_getindex
    @ ./broadcast.jl:620 [inlined]
  [8] getindex
    @ ./broadcast.jl:575 [inlined]
  [9] macro expansion
    @ ./broadcast.jl:984 [inlined]
 [10] macro expansion
    @ ./simdloop.jl:77 [inlined]
 [11] copyto!
    @ ./broadcast.jl:983 [inlined]
 [12] copyto!
    @ ./broadcast.jl:936 [inlined]
 [13] materialize!
    @ ./broadcast.jl:894 [inlined]
 [14] materialize!
    @ ./broadcast.jl:891 [inlined]
 [15] f(du::Vector{Num}, u::Vector{Num}, p::Num, t::Num)
    @ Main ./REPL[2]:2
 [16] f(du::Vector{Num}, u::Vector{Num}, p::Vector{Num}, t::Num)
    @ Main ./REPL[15]:1
 [17] (::ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing})(::Vector{Num}, ::Vararg{Any, N} where N)
    @ SciMLBase ~/.julia/packages/SciMLBase/070vZ/src/scimlfunctions.jl:329
 [18] modelingtoolkitize(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, typeof(my_parameter), ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/7KM6m/src/systems/diffeqs/modelingtoolkitize.jl:43
 [19] top-level scope
    @ REPL[16]:1

Trying to ignore the function parameter by setting something like

julia> Base.eachindex(p::typeof(my_parameter)) = Base.OneTo(0)

julia> Base.size(p::typeof(my_parameter)) = (0,)

julia> Base.length(p::typeof(my_parameter)) = 0

julia> Base.iterate(::typeof(my_parameter), x=1) = nothing

will not really help either, I think, since my_parameter will not be passed to f.

[ChrisRackauckas]

I think that needs a better form of tracing, like https://github.com/MasonProtter/SymbolicTracing.jl, and the function p would need to be @registered.