Faster/performant interpolation of very big arrays

General Usage Performance
1

I have two separate problems relating to interpolation of some dimensions of relatively big arrays.

  1. Creating interpolations objects using `Interpolations.jl’ is a little slow and uses a lot of memory when you create a lot of them.
  2. I can’t find out a way to re-use the interpolation weights instead of finding them every time. (frankly, the interpolations are so fast that this isn’t really a big bottleneck)

Hopefully this example illustrates:

## Small example
using Interpolations
using BenchmarkTools
xgrd = [1, 1.5, 2.5];nx=length(xgrd) #Not evenlye spaced
ygrd = [1, 1.5, 2.5];ny=length(ygrd)
Nd1 = 5
Nd2 = 5
A = [xgrd[ix]*ygrd[iy]+id1+id2 for id1=1:Nd1,id2=1:Nd2,iy=1:nx,ix=1:nx] 
    # Nd1*Nd2*3*3 array,  we will interpolate over the last two dimensions

# function to create interpolation object
function createitp(foo,xgrd,ygrd,Nd1,Nd2) 
    itp = [linear_interpolation((xgrd,ygrd,),foo[id1,id2,:,:],extrapolation_bc=Line()) 
        for id1=1:Nd1,id2=1:Nd2]
    # Maybe there is a way to do this faster/using less memory?
end

function evalute(fooitp,Nd1,Nd2,xval,yval)
    value = 0
    for id1=1:Nd1
        for id2=1:Nd2
            value += fooitp[id1,id2](xval,yval)
            # Presumably there is way to just reuse weights here, 
            # since the interpolation weights are the same. i.e.,:
            # w11*A[id1,id2,x1,y1] + w12*A[id1,id2,x1,y2] + 
            # w21*A[id1,id2,x2,y1] + w22*A[id1,id2,x2,y2] 
        end
    end

    return value
end


@btime Aitp = createitp(A,xgrd,ygrd,Nd1,Nd2); # Here I would like to use less allocations
@btime evalute(Aitp,Nd1,Nd2,2.0,2.0) # Here I would like to re-use weights to save time within the interpolation

Back in my previous Fortran days I wrote my own function to create interpolation weights but I was hoping to avoid doing that. There is no “real” reason for me to create the interpolation objects in `createitp’, since in the end I just want to interpolate the last two dimensions of A at different nodes of the first two dimension, except that it seems necessary to use Interpolations.

Here are a few related links: