Cannot evaluate spline parameter directly in an ODESystem

[hersle]

I want to evaluate a spline and its derivative(s) through an ODESystem. This works beautifully:

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

ts = 0.0:0.1:1.0
Fspline_t² = CubicSpline(ts .^ 2, ts) # spline for F(t) = t²

@variables F(t) F′(t)
@named sys = ODESystem([
    F ~ Fspline_t²(t)
    F′ ~ D(F)
], t)
ssys = structural_simplify(sys)

prob = ODEProblem(ssys, [], (0.0, 1.0), [])
sol = solve(prob)
@test sol(0.5, idxs=F′) ≈ 2 * 0.5 # F′(t) = 2t

However, to pass different splines, I want to make it a parameter:

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

ts = 0.0:0.1:1.0
Fspline_t² = CubicSpline(ts .^ 2, ts) # spline for F(t) = t²

@parameters Fspline::CubicSpline = Fspline_t²
@variables F(t) F′(t)
@named sys = ODESystem([
    F ~ Fspline(t)
    F′ ~ D(F)
], t)
ssys = structural_simplify(sys)

prob = ODEProblem(ssys, [], (0.0, 1.0), [])
sol = solve(prob)
@test sol(0.5, idxs=F′) ≈ 2 * 0.5 # F′(t) = 2t

This fails with ERROR: LoadError: Sym Fspline is not callable. Use @syms Fspline(var1, var2,...) to create it as a callable. Accordingly, I tried to declare @parameters Fspline(t)::CubicSpline = Fspline_t² instead, but it does not help. Is there a bug here?

I have worked around half the problem by evaluating the spline like F ~ spleval(x, Fspline) through a proxy function

@register_symbolic spleval(x, spline::CubicSpline)
spleval(x, spline) = spline(x)

This works for F. But it does not work for F′, which then evaluates to e.g. Differential(t)(1.0). Also, it feels unnecessarily complicated, because I believe these derivative rules are already defined in DataInterpolations.jl.

I really wish that the second example worked like the first, without the complicating workarounds. Is that possible?

Thanks!

[ChrisRackauckas]

I think a few of these are effectively the same issue that we need to make p::LinearInterpolation(t) work

[AayushSabharwal]

Roadmap for supporting this is in #2910

[hersle]

Thanks! For future reference, after #2995, evaluating the spline directly works great:

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

ts = 0.0:0.1:1.0
Fspline_t² = CubicSpline(ts .^ 2, ts) # spline for F(t) = t²

@parameters (Fspline::CubicSpline)(..)
@variables F(t) F′(t)
@named sys = ODESystem([
    F ~ Fspline(t)
    #F′ ~ D(F) # TODO: want to evaluate spline derivative
    #F′ ~ ModelingToolkit.derivative(Fspline, t) # TODO: want to evaluate spline derivative
], t)
ssys = structural_simplify(sys)

prob = ODEProblem(ssys, [], (0.0, 1.0), [Fspline => Fspline_t²])
sol = solve(prob)
@test sol(0.5, idxs=F) ≈ 0.5^2
#@test sol(0.5, idxs=F′) ≈ 2 * 0.5 # TODO: want to evaluate spline derivative

Is there also a way to evaluate the derivative of the spline, perhaps in a way that resembles any of the 3 commented lines?

[ChrisRackauckas]

DataInterpolations.derivative(Fspline, t) should work?

[hersle]

Not quite, unfortunately. Doing

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

ts = 0.0:0.1:1.0
Fspline_t² = CubicSpline(ts .^ 2, ts) # spline for F(t) = t²

@parameters (Fspline::CubicSpline)(..)
@variables F(t) F′(t)
@named sys = ODESystem([
    F ~ Fspline(t)
    F′ ~ DataInterpolations.derivative(Fspline, t)
], t)
ssys = structural_simplify(sys)

prob = ODEProblem(ssys, [], (0.0, 1.0), [Fspline => Fspline_t²])
sol = solve(prob)
@test sol(0.5, idxs=F) ≈ 0.5^2
@test sol(0.5, idxs=F′) ≈ 2 * 0.5

fails with

ERROR: LoadError: type CallWithMetadata has no field t
Stacktrace:
 [1] getproperty
   @ .\Base.jl:37 [inlined]
 [2] derivative(A::Symbolics.CallWithMetadata{SymbolicUtils.FnType{…}, Base.ImmutableDict{…}}, t::Num, order::Int64)
   @ DataInterpolations C:\Users\herma\.julia\packages\DataInterpolations\iUe8V\src\derivatives.jl:2
 [3] top-level scope
   @ C:\Users\herma\.julia\packages\ModelingToolkit\Vsl3C\src\systems\abstractsystem.jl:1987
in expression starting at C:\Users\herma\.julia\dev\SymBoltz\bug.jl:9 
Some type information was truncated. Use `show(err)` to see complete types.

Could the issue be that it is illegal to pass the argument-free Fspline to DataInterpolations.derivative, because it is defined as a callable parameter? I think this example shows it is necessary to allow using it both as Fspline and Fspline(t).

It seems to call the wrong derivative function instead of the desired Symbolics derivative extension?

[hersle]

Actually, this works 100%:

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

ts = 0.0:0.1:1.0
Fspline_t² = CubicSpline(ts .^ 2, ts) # spline for F(t) = t²

value(s, t) = s(t)
derivative(s, t) = DataInterpolations.derivative(s, t)
@register_symbolic value(s::CubicSpline, t)
@register_symbolic derivative(s::CubicSpline, t)

