Cache parameters

[SouthEndMusic]

Fixes #271 (the other half)

Checklist

• Appropriate tests were added
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• Any new documentation only uses public API

Additional context

Add any other context about the problem here.

[SouthEndMusic]

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]

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]

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]

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]

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)

[sathvikbhagavan]

Looks good! Can you also add a few tests to ensure all cache parameters are set properly?

[sathvikbhagavan]

I think it would be better to explain why we need safetycopy

[sathvikbhagavan]

Shouldn't this be allowed only if safetycopy is false?

[SouthEndMusic]

I don't think I fully understand the reasoning here. If safetycopy is true, then the u, t vectors within the interpolation object can safely be assumed to be only used within the context of this package. Surely then it is safe to add data to it? Only if u, t are also used in code outside the DataInterpolations context do you run the risk of breaking that code by modifying u, t.

[sathvikbhagavan]

I think I understand your point. It comes down what should be the semantics of push! for interpolation objects. Currently it mutates the original array. Maybe we should handle both cases? If it is true, it would push to the copy, if it is false it would mutate the original array. For both we have to compute the parameters again.
@ChrisRackauckas, your thoughts?

[SouthEndMusic]

@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.

[sathvikbhagavan]

Can you also add a line describing which interpolations support caching?

[sathvikbhagavan]

Is this needed? Shouldn't we allow for both?

[SouthEndMusic]

That is fine by me, I think it is clear enough that this leads to mutation of the input u, t if someone explicitly sets safetycopy = false.

[sathvikbhagavan]

Why was this changed?

[SouthEndMusic]

I added this function to conveniently test whether the implementation of a new interpolation type is complete: https://github.com/SciML/DataInterpolations.jl/pull/274/files#diff-f66955aae2aee262495ff068999fe35f8501cd1cbef7b1057aa26d2f1099a4f3R6-R17

To do this, I moved throwing an error when integration is not defined from specialised methods of integral to methods of _integral. However, that now means that the extrapolation error is thrown before the 'integral not found' error. So, to still test the integral not found error, I had to make sure that the test for this does not try to extrapolate.

[sathvikbhagavan]

ok, that should be fine.

[sathvikbhagavan]

Almost there. Just a few cleanups.

[sathvikbhagavan]

ok, that should be fine.

[sathvikbhagavan]

LGTM! Nice work! Thanks!