Altering the type of interpolation used with SampledData blocks

[aditya-sengupta]

I'd like to use any of the DataInterpolations.jl interpolation types with SampledData blocks. In particular I'm interested in the cubic spline interpolation, but anything past linear would be useful. I think this would involve changing this call to linear_interpolation to something from DataInterpolations.jl or to one of a few new functions.

[aditya-sengupta]

Possibly unrelated points, but I wasn't sure whether they were worth opening a new issue:

• do SampledData blocks support non-uniform sampling? As far as I can tell they just take a single sample_time instead of a vector of them.
• is it possible to use pre-transformed derived functions (that have been @register_symbolicd) on sampled data at the system creation level, so that structural_simplify will automatically apply it and correctly update the sampling/interpolation for the derived data? (This is probably possible to do manually, it'd just be nice to offload + relatively easy to automate.)

[ChrisRackauckas]

Yes that effectively would just be allowing the user to have the option to pass a different interpolation function.

do SampledData blocks support non-uniform sampling? As far as I can tell they just take a single sample_time instead of a vector of them.

They could, but the ODE solver only ever calls them at a single time point at a time so I'm not sure what the optimization would be here?

is it possible to use pre-transformed derived functions (that have been @register_symbolicd) on sampled data at the system creation level, so that structural_simplify will automatically apply it and correctly update the sampling/interpolation for the derived data? (This is probably possible to do manually, it'd just be nice to offload + relatively easy to automate.)

I'm not sure if it's implemented, but doing so should effectively work just by using observed.

[aditya-sengupta]

They could, but the ODE solver only ever calls them at a single time point at a time so I'm not sure what the optimization would be here?

Makes sense, I looked more into how it works and it seems like it'd be fine if I just resampled my data before passing it into the solver.

I'm not sure if it's implemented, but doing so should effectively work just by using observed.

Could you elaborate a bit on how this would work?

[ChrisRackauckas]

Could you elaborate a bit on how this would work?

Actually, rereading what you wrote, now I'm confused. Can you try describing it another way?

[aditya-sengupta]

Sure - I have sampled data that I'd like to feed into source terms for an ODE. The (data) -> (source term) transformation is a function that I can @register_symbolic, and I'd like my ODESystem to automatically apply that transformation and interpolate accurately on it to get the source terms. My concerns are avoiding repeated recomputation of the data-to-source transformation (i.e. making sure structural_simplify transforms that away) and making sure the interpolation on the source term is reasonably within the range of the transformation.

Currently, I can just apply the transformation myself and pass that in as a SampledData object, but that adds more steps before I'm able to remake the system each time, and it feels like the kind of thing MTK could automate.

[bradcarman]

If I understand correctly, what you'd like to do is shown here: SampleData demo. You can see that the get_prob() function is using remake to update the data given to the system without having to re-run structural_simplify(). And yes, SampledData component only supports sampled data computing a linear interpolation in time. But this could be improved, if you can post a MWE to work off of that would be really helpful.

[aditya-sengupta]

I changed the example you linked a bit to generate my MWE.

using ModelingToolkit
using ModelingToolkitStandardLibrary.Blocks: SampledData, Parameter
using OrdinaryDiffEq
using Plots
using DataInterpolations

@parameters t
D = Differential(t)

# Imagine we have a physical process/another simulation with few available data points
# but an expectation that the data should be ~smoothly varying.
# The "data" fields here are being used as a forcing function, as a stand-in for an actual science case.
times_coarse = 0:1e-2:0.1 # 11 data points
times_fine = 0:5e-5:0.1 # 2001 data points
data_coarse = sin.(2 * pi * times_coarse * 30)
data_fine = sin.(2 * pi * times_fine * 30)
# This is my simulated version of what I'd like MTK to do.
# This does a finer spline interpolation and MTK then linearly interpolates on top of that, which is probably fine
# but it'd be nicer if MTK could generate the spline and dynamically sample from it at whichever times it ends up needing.
data_resampled = (CubicSpline(data_coarse, times_coarse)).(times_fine) 

src = SampledData(Float64; name=:system)
vars = @variables f(t)=0 x(t)=0 dx(t)=0 ddx(t)=0
pars = @parameters mass=10 k=1000 d=1
eqs = [f ~ src.output.u
    ddx * 10 ~ k * x + d * dx + f
    D(x) ~ dx
    D(dx) ~ ddx]
system = ODESystem(eqs, t, vars, pars; systems = [src], name=:system)
sys = structural_simplify(system)
prob = ODEProblem(sys, [], (0, times_coarse[end]); tofloat = false)

function get_prob_and_solve(data, dt)
    solve(
        remake(prob, p = [prob.p[1], [Parameter(data, dt)]]),
        ImplicitEuler()
    )
end

sol_coarse = get_prob_and_solve(data_coarse, minimum(diff(times_coarse)))
sol_fine = get_prob_and_solve(data_fine, minimum(diff(times_fine)))
sol_resampled = get_prob_and_solve(data_resampled, minimum(diff(times_fine)))
# This plot shows the resampled solution getting closer to the fine one
# which demonstrates increased accuracy due to the resampling.
# This may increase slightly if the spline were used throughout within MTK.
begin
    plot(sol_coarse.t, sol_coarse[x], label="Coarse", xlabel="t", ylabel="x")
    plot!(sol_fine.t, sol_fine[x], label="Fine")
    plot!(sol_resampled.t, sol_resampled[x], label="Coarse data with spline interpolation")
end