Support for callable parameters in `IndexCache`

[SebastianM-C]

Would it be possible to support callable symbolic parameters? I tried the following

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using DataInterpolations

u = rand(5)
x = 0:4
interp = LinearInterpolation(u, x)

@variables y(t) = 0
@parameters i(..) = interp

eqs = [D(y) ~ i(t)]

@named model = ODESystem(eqs, t, [y], [i])
sys = structural_simplify(model)

which gives

ERROR: MethodError: Cannot `convert` an object of type 
  Symbolics.CallWithMetadata{SymbolicUtils.FnType{Tuple, Real}, Base.ImmutableDict{DataType, Any}} to an object of type 
  SymbolicUtils.BasicSymbolic
The function `convert` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  convert(::Type{T}, ::T) where T
   @ Base Base.jl:126

Stacktrace:
  [1] setindex!(h::Dict{…}, v0::Nothing, key0::Symbolics.CallWithMetadata{…})
    @ Base ./dict.jl:346
  [2] push!(s::Set{…}, x::Symbolics.CallWithMetadata{…})
    @ Base ./set.jl:137
  [3] (::ModelingToolkit.var"#insert_by_type!#94")(buffers::Dict{…}, sym::Symbolics.CallWithMetadata{…})
    @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/systems/index_cache.jl:103
  [4] ModelingToolkit.IndexCache(sys::ODESystem)
    @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/systems/index_cache.jl:245
  [5] macro expansion
    @ ~/.julia/packages/Setfield/PdKfV/src/sugar.jl:197 [inlined]
  [6] complete(sys::ODESystem; split::Bool)
    @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/systems/abstractsystem.jl:733
  [7] complete
    @ ~/.julia/dev/ModelingToolkit/src/systems/abstractsystem.jl:731 [inlined]
  [8] macro expansion
    @ ~/.julia/packages/Setfield/PdKfV/src/sugar.jl:197 [inlined]
  [9] structural_simplify(sys::ODESystem, io::Nothing; simplify::Bool, split::Bool, allow_symbolic::Bool, allow_parameter::Bool, conservative::Bool, fully_determined::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/systems/systems.jl:41
 [10] structural_simplify
    @ ~/.julia/dev/ModelingToolkit/src/systems/systems.jl:20 [inlined]
 [11] structural_simplify(sys::ODESystem)
    @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/systems/systems.jl:20
 [12] top-level scope
    @ ~/juliastuff/dev/parametrized_interpolation.jl:50
Some type information was truncated. Use `show(err)` to see complete types.

Should I specify the type for the callable parameter?

Also, if I don't specify the unknowns and parameters explicitly, I get

ERROR: MethodError: no method matching iterate(::LinearInterpolation{Vector{Float64}, UnitRange{Int64}, Vector{Float64}, Vector{Float64}, Float64})
The function `iterate` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  iterate(::Pkg.Registry.RegistryInstance)
   @ Pkg ~/.julia/juliaup/julia-1.11.0-rc1+0.x64.linux.gnu/share/julia/stdlib/v1.11/Pkg/src/Registry/registry_instance.jl:460
  iterate(::Pkg.Registry.RegistryInstance, ::Any)
   @ Pkg ~/.julia/juliaup/julia-1.11.0-rc1+0.x64.linux.gnu/share/julia/stdlib/v1.11/Pkg/src/Registry/registry_instance.jl:461
  iterate(::Tables.DictRowTable)
   @ Tables ~/.julia/packages/Tables/8p03y/src/dicts.jl:122
  ...

Stacktrace:
  [1] _foldl_impl(op::Base.BottomRF{…}, init::Set{…}, itr::LinearInterpolation{…})
    @ Base ./reduce.jl:56
  [2] foldl_impl(op::Base.BottomRF{…}, nt::Set{…}, itr::LinearInterpolation{…})
    @ Base ./reduce.jl:48
  [3] mapfoldl_impl(f::typeof(identity), op::ModelingToolkit.var"#54#55"{…}, nt::Set{…}, itr::LinearInterpolation{…})
    @ Base ./reduce.jl:44
  [4] mapfoldl(f::Function, op::Function, itr::LinearInterpolation{…}; init::Set{…})
    @ Base ./reduce.jl:175
  [5] foldl(op::Function, itr::LinearInterpolation{…}; kw::@Kwargs{…})
    @ Base ./reduce.jl:198
  [6] vars(exprs::LinearInterpolation{Vector{Float64}, UnitRange{Int64}, Vector{Float64}, Vector{Float64}, Float64}; op::Type)
    @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/utils.jl:374
  [7] collect_vars!(unknowns::OrderedCollections.OrderedSet{…}, parameters::OrderedCollections.OrderedSet{…}, expr::LinearInterpolation{…}, iv::SymbolicUtils.BasicSymbolic{…}; op::Type)
    @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/utils.jl:470
  [8] collect_vars!(unknowns::OrderedCollections.OrderedSet{…}, parameters::OrderedCollections.OrderedSet{…}, expr::LinearInterpolation{…}, iv::SymbolicUtils.BasicSymbolic{…})
    @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/utils.jl:466
  [9] collect_var!(unknowns::OrderedCollections.OrderedSet{…}, parameters::OrderedCollections.OrderedSet{…}, var::SymbolicUtils.BasicSymbolic{…}, iv::SymbolicUtils.BasicSymbolic{…})
    @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/utils.jl:489
 [10] collect_vars!(unknowns::OrderedCollections.OrderedSet{…}, parameters::OrderedCollections.OrderedSet{…}, expr::SymbolicUtils.BasicSymbolic{…}, iv::SymbolicUtils.BasicSymbolic{…}; op::Type)
    @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/utils.jl:474
 [11] collect_vars!(unknowns::OrderedCollections.OrderedSet{…}, parameters::OrderedCollections.OrderedSet{…}, expr::SymbolicUtils.BasicSymbolic{…}, iv::SymbolicUtils.BasicSymbolic{…})
    @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/utils.jl:466
 [12] ODESystem(eqs::Vector{Equation}, iv::Num; kwargs::@Kwargs{name::Symbol})
    @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/systems/diffeqs/odesystem.jl:301
 [13] top-level scope
    @ ~/.julia/dev/ModelingToolkit/src/systems/abstractsystem.jl:1810
Some type information was truncated. Use `show(err)` to see complete types.

[hersle]

See #2823, #2646, #1615 and #619 for context/discussion/workarounds.

[hersle]

This works:

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using DataInterpolations
using DifferentialEquations

u = rand(5)
x = 0:4
interp = CubicSpline(u, x)

@variables y(t) = 0
@parameters i::CubicSpline = interp
@register_symbolic ifunc(t, ispline::CubicSpline)
ifunc(t, ispline) = ispline(t)

eqs = [D(y) ~ ifunc(t, i)]

@named model = ODESystem(eqs, t, [y], [i])
sys = structural_simplify(model)
prob = ODEProblem(sys, [], extrema(x), [])
sol = solve(prob)

But (for reasons I do not know)

• it does not work with LinearInterpolation instead of CubicSpline
• the spline is evaluated through the proxy function ifunc(t, ispline), but I wish it could be done directly in a way that is similar to your example

[SebastianM-C]

I found a (better?) solution that should work without callable parameters in https://github.com/SciML/ModelingToolkitStandardLibrary.jl/pull/314

[hersle]

Beautiful! I would love to use this.

[AayushSabharwal]

Roadmap for callable parameters support:

1. Add an additional type parameter to SymbolicUtils.FnType to store the type of the function (FnType{Tuple{argtypes…}, rettype, functiontype}) and make it default to a sentinel value (Nothing?) if not provided (https://github.com/JuliaSymbolics/SymbolicUtils.jl/pull/648)
2. Change the @variables macro parsing to allow the syntax fn::functiontype(..)::rettype. .. may also be replaced by ::argtype1, argtype2.... Currently (fn::T)(..) is identical to fn(..)::T so the former can be made to do something useful, as it's arguably not syntax anyone would be using (https://github.com/JuliaSymbolics/Symbolics.jl/pull/1270)
3. Allow callable parameters in MTK, which are stored as type functiontype if provided, and a FunctionWrapper if not.