@parameters Fspline::CubicSpline
@variables F(t) F′(t)
@named sys = ODESystem([
    F ~ value(Fspline, t)
    F′ ~ derivative(Fspline, t)
], t)
ssys = structural_simplify(sys)

prob = ODEProblem(ssys, [], (0.0, 1.0), [Fspline => Fspline_t²])
sol = solve(prob)
@test sol(0.5, idxs=F) ≈ 0.5^2
@test sol(0.5, idxs=F′) ≈ 2*0.5

It feels like a workaround, though. I lose the callable spline, but register two proxy functions for evaluating its value and derivative. The latter must also be typed; just @register_symbolic derivative(s, t) does not work.

[AayushSabharwal]

It seems to call the wrong derivative function instead of the desired Symbolics derivative extension?

Yes, DataInterpolations.derivative and Symbolics.derivative are two different functions. I think codegen in general doesn't understand how to handle derivatives of callable symbolics with respect to one of their arguments. It is possible to symbolically represent calling a spline and taking its derivative:

julia> Fspline_t²(t)
CubicSpline{Vector{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}([0.0, 0.010000000000000002, 0.04000000000000001, 0.09, 0.16000000000000003, 0.25, 0.36, 0.48999999999999994, 0.6400000000000001, 0.81, 1.0], 0.0:0.1:1.0, Float64[], DataInterpolations.CubicSplineParameterCache{Vector{Float64}}(Float64[], Float64[]), [0.0, 0.1, 0.1, 0.09999999999999998, 0.10000000000000003, 0.09999999999999998, 0.09999999999999998, 0.09999999999999998, 0.10000000000000009, 0.09999999999999998, 0.09999999999999998], [0.0, 2.5359116022099464, 1.8563535911602187, 2.038674033149176, 1.9889502762430864, 2.0055248618784605, 1.988950276243084, 2.03867403314919, 1.856353591160198, 2.535911602209954, 0.0], false, FindFirstFunctions.Guesser{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}(0.0:0.1:1.0, Base.RefValue{Int64}(1), true), false, false)(t)

julia> D(Fspline_t²(t))
Differential(t)(CubicSpline{Vector{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}([0.0, 0.010000000000000002, 0.04000000000000001, 0.09, 0.16000000000000003, 0.25, 0.36, 0.48999999999999994, 0.6400000000000001, 0.81, 1.0], 0.0:0.1:1.0, Float64[], DataInterpolations.CubicSplineParameterCache{Vector{Float64}}(Float64[], Float64[]), [0.0, 0.1, 0.1, 0.09999999999999998, 0.10000000000000003, 0.09999999999999998, 0.09999999999999998, 0.09999999999999998, 0.10000000000000009, 0.09999999999999998, 0.09999999999999998], [0.0, 2.5359116022099464, 1.8563535911602187, 2.038674033149176, 1.9889502762430864, 2.0055248618784605, 1.988950276243084, 2.03867403314919, 1.856353591160198, 2.535911602209954, 0.0], false, FindFirstFunctions.Guesser{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}(0.0:0.1:1.0, Base.RefValue{Int64}(1), true), false, false)(t))

julia> D(Fspline_t²(t)) |> expand_derivatives
DataInterpolations.derivative(CubicSpline{Vector{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}([0.0, 0.010000000000000002, 0.04000000000000001, 0.09, 0.16000000000000003, 0.25, 0.36, 0.48999999999999994, 0.6400000000000001, 0.81, 1.0], 0.0:0.1:1.0, Float64[], DataInterpolations.CubicSplineParameterCache{Vector{Float64}}(Float64[], Float64[]), [0.0, 0.1, 0.1, 0.09999999999999998, 0.10000000000000003, 0.09999999999999998, 0.09999999999999998, 0.09999999999999998, 0.10000000000000009, 0.09999999999999998, 0.09999999999999998], [0.0, 2.5359116022099464, 1.8563535911602187, 2.038674033149176, 1.9889502762430864, 2.0055248618784605, 1.988950276243084, 2.03867403314919, 1.856353591160198, 2.535911602209954, 0.0], false, FindFirstFunctions.Guesser{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}(0.0:0.1:1.0, Base.RefValue{Int64}(1), true), false, false), t, 1)

But codegen doesn't understand that D(Fspline(t)) should turn into Symbolics.derivative(Fspline, t, 1). In fact, I'm not entirely sure how to get codegen to work since the only difference between Fspline(t) and F in the code above is their metadata.

[AayushSabharwal]

It's also a little annoying that this doesn't work:

julia> @register_symbolic DataInterpolations.derivative(interp::DataInterpolations.AbstractInterpolation, t)
julia> DataInterpolations.derivative(Fspline, t)
ERROR: type CallWithMetadata has no field t

Because the registered function expects a BasicSymbolic{AbstractInterpolation} not a CallWithMetadata

@ChrisRackauckas should we add a method to @register_symbolic which does this? Maybe then at least we can enable doing

F′ ~ DataInterpolations.derivative(Fspline, t)

[ChrisRackauckas]

Yeah seems like we're missing a dispatch here.

[AayushSabharwal]

Should all nonnumeric arguments to the registered function have an additional method where they're callable, or do we want an opt-in (Symbolics.is_callable_type or something) for which we generate the additional method, or different macro syntax?

[AayushSabharwal]

The first would potentially cause an explosion in the number of methods, the second is a little weird API imo. I like the third but don't have a desirable syntax off the top of my head.