# Pre allocating interpolant with Interpolations.jl

**URL:** <https://discourse.julialang.org/t/pre-allocating-interpolant-with-interpolations-jl/112712>\
**Category:** Numerics\
**Tags:** package, interpolations\
**Created:** [April 9, 2024, 8:25am UTC](https://discourse.julialang.org/t/pre-allocating-interpolant-with-interpolations-jl/112712 "2024-04-09T08:25:03Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Iddingsite](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/iddingsite/32/31840_2.png) [@Iddingsite](https://discourse.julialang.org/u/Iddingsite)\
**Post date:** [April 9, 2024, 8:25am UTC](https://discourse.julialang.org/t/pre-allocating-interpolant-with-interpolations-jl/112712/1 "2024-04-09T08:25:03Z")

</div>

Hi,

I have a loop where I need to create new interpolants each time. I would like to know if I can somehow pre allocate the interpolant in Interpolations.jl so I don’t have to create new memory allocation at each iteration of the loop. The loop is iterating through the same grid each time with only the values to interpolate that changes.

Here is a MWE taking the example in Interpolations.jl:

```julia
xs = 1:0.2:5
A = log.(xs)

# Create linear interpolation object without extrapolation
for _ in 1:5
    interp_linear = linear_interpolation(xs, A)
end

```

I tried to use the type of inter\_linear to pre allocate but it doesn’t seem to work:

```julia
interp_linear = Interpolations.Extrapolation{Float64, 1, ScaledInterpolation{Float64, 1, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, BSpline{Linear{Throw{OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}, BSpline{Linear{Throw{OnGrid}}}, Throw{Nothing}}

for _ in 1:5
    interp_linear .= linear_interpolation(xs, A)
end

```

error with `ERROR: CanonicalIndexError: setindex! not defined for Interpolations.Extrapolation{Float64, 1, ScaledInterpolation{Float64, 1, Interpolations.BSplineInterpolation{Float64, 1, Vector{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, BSpline{Linear{Throw{OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}, BSpline{Linear{Throw{OnGrid}}}, Throw{Nothing}}`

---

<div class="post-metadata">

**Author:** ![mkitti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkitti/32/12459_2.png) [@mkitti](https://discourse.julialang.org/u/mkitti)\
**Post date:** [April 11, 2024, 10:11am UTC](https://discourse.julialang.org/t/pre-allocating-interpolant-with-interpolations-jl/112712/2 "2024-04-11T10:11:13Z")

</div>

Sorry, this is not quite how Julia works in general. Specifically, storing a type in a variable is not a way to preallocatw memory. Also, trying to do broadcast assignment into the type will not perform an in-place operation.

What night help are functions that end with an excalamation mark, such as [`Interpolations.interpolate!`](https://juliamath.github.io/Interpolations.jl/stable/api/#Interpolations.interpolate!-Union%7BTuple%7BIT%7D,%20Tuple%7BTWeights%7D,%20Tuple%7BType%7BTWeights%7D,%20AbstractArray,%20IT%7D%7D%20where%20%7BTWeights,%20IT%3C:Union%7BNoInterp,%20Tuple%7BVararg%7BUnion%7BNoInterp,%20BSpline%7D,%20N%7D%20where%20N%7D,%20BSpline%7D%7D) however that does not do what you would like.

There are many interpolation packages. BSplineKit.jl’s [`interpolate!`](https://jipolanco.github.io/BSplineKit.jl/stable/interpolation/#BSplineKit.SplineInterpolations.interpolate!) looks like it my do what you would like.

---

<div class="post-metadata">

**Author:** ![Iddingsite](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/iddingsite/32/31840_2.png) [@Iddingsite](https://discourse.julialang.org/u/Iddingsite)\
**Post date:** [April 11, 2024, 1:29pm UTC](https://discourse.julialang.org/t/pre-allocating-interpolant-with-interpolations-jl/112712/3 "2024-04-11T13:29:09Z")

</div>

Thx for the answer! I will have a look on other packages then.

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [April 11, 2024, 1:49pm UTC](https://discourse.julialang.org/t/pre-allocating-interpolant-with-interpolations-jl/112712/4 "2024-04-11T13:49:14Z")

</div>

If you are doing linear interpolation, you can precompute the interpolation weights (even with Interpolations.jl although you need to dig into the internals a bit).

Beyond linear, I think you need the observed `y` values, so doesn’t the interpolation change every iteration anyway? In this case the best you can do is to try to reuse memory via in-place operations.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [April 11, 2024, 4:48pm UTC](https://discourse.julialang.org/t/pre-allocating-interpolant-with-interpolations-jl/112712/5 "2024-04-11T16:48:06Z")

</div>

Perhaps you are trying to solve a problem that we already solved in GeoStats.jl? Our interpolation routines are efficient: they preallocate a buffer for the results, and fill the buffer during traversal of large geospatial domains.

---

<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:** [April 11, 2024, 6:49pm UTC](https://discourse.julialang.org/t/pre-allocating-interpolant-with-interpolations-jl/112712/6 "2024-04-11T18:49:05Z")

</div>

While it may not be a good idea, you _can_ change the array your interpolating over if you’re using gridded interpolation:

```Julia
xs = collect(1:0.2:5)
ny = 10
A = [log(x)+y for x in xs, y in 1:ny]

itp = linear_interpolation(xs, A[:,1]) # Lets call this the placeholder interpolation object

# Replace the `coefs` field of the interpolation object with the new data
for _ in 1:ny
    for iy = 1:ny
        @views itp.itp.coefs[:] = A[:,iy]
        @assert itp(1.) == iy "$(itp(1.0)) == $(iy)"
    end
end

```

I do something like this in my code, it seemed to help on performance as I was re-creating some tens of thousands of interpolation objects. But given @mkitti reply above, I’m sure he’s right and if you’re doing this intelligently from the beginning (unlike me 🙂) it won’t help.

Regarding weights, I was never able to make it work, though this issue has some details: [Reuse weights with extrapolation · Issue #484 · JuliaMath/Interpolations.jl · GitHub](https://github.com/JuliaMath/Interpolations.jl/issues/484)