SciML / DataInterpolations.jl

A library of data interpolation and smoothing functions
MIT License
203 stars 43 forks source link

Invert interpolation integrals #282

Closed SouthEndMusic closed 4 days ago

SouthEndMusic commented 1 week ago

Fixes https://github.com/SciML/DataInterpolations.jl/issues/277

Checklist

SouthEndMusic commented 1 week ago

Integral inversion POC:

using DataInterpolations
using Plots

t = collect(1:5)
u = [1.0, 1.0, 2.0, 4.0, 3.0]
A = LinearInterpolation(u, t)

ts = range(first(t), last(t), length = 100)
Is = DataInterpolations.integral.(Ref(A), ts)

A_intinv = DataInterpolations.invert_integral(A)
ts_ = A_intinv.(Is)

plot(ts, Is, label = "Original integral")
plot!(ts_, Is, label = "Integral inverse")

plot

It handles the edge case of the interpolation being constant perfectly 🥳

SouthEndMusic commented 6 days ago

Speed boost when using cached integral values:

using DataInterpolations
using BenchmarkTools
using Random
Random.seed!(42)

N = 100
t = cumsum(rand(N))
ddu = rand(N)
du = rand(N)
u = rand(N)
A = QuinticHermiteSpline(ddu, du, u, t)

@btime DataInterpolations.integral(A, last(t))

# old: 1.320 μs (2 allocations: 32 bytes)
# new: 66.395 ns (2 allocations: 32 bytes)
SouthEndMusic commented 6 days ago

@sathvikbhagavan how shall I document this new feature?

ChrisRackauckas commented 5 days ago

Make a new docs page on it with a short tutorial and then the docstrings. And make a list of the methods which are compatible with it on there.

SouthEndMusic commented 5 days ago

I'm wondering what to do with the docstrings of the LinearInterpolationIntInv and ConstantInterpolationIntInv structs. For the other interpolation methods this describes the public constructor, but I've set up my code so that users will probably always create these objects via DataInterpolations.invert_integral.