Weights and indices for linear interpolation/extrapolation

General Usage Performance
,
13

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 :roll_eyes:).
  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 :thinking:

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