# Faster/performant interpolation of very big arrays

**URL:** <https://discourse.julialang.org/t/faster-performant-interpolation-of-very-big-arrays/99996>\
**Category:** Performance\
**Created:** [June 7, 2023, 11:49am UTC](https://discourse.julialang.org/t/faster-performant-interpolation-of-very-big-arrays/99996 "2023-06-07T11:49:34Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![eirikeb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eirikeb/32/11674_2.png) [@eirikeb](https://discourse.julialang.org/u/eirikeb)\
**Post date:** [June 7, 2023, 11:49am UTC](https://discourse.julialang.org/t/faster-performant-interpolation-of-very-big-arrays/99996/1 "2023-06-07T11:49:34Z")

</div>

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:

```julia
## 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:

> <https://github.com/JuliaMath/Interpolations.jl/issues/484>
>
> I'm finding the \[documentation on how to reuse weights\](http://juliamath.github.…io/Interpolations.jl/latest/devdocs/#Interpolant-usage) quite difficult to comprehend. But I would like to compute weights for an extrapolant/interpolant. My best attempt at replicating the code in the documentation but with my own interpolation object (below) results in a \`MethodError\`.
> 
> \`\`\`julia-repl
> julia\> using Interpolations
> julia\> x = collect(range(0.1, 0.9, 100));
> julia\> v = 2x;
> julia\> itp = LinearInterpolation(x, v; extrapolation\_bc=((Flat(), Line()),));
> julia\> itp(0.5) # Works just fine :)
> 1.0
> julia\> wis = Interpolations.weightedindexes((Interpolations.value\_weights,), Interpolations.itpinfo(itp)..., 0.5);
> ERROR: MethodError: no method matching weightedindexes(::Tuple{typeof(Interpolations.value\_weights)}, ::Tuple{Gridded{Linear{Throw{OnGrid}}}}, ::Tuple{Base.OneTo{Int64}}, ::Float64)
> Closest candidates are:
> weightedindexes(::F, ::Tuple{Vararg{Interpolations.Flag, N}}, ::Tuple{Vararg{AbstractVector, N}}, ::Tuple{Vararg{Number, N}}) where {F, N} at C:\\...\\.julia\\packages\\Interpolations\\Glp9h\\src\\b-splines\\indexing.jl:64
> \`\`\`

> [@Weights and indices for linear interpolation/extrapolation](https://discourse.julialang.org/t/weights-and-indices-for-linear-interpolation-extrapolation/77580):
>
> A performance-critical function of a bigger program that I’m writing consists of locating points on a grid. It amounts to finding the indices and weights needed for a linear interpolation/extrapolation—without doing the actual interpolation/extrapolation. What I have written as of yet is below. I’ve been coding Julia for about a week and would very much appreciate some pointers on how I might improve what I’ve written. Perhaps this functionality can be found in the Interpolations.jl package, b…