# Spline fit from DataInterpolations erroring out with ModelingToolkit symbolic derivative

**URL:** <https://discourse.julialang.org/t/spline-fit-from-datainterpolations-erroring-out-with-modelingtoolkit-symbolic-derivative/95417>\
**Category:** Modelling & Simulations\
**Tags:** question, sciml, modelingtoolkit\
**Created:** [March 1, 2023, 11:49pm UTC](https://discourse.julialang.org/t/spline-fit-from-datainterpolations-erroring-out-with-modelingtoolkit-symbolic-derivative/95417 "2023-03-01T23:49:50Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![aditya-sengupta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aditya-sengupta/32/47348_2.png) [@aditya-sengupta](https://discourse.julialang.org/u/aditya-sengupta)\
**Post date:** [March 1, 2023, 11:49pm UTC](https://discourse.julialang.org/t/spline-fit-from-datainterpolations-erroring-out-with-modelingtoolkit-symbolic-derivative/95417/1 "2023-03-01T23:49:50Z")

</div>

I’m trying to take a ModelingToolkit symbolic derivative of a variable whose values I’m getting with a DataInterpolations spline, and it’s unable to evaluate the derivative of a constant so I can’t solve the system. MWE:

```julia

using ModelingToolkit
using OrdinaryDiffEq
using DataInterpolations
using DataDrivenDiffEq

itp_method = InterpolationMethod(CubicSpline)
xr = 1.0:5.0
@register spl(z)
# originally had just spl = CubicSpline(xr,xr), with the same results
spl(z) = itp_method(xr, xr)(z)

@variables z
@variables v_terminal(z)
Dz = Differential(z)
@variables coll_rate(z)

@named system = ODESystem(
    [
        v_terminal ~ spl(z),
        Dz(coll_rate) ~ Dz(v_terminal)
        # this example is contrived, in the real system I do need Dz(v_terminal)
        # and anything simpler than this gets transformed away by structural_simplify
    ],
    z, [v_terminal, coll_rate], [],
)
system = structural_simplify(system)
prob = ODEProblem(system, [coll_rate => 1.0], (0.0, 1.0))
prob.f(prob.u0, prob.p, 0.0) 
# gives 1-element Vector{Term{Float64, Nothing}}: Differential(z)(0.0)
solve(prob, Tsit5()) # errors out because it can't convert Term{Float64,Nothing} to Float64

```

I attempted to use the wrapper from [DataDrivenDiffEq.jl](https://docs.sciml.ai/DataDrivenDiffEq/dev/utils/) but it doesn’t seem to have done anything. I tried to fix this example by adding in a variable like dv\_dz that I explicitly describe with `derivative(spl, z)`, but that throws `MethodError: isless(::Float64, ::Num) is ambiguous.` In case that does work, I still have a more complex function for v\_terminal whose derivative I need, and I don’t want to have to take its derivative by hand to plug it into the symbolic derivative terms. Is there a better way to handle this situation?

---

<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:** [March 2, 2023, 10:01pm UTC](https://discourse.julialang.org/t/spline-fit-from-datainterpolations-erroring-out-with-modelingtoolkit-symbolic-derivative/95417/2 "2023-03-02T22:01:30Z")

</div>

This looks tricky. I though I had a solution, but it fails in a place I don’t know how to fix with `ERROR: Differentiation with array expressions is not yet supported`.

Not defining a new function and instead using the `CubicInterpolation` struct gets a little farther. This allows a symbolic derivative to be registered and the derivative function registered as opaque as well. Unfortunately, `<:AbstractArray` means something different to Symbolics than it does to DataInterpolations, so I’m to sure where to go next.

```julia
using ModelingToolkit
using OrdinaryDiffEq
using DataInterpolations
using Symbolics

xr = 1.0:5.0
spl = CubicSpline(xr,xr)

# Define a symbolic derivative of the cubic spline
Symbolics.derivative(s::typeof(spl), args::NTuple{1, Any}, ::Val{1}) = DataInterpolations.derivative(s, args[1])

# register the derivative function so Symbolics does not attempt to trace it.
@register_symbolic DataInterpolations.derivative(itp::DataInterpolations.AbstractInterpolation, val)

@variables z
@variables v_terminal(z)
Dz = Differential(z)
@variables coll_rate(z)

# Now differentiating works...
expand_derivatives(Dz(spl(z)))
# DataInterpolations.derivative([1.0, 2.0, 3.0, 4.0, 5.0, 1.0, 2.0, 3.0, 4.0, 5.0], z)
substitute(expand_derivatives(Dz(spl(z))), z=>1.0)
# 1.0

@named system = ODESystem(
    [
        v_terminal ~ spl(z),
        Dz(coll_rate) ~ Dz(v_terminal)
        # this example is contrived, in the real system I do need Dz(v_terminal)
        # and anything simpler than this gets transformed away by structural_simplify
    ],
    z, [v_terminal, coll_rate], [],
)

# But now errors here because CubicSpline <: AbstractArray
system = structural_simplify(system)

prob = ODEProblem(system, [coll_rate => 1.0], (0.0, 1.0))

prob.f(prob.u0, prob.p, 0.0)

solve(prob, Tsit5())

```

> **The new error**
>
> julia\> system = structural\_simplify(system)  
> ERROR: Differentiation with array expressions is not yet supported  
> Stacktrace:  
> [1] error(s::String)  
> @ Base ./error.jl:35  
> [2] occursin\_info  
> @ ~/.julia/packages/Symbolics/UrqtQ/src/diff.jl:59 [inlined]  
> [3] (::Symbolics.var"#210#212"{Term{Real, Nothing}, Term{Real, Nothing}})(a::CubicSpline{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Float64}, Vector{Float64}, true, Float64})  
> @ Symbolics ~/.julia/packages/Symbolics/UrqtQ/src/diff.jl:92  
> [4] iterate  
> @ ./generator.jl:47 [inlined]  
> [5] \_collect(c::Vector{Any}, itr::Base.Generator{Vector{Any}, Symbolics.var"#210#212"{Term{Real, Nothing}, Term{Real, Nothing}}}, #unused#::Base.EltypeUnknown, isz::Base.HasShape{1})  
> @ Base ./array.jl:807  
> [6] collect\_similar(cont::Vector{Any}, itr::Base.Generator{Vector{Any}, Symbolics.var"#210#212"{Term{Real, Nothing}, Term{Real, Nothing}}})  
> @ Base ./array.jl:716  
> [7] map(f::Function, A::Vector{Any})  
> @ Base ./abstractarray.jl:2933  
> [8] occursin\_info(x::Term{Real, Nothing}, expr::Term{Real, Nothing}, fail::Bool)  
> @ Symbolics ~/.julia/packages/Symbolics/UrqtQ/src/diff.jl:92  
> [9] (::Symbolics.var"#210#212"{Term{Real, Nothing}, SymbolicUtils.Add{Real, Int64, Dict{Any, Number}, Nothing}})(a::Term{Real, Nothing})  
> @ Symbolics ~/.julia/packages/Symbolics/UrqtQ/src/diff.jl:92  
> [10] iterate  
> @ ./generator.jl:47 [inlined]  
> [11] collect\_to!(dest::Vector{Term{Real, Nothing}}, itr::Base.Generator{Vector{SymbolicUtils.Symbolic{Real}}, Symbolics.var"#210#212"{Term{Real, Nothing}, SymbolicUtils.Add{Real, Int64, Dict{Any, Number}, Nothing}}}, offs::Int64, st::Int64)  
> @ Base ./array.jl:845  
> [12] collect\_to\_with\_first!(dest::Vector{Term{Real, Nothing}}, v1::Term{Real, Nothing}, itr::Base.Generator{Vector{SymbolicUtils.Symbolic{Real}}, Symbolics.var"#210#212"{Term{Real, Nothing}, SymbolicUtils.Add{Real, Int64, Dict{Any, Number}, Nothing}}}, st::Int64)  
> @ Base ./array.jl:823  
> [13] \_collect(c::Vector{SymbolicUtils.Symbolic{Real}}, itr::Base.Generator{Vector{SymbolicUtils.Symbolic{Real}}, Symbolics.var"#210#212"{Term{Real, Nothing}, SymbolicUtils.Add{Real, Int64, Dict{Any, Number}, Nothing}}}, #unused#::Base.EltypeUnknown, isz::Base.HasShape{1})  
> @ Base ./array.jl:817  
> [14] collect\_similar(cont::Vector{SymbolicUtils.Symbolic{Real}}, itr::Base.Generator{Vector{SymbolicUtils.Symbolic{Real}}, Symbolics.var"#210#212"{Term{Real, Nothing}, SymbolicUtils.Add{Real, Int64, Dict{Any, Number}, Nothing}}})  
> @ Base ./array.jl:716  
> [15] map(f::Function, A::Vector{SymbolicUtils.Symbolic{Real}})  
> @ Base ./abstractarray.jl:2933  
> [16] occursin\_info(x::Term{Real, Nothing}, expr::SymbolicUtils.Add{Real, Int64, Dict{Any, Number}, Nothing}, fail::Bool)  
> @ Symbolics ~/.julia/packages/Symbolics/UrqtQ/src/diff.jl:92  
> [17] occursin\_info(x::Term{Real, Nothing}, expr::SymbolicUtils.Add{Real, Int64, Dict{Any, Number}, Nothing})  
> @ Symbolics ~/.julia/packages/Symbolics/UrqtQ/src/diff.jl:57  
> [18] expand\_derivatives(O::Term{Real, Nothing}, simplify::Bool; occurances::Nothing)  
> @ Symbolics ~/.julia/packages/Symbolics/UrqtQ/src/diff.jl:169  
> [19] expand\_derivatives(O::Term{Real, Nothing}, simplify::Bool)  
> @ Symbolics ~/.julia/packages/Symbolics/UrqtQ/src/diff.jl:163  
> [20] jacobian(ops::Vector{SymbolicUtils.Add{Real, Int64, Dict{Any, Number}, Nothing}}, vars::Vector{Any}; simplify::Bool)  
> @ Symbolics ~/.julia/packages/Symbolics/UrqtQ/src/diff.jl:443  
> [21] jacobian(ops::Vector{SymbolicUtils.Add{Real, Int64, Dict{Any, Number}, Nothing}}, vars::Vector{Any})  
> @ Symbolics ~/.julia/packages/Symbolics/UrqtQ/src/diff.jl:440  
> [22] (::ModelingToolkit.StructuralTransformations.var"#141#144"{TearingState{ODESystem}})(eqs::Vector{Int64}, vars::Vector{Int64})  
> @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/jCQlF/src/structural\_transformation/symbolics\_tearing.jl:729  
> [23] dummy\_derivative\_graph!(structure::SystemStructure, var\_eq\_matching::ModelingToolkit.BipartiteGraphs.Matching{ModelingToolkit.BipartiteGraphs.Unassigned, Vector{Union{ModelingToolkit.BipartiteGraphs.Unassigned, Int64}}}, jac::ModelingToolkit.StructuralTransformations.var"#141#144"{TearingState{ODESystem}}, ::Tuple{ModelingToolkit.AliasGraph, Nothing}, state\_priority::Function)  
> @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/jCQlF/src/structural\_transformation/partial\_state\_selection.jl:220  
> [24] dummy\_derivative\_graph!(state::TearingState{ODESystem}, jac::Function, ::Tuple{ModelingToolkit.AliasGraph, Nothing}; state\_priority::Function, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})  
> @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/jCQlF/src/structural\_transformation/partial\_state\_selection.jl:159  
> [25] dummy\_derivative(sys::ODESystem, state::TearingState{ODESystem}, ag::ModelingToolkit.AliasGraph; simplify::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})  
> @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/jCQlF/src/structural\_transformation/symbolics\_tearing.jl:747  
> [26] \_structural\_simplify!(state::TearingState{ODESystem}, io::Nothing; simplify::Bool, check\_consistency::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})  
> @ ModelingToolkit.SystemStructures ~/.julia/packages/ModelingToolkit/jCQlF/src/systems/systemstructure.jl:542  
> [27] structural\_simplify!(state::TearingState{ODESystem}, io::Nothing; simplify::Bool, check\_consistency::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})  
> @ ModelingToolkit.SystemStructures ~/.julia/packages/ModelingToolkit/jCQlF/src/systems/systemstructure.jl:496  
> [28] structural\_simplify(sys::ODESystem, io::Nothing; simplify::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})  
> @ ModelingToolkit ~/.julia/packages/ModelingToolkit/jCQlF/src/systems/systems.jl:39  
> [29] structural\_simplify (repeats 2 times)  
> @ ~/.julia/packages/ModelingToolkit/jCQlF/src/systems/systems.jl:19 [inlined]  
> [30] top-level scope  
> @ REPL[51]:1

