SciML / DataInterpolations.jl

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

Cache parameters #274

Closed SouthEndMusic closed 5 days ago

SouthEndMusic commented 2 weeks ago

Fixes https://github.com/SciML/DataInterpolations.jl/issues/271 (the other half)

Checklist

Additional context

Add any other context about the problem here.

SouthEndMusic commented 2 weeks ago

That was a bit fiddly to get right (or at least pass the LinearInterpolation tests). Quick benchmark:

using DataInterpolations
using BenchmarkTools

N = 1000
u = rand(N)
t = cumsum(rand(N))

itp = LinearInterpolation(u, t)

t_eval = range(first(t), last(t), length = 10000)

@btime itp($t_eval)

# old: 162.200 μs (3 allocations: 78.23 KiB)
# new: 119.900 μs (3 allocations: 78.23 KiB)
SouthEndMusic commented 1 week ago

Another benchmark:

using DataInterpolations
using BenchmarkTools

N = 1000
u = rand(N)
t = cumsum(rand(N))

itp = QuadraticInterpolation(u, t)

@btime DataInterpolations.integral(itp, t[end])

# old: 18.400 μs (3 allocations: 48 bytes)
# new:  9.000 μs (3 allocations: 48 bytes)
SouthEndMusic commented 1 week ago

I'll leave it at caching parameters for LinearInterpolation, QuadraticInterpolation, QuadraticSplineInterpolation and CubicSplineInterpolation for this PR.

Interesting ones to cache parameters for next are BsplineInterpolation and BsplineApprox, as they now allocate at every call in spline_coefficients:

using DataInterpolations
using BenchmarkTools

N = 1000
u = rand(N)
t = cumsum(rand(N))

t = [0, 62.25, 109.66, 162.66, 205.8, 252.3]
u = [14.7, 11.51, 10.41, 14.95, 12.24, 11.22]
A1 = QuadraticInterpolation(u,t)
A2 = BSplineInterpolation(u, t, 2, :Uniform, :Uniform)
A3 = BSplineApprox(u,t, 2, 4, :Uniform, :Uniform)

ts = collect(range(first(t), last(t), length = 10000))

@btime A1.(ts)
# 91.100 μs (3 allocations: 78.25 KiB)

@btime A2.(ts)
# 394.500 μs (10003 allocations: 1.14 MiB)

@btime A3.(ts)
# 387.400 μs (10003 allocations: 1015.77 KiB)

Edit: I made the non-ForwardDiff calls of BSplineInterpolation and BSplineApprox non-allocating by computing the B-spline coefficients in preallocated memory. This gives problems with ForwardDiff as expected, which could be solved using PreallocationTools, but I'm not sure that's desirable.

@btime A2.(ts)
# 250.600 μs (3 allocations: 78.28 KiB)

@btime A3.(ts)
# 244.500 μs (3 allocations: 78.28 KiB)
SouthEndMusic commented 1 week ago

Refactor of LagrangeInterpolation to be non-allocating and using idx_prev:

using DataInterpolations
using BenchmarkTools

N = 10
t = cumsum(rand(N))
u = rand(N)
itp = LagrangeInterpolation(u, t, N-1)

ts = range(first(t), last(t), length = 10000)

@btime itp($ts)
# old: 1.806 ms (20003 allocations: 2.82 MiB)
# new: 465.300 μs (3 allocations: 78.23 KiB)
SouthEndMusic commented 1 week ago

Using knowledge of which spline coefficients are non-zero:

using DataInterpolations
using BenchmarkTools

t = cumsum(rand(100))
u = rand(100)
A = BSplineInterpolation(u, t, 2, :Uniform, :Uniform)

ts = collect(range(first(t), last(t), length = 10000))

@btime A($ts)
# old: 1.161 ms (2 allocations: 78.17 KiB)
# new: 663.900 μs (2 allocations: 78.17 KiB)
SouthEndMusic commented 1 week ago

@sathvikbhagavan I made it so that push! and append! throw an error when safetycopy = true, but the error I now get in online_test.jl I think illustrates my point from earlier. u1, t1 are modified in the first iteration of the loop, which causes problems in the later iterations. Or would you say that this is bad use of the interpolation objects?

If this is indeed the behavior you want, I still need some help explaining the use of safetycopy in the docs.