# Weights and indices for linear interpolation/extrapolation

**URL:** <https://discourse.julialang.org/t/weights-and-indices-for-linear-interpolation-extrapolation/77580>\
**Category:** Performance\
**Tags:** question, interpolations\
**Created:** [March 8, 2022, 12:10pm UTC](https://discourse.julialang.org/t/weights-and-indices-for-linear-interpolation-extrapolation/77580 "2022-03-08T12:10:05Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![fredrikpaues](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fredrikpaues/32/34080_2.png) [@fredrikpaues](https://discourse.julialang.org/u/fredrikpaues)\
**Post date:** [March 8, 2022, 12:10pm UTC](https://discourse.julialang.org/t/weights-and-indices-for-linear-interpolation-extrapolation/77580/1 "2022-03-08T12:10:05Z")

</div>

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, but I haven’t been able to figure out how. The `gridpoints` vectors will typically be irregularly spaced.

```julia
function locate(
    value::AbstractFloat,
    gridpoints::AbstractVector{<:AbstractFloat};
    locate_below::Bool=false,
    locate_above::Bool=false,
)::Tuple{Integer, AbstractVector{<:AbstractFloat}}
    if value <= gridpoints[1]
        index = 1

        if !locate_below
            weights = [1.0, 0.0]
            return index, weights
        end
    elseif value >= gridpoints[end]
        index = length(gridpoints) - 1

        if !locate_above
            weights = [0.0, 1.0]
            return index, weights
        end
    else
        index = findlast(gridpoints .<= value)
    end

    weights = (
        [
            gridpoints[index + 1] - value,
            value - gridpoints[index]
        ] ./ (
            gridpoints[index + 1] - gridpoints[index]
        )
    )

    return index, weights
end
function locate(
    values::AbstractArray{<:AbstractFloat},
    gridpoints::AbstractVector{<:AbstractFloat};
    locate_below::Bool=false,
    locate_above::Bool=false,
)::Tuple{AbstractArray{Integer}, AbstractArray{<:AbstractFloat}}
    indices = ones(Integer, size(values))
    weights = zeros(2, size(values)...)

    for ix in CartesianIndices(values)
        indices[ix], weights[:, ix] = locate(
            values[ix],
            gridpoints;
            locate_below=locate_below,
            locate_above=locate_above
        )
    end

    return indices, weights
end

```

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 8, 2022, 1:31pm UTC](https://discourse.julialang.org/t/weights-and-indices-for-linear-interpolation-extrapolation/77580/2 "2022-03-08T13:31:44Z")

</div>

> [@fredrikpaues](#):
>
> 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.

Have you tried one of the existing interpolations packages, like [LinearInterpolations.jl](https://github.com/jw3126/LinearInterpolations.jl) or the other interpolations packages referenced on that page?

Some general tips:

1. Using things like `[1.0, 0.0]` as a return value, or anywhere in your code, allocates a small array on the heap. This is not a good idea in performance-critical code. e.g. use a tuple instead for a small fixed-length container (or use StaticArrays if you need other array operations). Doing `./` on arrays of length 2 is far slower than just doing two scalar divisions (unless you use StaticArrays) because of the heap allocations among other things, but it is fine on a tuple.
2. Your type annotations accomplish nothing for performance. (I would especially omit the return-type annotation.) `value` is probably overly restrictively typed; why not accept any `Real` value? See [Functions · The Julia Language](https://docs.julialang.org/en/v1/manual/functions/#Argument-type-declarations).
3. `findlast(gridpoints .<= value)` allocates large new array to hold `gridpoints .<= value` — not a good idea to repeat many times in an inner loop. You could do `findlast(<=(value), gridpoints)` to avoid this.
4. `findlast` is a linear search. Far better to require that the grid points be sorted and use `searchsortedfirst` to have a binary search. (If you have a regular grid, of course, you can avoid searching at all and just compute `(value - start) / Δx`.)

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [March 8, 2022, 4:10pm UTC](https://discourse.julialang.org/t/weights-and-indices-for-linear-interpolation-extrapolation/77580/3 "2022-03-08T16:10:13Z")

</div>

> [@fredrikpaues](#):
>
> ```julia
> function locate(
> values::AbstractArray{<:AbstractFloat},
> gridpoints::AbstractVector{<:AbstractFloat};
> locate_below::Bool=false,
> locate_above::Bool=false,
> )::Tuple{AbstractArray{Integer}, AbstractArray{<:AbstractFloat}}
> 
> ```

It is a really unfortunate habit you’ve picked up with the output annotations. You should absolutely stop doing it. In this case, it doesn’t just ‘do nothing’, it can actually significantly harm performance, since you are forcing your function to return an array with an abstract element type, `AbstractArray{Integer}`.

Here’s an example:

```julia
foo(x)::AbstractArray{Integer} = reverse(x)
bar(x) = reverse(x) # same but no annotation

1.8.0-beta1> x = rand(0:9, 1000);

1.8.0-beta1> @btime foo($x);
  5.133 μs (2 allocations: 15.88 KiB)

1.8.0-beta1> @btime bar($x);
  1.120 μs (1 allocation: 7.94 KiB)

```

If you then use this vector later in your calculations, it is _very bad_ for performance.

You should use output annotations only when there is a very good reason. Right now, they make the code less readable, confusing, restrictive, error prone, and slow.

(In fact, I suggest you start writing your code with zero type annotations, and then just add them when necessary. The exception being if you define new structs.)

Good luck! 😉

---

<div class="post-metadata">

**Author:** ![markmbaum](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/markmbaum/32/32745_2.png) [@markmbaum](https://discourse.julialang.org/u/markmbaum)\
**Post date:** [March 8, 2022, 4:33pm UTC](https://discourse.julialang.org/t/weights-and-indices-for-linear-interpolation-extrapolation/77580/4 "2022-03-08T16:33:38Z")

</div>

There is a simple bisection search function inside of [BasicInterpolators.jl](https://github.com/markmbaum/BasicInterpolators.jl) called `findcell` if that suits your particular need. You could copy it from [here](https://github.com/markmbaum/BasicInterpolators.jl/blob/6986d0418913b39207ccc7825a7f00822e434ef4/src/base.jl#L30) or just install the package and put `using BasicInterpolators: findcell` in your code.

The function has a [docstring](https://markmbaum.github.io/BasicInterpolators.jl/dev/misc/#BasicInterpolators.findcell), but it simply finds the index of the grid cell containing some value along a sorted vector of coordinates. Once you have that index, you can do a linear interpolation pretty easily.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [March 8, 2022, 4:44pm UTC](https://discourse.julialang.org/t/weights-and-indices-for-linear-interpolation-extrapolation/77580/5 "2022-03-08T16:44:56Z")

</div>

Just for fun, here’s an even more catastrophic example of the effects of that annotation:

```julia
1.8.0-beta1> x = 1:1_000_000;

1.8.0-beta1> @btime bar($x);
  7.900 ns (0 allocations: 0 bytes)

1.8.0-beta1> @btime foo($x);
  7.857 ms (999491 allocations: 22.88 MiB)

```

Note nanosecond vs millisecond. The type annotated version is literally a _million_ times slower, because `bar` does basically no work at all (independent of vector length), while `foo` scales linearly and allocates heavily.

---

<div class="post-metadata">

**Author:** ![fredrikpaues](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fredrikpaues/32/34080_2.png) [@fredrikpaues](https://discourse.julialang.org/u/fredrikpaues)\
**Post date:** [March 9, 2022, 5:36am UTC](https://discourse.julialang.org/t/weights-and-indices-for-linear-interpolation-extrapolation/77580/6 "2022-03-09T05:36:01Z")

</div>

Thank you!

I haven’t look into other interpolation packages except Interpolations.jl. Thanks for pointing out that there are alternatives.

1. In order to address this, it seems that I need to change how I use the function (as `i, w .= locate(value, gridpoints)` throws an error). I will be sure to look into it.

2. I mostly put in the type notation to discipline myself while coding and thought that they at worst didn’t affect performance at all. But after this and the replies by @DNF, I will be sure to stop with this unfortunate habit.

3. This improved performance of that line a lot!

4. Using the `searchsortedlast` function improved performance even more!

---

<div class="post-metadata">

**Author:** ![fredrikpaues](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fredrikpaues/32/34080_2.png) [@fredrikpaues](https://discourse.julialang.org/u/fredrikpaues)\
**Post date:** [March 9, 2022, 5:41am UTC](https://discourse.julialang.org/t/weights-and-indices-for-linear-interpolation-extrapolation/77580/7 "2022-03-09T05:41:42Z")

</div>

Thank you for pointing this out! As I wrote in my last reply I hadn’t considered that the type annotations could actually **hurt** performance. I will immediately stop with this.

1. But I remember having read that type annotations sometimes can help the compiler. How do I figure out if a particular variable would benefit from a type annotation?

2. If I hypothetically switched to concrete types, could that improve performance? Or would it at best not hurt it?

3. In the code in my original post, I dabbled in using multiple dispatch. To do that, I’ve gotten the impression that I **need** to use type annotations. Should I just use `Number` for the first/scalar case and `AbstractArray` for the second? (Could I somehow write my function so that I could call it with dot notation, i.e. `locate.(values, gridpoints)`?)

---

<div class="post-metadata">

**Author:** ![fredrikpaues](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fredrikpaues/32/34080_2.png) [@fredrikpaues](https://discourse.julialang.org/u/fredrikpaues)\
**Post date:** [March 9, 2022, 5:47am UTC](https://discourse.julialang.org/t/weights-and-indices-for-linear-interpolation-extrapolation/77580/8 "2022-03-09T05:47:31Z")

</div>

I will try this, but I’m a bit skeptical. As the `findcell` function only returns the index, it seems that I would have to check the boundaries another time to compute the weights. Perhaps I could improve performance by computing the index and the weights myself using the bisection method as `findcell` 🤔

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [March 9, 2022, 8:25am UTC](https://discourse.julialang.org/t/weights-and-indices-for-linear-interpolation-extrapolation/77580/9 "2022-03-09T08:25:26Z")

</div>

> [@fredrikpaues](#):
>
> But I remember having read that type annotations sometimes can help the compiler. How do I figure out if a particular variable would benefit from a type annotation?

I think that as a first approximation you can assume ‘never’.There may be some few cases, it used to be that things like

```julia
x += boolean_expression::Int

```

where `x` is an `Int`, but I don’t think that’s been needed for a long time. There’s also a long-standing issue with captured variables in closures, that I think is still ongoing.

But if you yourself can predict the output type based solely on the input types, then the compiler can too, virtually always, and better than a human. But if you have an actual type unstable function, and you have privileged knowledge about the computation from outside the type domain, then you can try. But first you should verify that there is an _actual_ instability, using e.g. `@code_warntype`, otherwise, there is _definitely_ no point.

> [@fredrikpaues](#):
>
> If I hypothetically switched to concrete types, could that improve performance? Or would it at best not hurt it?

In function signatures, no. In output signatures, concrete is better than abstract, apparently, but it’s better to leave it off, unless you discover a _bona fide_ instability.

In struct definitions, yes, you should use concrete or parametric types.

> [@fredrikpaues](#):
>
> In the code in my original post, I dabbled in using multiple dispatch. To do that, I’ve gotten the impression that I **need** to use type annotations.

Yes, you need type annotations to define different methods of the same function for different input types. This is actually the primary purpose of function type signatures.

> [@fredrikpaues](#):
>
> Should I just use `Number` for the first/scalar case and `AbstractArray` for the second? (Could I somehow write my function so that I could call it with dot notation, i.e. `locate.(values, gridpoints)` ?)

Yes, from a very brief look, I think you need only one method. Just keep the first one, you could give it the signature

```julia
locate(::Real, ::AbstractVector, ::Bool, ::Bool)

```

or just drop the types completely. And then, when you have a vector of values, you call

```julia
locate.(values, Ref(gridpoints), locate_below, locate_above)

```

The `Ref()` protects `gridpoints` against getting broadcasted over. This will give you a somewhat different output type, something like `Vector{Tuple{Int, Vector}}` (don’t annotate your function with it, though 😉 ) I think someone already mentioned that `weights` ought to be a tuple instead of a vector, if you change that, the output would be, approximately, `Tuple{Int, Tuple{Float64, Float64}}`.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 9, 2022, 12:33pm UTC](https://discourse.julialang.org/t/weights-and-indices-for-linear-interpolation-extrapolation/77580/10 "2022-03-09T12:33:36Z")

</div>

> [@fredrikpaues](#):
>
> In order to address this, it seems that I need to change how I use the function (as `i, w .= locate(value, gridpoints)` throws an error). I will be sure to look into it.

If you use a tuple for `w` instead of an array, you can just use `i, w = locate(...)`.

---

<div class="post-metadata">

**Author:** ![markmbaum](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/markmbaum/32/32745_2.png) [@markmbaum](https://discourse.julialang.org/u/markmbaum)\
**Post date:** [March 11, 2022, 10:45pm UTC](https://discourse.julialang.org/t/weights-and-indices-for-linear-interpolation-extrapolation/77580/11 "2022-03-11T22:45:51Z")

</div>

Here is an example script that wraps the `findcell` function to handle boundaries with extrapolation or not, plotting a test case with with the `PyPlot` package.

```julia
using PyPlot

pygui(true)

##

function findcell(xᵢ, X)
    #number of grid cells
    n = length(X)
    #handle boundaries
    @inbounds (xᵢ <= X[1]) && return(1)
    @inbounds (xᵢ >= X[n]) && return(n-1)
    #bisection search for the containing cell
    L = 0
    H = n
    while H - L > 1
        M = (H + L) ÷ 2
        @inbounds if X[M] > xᵢ
            H = M
        else
            L = M
        end
    end
    return L
end

function locate(value, gridpoints, locate_below=false, locate_above=true)
    #get the cell index or boundaries
    index = findcell(value, gridpoints)
    #width of target cell
    Δ = gridpoints[index+1] - gridpoints[index]
    #get interpolation weights, handling boundaries
    weights = if (value < gridpoints[1]) & !locate_below
        (1.0, 0.0)
    elseif (value > gridpoints[end]) & !locate_above
        (0.0, 1.0)
    else
        (
            (gridpoints[index+1] - value)/Δ,
            (value - gridpoints[index])/Δ
        )
    end
    return index, weights
end

##

#grid coordinates
gridpoints = [0.5, 1, 2, 2.5, 4]
#values to interpolate
gridvalues = 2*rand(5) .- 1

#coordinates to interpolate
x = LinRange(-1, 6, 1000)
#interpolation values
y = zeros(length(x))

figure()
#plot the interpolation points
plot(gridpoints, gridvalues, ".")
#interpolate and plot with different boundary behaviors
for extrapolate ∈ (true, false)
    for i ∈ 1:length(x)
        index, weights = locate(x[i], gridpoints, extrapolate, extrapolate)
        y[i] = sum(weights .* gridvalues[index:index+1])
    end
    plot(x, y, label="extrapolate=$extrapolate")
end
legend()

```

Hopefully this is helpful. There shouldn’t be any need to annotate any of the function inputs or outputs here.

As you did in the initial post, you can write another function to wrap the `locate` function if you need to store a bunch of indices and weights.

---

<div class="post-metadata">

**Author:** ![fredrikpaues](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fredrikpaues/32/34080_2.png) [@fredrikpaues](https://discourse.julialang.org/u/fredrikpaues)\
**Post date:** [March 12, 2022, 10:07am UTC](https://discourse.julialang.org/t/weights-and-indices-for-linear-interpolation-extrapolation/77580/12 "2022-03-12T10:07:32Z")

</div>

Thanks! I ended up doing something very similar, but that avoids the duplicated boundary handling (once in `findcell` and once in `locate`).

---

<div class="post-metadata">

**Author:** ![fredrikpaues](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fredrikpaues/32/34080_2.png) [@fredrikpaues](https://discourse.julialang.org/u/fredrikpaues)\
**Post date:** [March 12, 2022, 10:21am UTC](https://discourse.julialang.org/t/weights-and-indices-for-linear-interpolation-extrapolation/77580/13 "2022-03-12T10:21:25Z")

</div>

This is what I have right now.

1. Weights are tuples (as advised by @stevengj).
2. I have tried to avoid the type annotations (@stevengj and @DNF) except for multiple dispatch and the booleans. (Perhaps I should skip the `Bool`:s as well 🙄).
3. I use binary search for the interior points (@stevengj and @markmbaum).
4. I found that dot notation (`locate.()`) performed worse than just implementing the `values::AbstractArray{<:Real}` case myself. That way, I could also (a) control the structure of the output to something that was a little more tractable and (b) ensure that `gridpoints` was a tuple, which turned out to improve performance.

Perhaps performance could be improved by using `@inbounds` when checking the boundaries (as is done in `findcell` suggested by @markmbaum), but I haven’t managed to figure out what exactly that does 🤔

```julia
function locate(
    value::Real,
    gridpoints;
    locate_below::Bool=false,
    locate_above::Bool=false,
)
    if value <= gridpoints[1]
        index = 1

        if !locate_below
            weights = (1.0, 0.0)
            return index, weights
        end
    elseif value >= gridpoints[end]
        index = length(gridpoints) - 1

        if !locate_above
            weights = (0.0, 1.0)
            return index, weights
        end
    else
        # Implement binary search as
        # `index = searchsortedlast(gridpoints, value)`
        # but with support for tuple gridpoints
        low = 0
        high = length(gridpoints) + 1
        @inbounds while low < high - 1
            mid = low + ((high - low) >>> 0x01) # floor(Integer, (low + high) / 2)
            if value < gridpoints[mid]
                high = mid
            else
                low = mid
            end
        end

        index = low
    end

    weights = (
        (
            gridpoints[index + 1] - value,
            value - gridpoints[index]
        ) ./ (
            gridpoints[index + 1] - gridpoints[index]
        )
    )

    return index, weights
end
function locate(
    values::AbstractArray{<:Real},
    gridpoints;
    locate_below::Bool=false,
    locate_above::Bool=false,
)
    if !(typeof(gridpoints) <: Tuple{Vararg{<:Real}})
        gridpoints = Tuple{Vararg{<:Real}}(gridpoints)
    end

    indices = ones(Integer, size(values))
    weights = Array{Tuple, ndims(values)}(undef, size(values))

    for ix in eachindex(values)
        indices[ix], weights[ix] = locate(
            values[ix],
            gridpoints;
            locate_below=locate_below,
            locate_above=locate_above
        )
    end

    return indices, weights
end

```

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 12, 2022, 1:46pm UTC](https://discourse.julialang.org/t/weights-and-indices-for-linear-interpolation-extrapolation/77580/14 "2022-03-12T13:46:11Z")

</div>

> [@fredrikpaues](#):
>
> ```julia
> Implement binary search as
> # `index = searchsortedlast(gridpoints, value)`
> # but with support for tuple gridpoints
> 
> ```

I’m confused, why can’t you use `searchsortedlast`? You almost certainly don’t want to use a tuple to store a long array of grid points.

---

<div class="post-metadata">

**Author:** ![fredrikpaues](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fredrikpaues/32/34080_2.png) [@fredrikpaues](https://discourse.julialang.org/u/fredrikpaues)\
**Post date:** [March 12, 2022, 2:16pm UTC](https://discourse.julialang.org/t/weights-and-indices-for-linear-interpolation-extrapolation/77580/15 "2022-03-12T14:16:40Z")

</div>

I found it was faster when I tried it with `@btime`, which is why I change `gridpoints` to a tuple in the array case. Maybe I was unsuccessful in my thinking here 🤔 I will check again in a few hours and post my tests here.

---

<div class="post-metadata">

**Author:** ![fredrikpaues](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fredrikpaues/32/34080_2.png) [@fredrikpaues](https://discourse.julialang.org/u/fredrikpaues)\
**Post date:** [March 12, 2022, 7:29pm UTC](https://discourse.julialang.org/t/weights-and-indices-for-linear-interpolation-extrapolation/77580/16 "2022-03-12T19:29:10Z")

</div>

So these are the tests that I ran.

First, I ran these but without the part in `locate` that converts `gridpoints` to a tuple.

```julia-repl
julia> logrange(x1, x2, n) = (10^y for y in range(log10(x1), log10(x2); length=n));
julia> g_vector = collect(logrange(0.1, 0.9, 100));
julia> g_tuple = Tuple{Real}(g_vector);
julia> using StaticArrays;
julia> g_svector = SVector{length(g_vector), Real}(g_vector);
julia> values = rand(200);
julia> using BenchmarkTools;
julia> @btime locate(values, g_vector); @btime locate(values, g_tuple); @btime locate(values, g_svector);
3.612 μs (203 allocations: 9.81 KiB)
2.100 μs (203 allocations: 9.81 KiB)
187.400 μs (5183 allocations: 103.09 KiB)
julia> @btime locate(values, g_vector); @btime locate(values, g_tuple); @btime locate(values, g_svector);
3.600 μs (203 allocations: 9.81 KiB)
2.089 μs (203 allocations: 9.81 KiB)
188.200 μs (5183 allocations: 103.09 KiB)

```

And then I ran the same tests but **did** convert `gridpoints` to a tuple.

```julia-repl
julia> @btime locate(values, g_vector); @btime locate(values, g_tuple); @btime locate(values, g_svector);
32.100 μs (1115 allocations: 32.20 KiB)
20.200 μs (1009 allocations: 28.78 KiB)
46.300 μs (1519 allocations: 42.36 KiB)
julia> @btime locate(values, g_vector); @btime locate(values, g_tuple); @btime locate(values, g_svector);
32.400 μs (1115 allocations: 32.20 KiB)
20.200 μs (1009 allocations: 28.78 KiB)
46.800 μs (1519 allocations: 42.36 KiB)

```

So, tuples perform the best. But apparently my decision to convert was rather stupid. (I think that I originally tested with fewer nodes in my grids, but these number are more in line with what I will use the function for.)