Evalutate multiple interpolations/derivatives simultaneously

[bclyons12]

Is your feature request related to a problem? Please describe.

In a performance critical part of my code, I'm evaluating multiple interpolations at the same t location. Internally, this leads to lots of repeated calculations, since most of time in _interpolate is spent on lines that only depend on t and A.t. Here's a VScode profiling, with the blue bars showing the time spent on each line.

Describe the solution you’d like

It would be good to have some standard method for evaluating multiple interpolations with the same A.t grid at the same t location simultaneously so that the most computationally intensive lines are not unnecessarily recomputed. Likewise I am taking derivatives at the same point, which again reuses things like get_idx(), so having this standard method work with derivatives would be great as well.

Describe alternatives you’ve considered

• I could hack _interpolate myself, rewriting it for my specific case, but something more general would be appreciated.
• I could also use the FiniteElementHermite.jl package that I wrote a while back, which does exactly what I'd want: First evaulate the value of common basis functions, then evaluate the specific interpolation with those bases: https://github.com/ProjectTorreyPines/FiniteElementHermite.jl/blob/0677b81c29240cb7824bd66476b424dbe4374667/src/hermite.jl#L345-L355

[ChrisRackauckas]

with the same A.t grid

The interpolation can be an AbstractVector{T}, is there a reason to not just make these all one interpolation?

[bclyons12]

Maybe there's a way I could refactor it that way, but I see a couple potential problems:

1. I don't always want to compute all of these interpolations every time I evaluate one of them
2. In this context, I want the derivative of some quantities and the values of others.
3. I think if I combined them into a Matrix to evaluate them together, the result would allocate a small vector object every time, which wouldn't work for performance critical code.

[ChrisRackauckas]

There's nothing necessarily bad about adding an interface for this. But someone would have to take up the work to refactor all of the kernels which isn't small.

[cgarling]

with the same A.t grid

The interpolation can be an AbstractVector{T}, is there a reason to not just make these all one interpolation?

I have a similar use case where I'm evaluating multiple interpolations at the same t. Your suggestion of using one multi-valued interpolation works for some interpolations (e.g., CubicSpline), but for my use case I would rather use AkimaInterpolation or PCHIPInterpolation, neither of which seems to support the same type of multi-valued interpolation. Would you consider this a bug in these implementations, or is it expected that they do not support this behavior? The problem seems to be the use of methods (sign and abs) that do not have implementations for SVector inputs. The below code shows an example of this.

import StaticArrays: SVector
import DataInterpolations: CubicSpline, PCHIPInterpolation, AkimaInterpolation

x = 1:1_000
a,b,c = [rand(length(x)) for i in 1:3]
y = [SVector{3, Float64}(a[i], b[i], c[i]) for i in eachindex(a,b,c)]

cubic = CubicSpline(y, x)
cubic(2.5) # works, returns an SVector{3, Float64}
pchip = PCHIPInterpolation(y, x)
# pchip fails with ERROR: MethodError: no method matching sign(::SVector{3, Float64})
akima = AkimaInterpolation(y, x)
# akima fails with ERROR: MethodError: no method matching abs(::SVector{3, Float64})

[SouthEndMusic]

By the current design of DataInterpolations.jl, vector valued interpolations are always allocating.

I'd like to do a shameless plug of NDInterpolations.jl here, which is in the registration waiting period at the time of writing. This package can do in-place vector valued interpolation, and also caches basis function values. It currently assumes that you want to evaluate your interpolation at the same t values each time though. I think it's worth a look to see if this can be helpful to you.

[ChrisRackauckas]

The optimal thing here would be to do the same thing as the DiffEq interface, which is sol(out,t), where this can be used if the underlying elements are an array or if t is an array. And in the case where t is an array, it can pre-sort and then use the sorting assumption in order to decrease the amount of searching done https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/lib/OrdinaryDiffEqCore/src/dense/generic_dense.jl#L614-L619. We just never moved this part over to DataInterpolations.jl I think.