---

<div class="post-metadata">

**Author:** ![aditya-sengupta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aditya-sengupta/32/47348_2.png) [@aditya-sengupta](https://discourse.julialang.org/u/aditya-sengupta)\
**Post date:** [March 2, 2023, 10:19pm UTC](https://discourse.julialang.org/t/spline-fit-from-datainterpolations-erroring-out-with-modelingtoolkit-symbolic-derivative/95417/3 "2023-03-02T22:19:38Z")

</div>

Maybe wrapping the interpolation in a new struct that isn’t \<: AbstractArray would help? It wouldn’t be very elegant but if it works I don’t mind. I can try this out in a bit

---

<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:** [March 2, 2023, 10:32pm UTC](https://discourse.julialang.org/t/spline-fit-from-datainterpolations-erroring-out-with-modelingtoolkit-symbolic-derivative/95417/4 "2023-03-02T22:32:30Z")

</div>

Unfortunately that doesn’t work yet either, there is a [PR](https://github.com/JuliaSymbolics/Symbolics.jl/issues/806) though.

---

<div class="post-metadata">

**Author:** ![aditya-sengupta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aditya-sengupta/32/47348_2.png) [@aditya-sengupta](https://discourse.julialang.org/u/aditya-sengupta)\
**Post date:** [March 4, 2023, 6:40pm UTC](https://discourse.julialang.org/t/spline-fit-from-datainterpolations-erroring-out-with-modelingtoolkit-symbolic-derivative/95417/5 "2023-03-04T18:40:38Z")

</div>

For now I’ve been able to work around this by explicitly constructing the derivative as its own spline:

```julia
using ModelingToolkit
using OrdinaryDiffEq
using DataInterpolations
using Symbolics

xr = 1.0:5.0
spl = CubicSpline(xr,xr)
spl_d = QuadraticSpline(map(x -> DataInterpolations.derivative(spl, x), xr), xr)

@variables z
@variables v_terminal(z)
Dz = Differential(z)
@variables coll_rate(z)

@named system = ODESystem(
    [
        Dz(coll_rate) ~ spl_d(z)
    ],
    z, [coll_rate], [],
)

system = structural_simplify(system)
prob = ODEProblem(system, [coll_rate => 1.0], (0.0, 1.0))
prob.f(prob.u0, prob.p, 0.0)
solve(prob, Tsit5()) # runs as expected, gives coll_rate(t) = t + 1

```

This makes it difficult to describe complicated functions of splines like my main use case, but it should at least be possible by taking a lot of derivatives by hand.