Public API

ChainRulesCore.rruleMethod
ChainRulesCore.rrule(itp::AbstractInterpolation, x...)

ChainRulesCore.jl rrule for integration with automatic differentiation libraries. Note that it gives the gradient only with respect to the evaluation point x, and not the data inside itp.

source
Interpolations.ConstantInterpolationFunction
etp = ConstantInterpolation(knots, A; extrapolation_bc=Throw())

A shorthand for extrapolate(interpolate(knots, A, scheme), extrapolation_bc), where scheme is either BSpline(Constant()) or Gridded(Constant()) depending on whether knots are ranges or vectors.

source
Interpolations.CubicSplineInterpolationFunction
etp = CubicSplineInterpolation(knots, A; bc=Line(OnGrid()), extrapolation_bc=Throw())

A shorthand for extrapolate(interpolate(knots, A, BSpline(Cubic(bc))), extrapolation_bc).

source
Interpolations.LinearInterpolationFunction
etp = LinearInterpolation(knots, A; extrapolation_bc=Throw())

A shorthand for extrapolate(interpolate(knots, A, scheme), extrapolation_bc), where scheme is either BSpline(Linear()) or Gridded(Linear()) depending on whether knots are ranges or vectors.

source
Interpolations.boundsMethod
bounds(itp::AbstractInterpolation)

Return the bounds of the domain of itp as a tuple of (min, max) pairs for each coordinate. This is best explained by example:

julia> itp = interpolate([1 2 3; 4 5 6], BSpline(Linear()));

julia> bounds(itp)
((1, 2), (1, 3))

julia> data = 1:3;

julia> knots = ([10, 11, 13.5],);

julia> itp = interpolate(knots, data, Gridded(Linear()));

julia> bounds(itp)
((10.0, 13.5),)
source
Interpolations.extrapolateMethod

extrapolate(itp, scheme) adds extrapolation behavior to an interpolation object, according to the provided scheme.

The scheme can take any of these values:

  • Throw - throws a BoundsError for out-of-bounds indices
  • Flat - for constant extrapolation, taking the closest in-bounds value
  • Line - linear extrapolation (the wrapped interpolation object must support gradient)
  • Reflect - reflecting extrapolation (indices must support mod)
  • Periodic - periodic extrapolation (indices must support mod)

You can also combine schemes in tuples. For example, the scheme (Line(), Flat()) will use linear extrapolation in the first dimension, and constant in the second.

Finally, you can specify different extrapolation behavior in different direction. ((Line(),Flat()), Flat()) will extrapolate linearly in the first dimension if the index is too small, but use constant etrapolation if it is too large, and always use constant extrapolation in the second dimension.

source
Interpolations.extrapolateMethod

extrapolate(itp, fillvalue) creates an extrapolation object that returns the fillvalue any time the indexes in itp(x1,x2,...) are out-of-bounds.

source
Interpolations.interpolateMethod
itp = interpolate(A, interpmode, gridstyle, λ, k)

Interpolate an array A in the mode determined by interpmode and gridstyle with regularization following [1], of order k and constant λ. interpmode may be one of

  • BSpline(NoInterp())
  • BSpline(Linear())
  • BSpline(Quadratic(BC())) (see BoundaryCondition)
  • BSpline(Cubic(BC()))

It may also be a tuple of such values, if you want to use different interpolation schemes along each axis.

gridstyle should be one of OnGrid() or OnCell().

k corresponds to the derivative to penalize. In the limit λ->∞, the interpolation function is a polynomial of order k-1. A value of 2 is the most common.

λ is non-negative. If its value is zero, it falls back to non-regularized interpolation.

[1] https://projecteuclid.org/euclid.ss/1038425655.

source
Interpolations.interpolateMethod
itp = interpolate(A, interpmode)

Interpolate an array A in the mode determined by interpmode. interpmode may be one of

  • NoInterp()
  • BSpline(Constant())
  • BSpline(Linear())
  • BSpline(Quadratic(bc)) (see BoundaryCondition)
  • BSpline(Cubic(bc))

It may also be a tuple of such values, if you want to use different interpolation schemes along each axis.

source
Interpolations.interpolateMethod
itp = interpolate((nodes1, nodes2, ...), A, interpmode)

Interpolate an array A on a non-uniform but rectangular grid specified by the given nodes, in the mode determined by interpmode.

interpmode may be one of

  • NoInterp()
  • Gridded(Constant())
  • Gridded(Linear())

It may also be a tuple of such values, if you want to use different interpolation schemes along each axis.

source
Interpolations.knotsMethod
knots(itp::AbstractInterpolation)
knots(etp::AbstractExtrapolation)

Returns an iterator over knot locations for an AbstractInterpolation or AbstractExtrapolation.

Iterator will yield scalar values for interpolations over a single dimension, and tuples of coordinates for higher dimension interpolations. Iteration over higher dimensions is taken as the product of knots along each dimension.

ie. Iterator.product(knots on first dim, knots on 2nd dim,...)

Extrapolations with Periodic or Reflect boundary conditions, will produce an infinite sequence of knots.

Example

julia> using Interpolations;

julia> etp = LinearInterpolation([1.0, 1.2, 2.3, 3.0], rand(4); extrapolation_bc=Periodic());

julia> Iterators.take(knots(etp), 5) |> collect
5-element Vector{Float64}:
 1.0
 1.2
 2.3
 3.0
 3.2
source
Interpolations.knotsbetweenMethod
knotsbetween(iter; start, stop)
knotsbetween(iter, start, stop)

Iterate over all knots of iter such that start < k < stop.

iter can be an AbstractInterpolation, or the output of knots (ie. a KnotIterator or ProductIterator wrapping KnotIterator)

If start is not provided, iteration will start from the first knot. An ArgumentError will be raised if both start and stop are not provided.

If no such knots exists will return a KnotIterator with length 0

Example

julia> using Interpolations;

julia> etp = LinearInterpolation([1.0, 1.2, 2.3, 3.0], rand(4); extrapolation_bc=Periodic());

julia> knotsbetween(etp; start=38, stop=42) |> collect
6-element Vector{Float64}:
 38.3
 39.0
 39.2
 40.3
 41.0
 41.2
source
Interpolations.scaleMethod

scale(itp, xs, ys, ...) scales an existing interpolation object to allow for indexing using other coordinate axes than unit ranges, by wrapping the interpolation object and transforming the indices from the provided axes onto unit ranges upon indexing.

The parameters xs etc must be either ranges or linspaces, and there must be one coordinate range/linspace for each dimension of the interpolation object.

For every NoInterp dimension of the interpolation object, the range must be exactly 1:size(itp, d).

source
Interpolations.CardinalMonotonicInterpolationType
CardinalMonotonicInterpolation(tension)

Cubic cardinal splines, uses tension parameter which must be between [0,1] Cubin cardinal splines can overshoot for non-monotonic data (increasing tension reduces overshoot).

source
Interpolations.ConstantType

Constant b-splines are nearest-neighbor interpolations, and effectively return A[round(Int,x)] when interpolating.

source
Interpolations.CubicType
Cubic(bc::BoundaryCondition)

Indicate that the corresponding axis should use cubic interpolation.

Extended help

Assuming uniform knots with spacing 1, the ith piece of cubic spline implemented here is defined as follows.

y_i(x) = cm p(x-i) + c q(x-i) + cp q(1- (x-i)) + cpp p(1 - (x-i))

where

p(δx) = 1/6 * (1-δx)^3
q(δx) = 2/3 - δx^2 + 1/2 δx^3

and the values cX for X ∈ {m, _, p, pp} are the pre-filtered coefficients.

For future reference, this expands out to the following polynomial:

y_i(x) = 1/6 cm (1+i-x)^3 + c (2/3 - (x-i)^2 + 1/2 (x-i)^3) +
         cp (2/3 - (1+i-x)^2 + 1/2 (1+i-x)^3) + 1/6 cpp (x-i)^3

When we derive boundary conditions we will use derivatives y_0'(x) and y_0''(x)

source
Interpolations.CubicHermiteType
 CubicHermite

This type is purposely left undocumented since the interface is expected to radically change in order to make it conform to the AbstractInterpolation interface. Consider this API to be highly unstable and non-public, and that breaking changes to this code could be made in a point release.

source
Interpolations.FritschButlandMonotonicInterpolationType
FritschButlandMonotonicInterpolation

Monotonic interpolation based on Fritsch & Butland (1984), "A Method for Constructing Local Monotone Piecewise Cubic Interpolants", doi:10.1137/0905021.

Faster than FritschCarlsonMonotonicInterpolation (only requires one pass) but somewhat higher apparent "tension".

source
Interpolations.FritschCarlsonMonotonicInterpolationType
FritschCarlsonMonotonicInterpolation

Monotonic interpolation based on Fritsch & Carlson (1980), "Monotone Piecewise Cubic Interpolation", doi:10.1137/0717021.

Tangents are first initialized, then adjusted if they are not monotonic can overshoot for non-monotonic data

source
Interpolations.InPlaceType

InPlace(gt) is a boundary condition that allows prefiltering to occur in-place (it typically requires padding)

source
Interpolations.InPlaceQType

InPlaceQ(gt) is similar to InPlace(gt), but is exact when the values being interpolated arise from an underlying quadratic. It is primarily useful for testing purposes, allowing near-exact (to machine precision) comparisons against ground truth.

source
Interpolations.LanczosType
Lanczos{N}(a=4)

Lanczos resampling via a kernel with scale parameter a and support over N neighbors.

This form of interpolation is merely the discrete convolution of the samples with a Lanczos kernel of size a. The size is directly related to how "far" the interpolation will reach for information, and has O(N^2) impact on runtime. An alternative implementation matching lanczos4 from OpenCV is available as Lanczos4OpenCV.

source
Interpolations.Lanczos4OpenCVType
Lanczos4OpenCV()

Alternative implementation of Lanczos resampling using algorithm lanczos4 function of OpenCV: https://github.com/opencv/opencv/blob/de15636724967faf62c2d1bce26f4335e4b359e5/modules/imgproc/src/resize.cpp#L917-L946

source
Interpolations.LinearType
Linear()

Indicate that the corresponding axis should use linear interpolation.

Extended help

Assuming uniform knots with spacing 1, the ith piece of linear b-spline implemented here is defined as follows.

y_i(x) = c p(x) + cp p(1-x)

where

p(δx) = x

and the values cX for X ∈ {_, p} are the coefficients.

Linear b-splines are naturally interpolating, and require no prefiltering; there is therefore no need for boundary conditions to be provided.

Also, although the implementation is slightly different in order to re-use the framework built for general b-splines, the resulting interpolant is just a piecewise linear function connecting each pair of neighboring data points.

source
Interpolations.NoInterpType

NoInterp() indicates that the corresponding axis must use integer indexing (no interpolation is to be performed)

source
Interpolations.OnCellType

OnCell() indicates that the boundary condition applies a half-gridspacing beyond the first & last nodes

source
Interpolations.QuadraticType
Quadratic(bc::BoundaryCondition)

Indicate that the corresponding axis should use quadratic interpolation.

Extended help

Assuming uniform knots with spacing 1, the ith piece of quadratic spline implemented here is defined as follows:

y_i(x) = cm p(x-i) + c q(x) + cp p(1-(x-i))

where

p(δx) = (δx - 1)^2 / 2
q(δx) = 3/4 - δx^2

and the values for cX for X ∈ {m,_,p} are the pre-filtered coefficients.

For future reference, this expands to the following polynomial:

y_i(x) = cm * 1/2 * (x-i-1)^2 + c * (3/4 - x + i)^2 + cp * 1/2 * (x-i)^2

When we derive boundary conditions we will use derivatives y_1'(x-1) and y_1''(x-1)

source
Interpolations.SteffenMonotonicInterpolationType
SteffenMonotonicInterpolation

Monotonic interpolation based on Steffen (1990), "A Simple Method for Monotonic Interpolation in One Dimension", http://adsabs.harvard.edu/abs/1990A%26A...239..443S

Only one pass, results usually between FritschCarlson and FritschButland.

source

Internal API

Interpolations.boundstepFunction

Returns half the width of one step of the range.

This function is used to calculate the upper and lower bounds of OnCell interpolation objects.

source
Interpolations.count_interp_dimsMethod
n = count_interp_dims(ITP)

Count the number of dimensions along which type ITP is interpolating. NoInterp dimensions do not contribute to the sum.

source
Interpolations.deduplicate_knots!Method
Interpolations.deduplicate_knots!(knots; move_knots = false)

Makes knots unique by incrementing repeated but otherwise sorted knots using `nextfloat`.
If keyword `move_knots` is true, then `nextfloat` will be applied successively until knots
are unique. Otherwise, a warning will be issued.

# Example

```jldoctest
julia> knots = [-8.0, 0.0, 20.0, 20.0]
4-element Vector{Float64}:
-8.0
0.0
20.0
20.0

julia> Interpolations.deduplicate_knots!(knots)
4-element Vector{Float64}:
-8.0
0.0
20.0
20.000000000000004

julia> Interpolations.deduplicate_knots!([1.0, 1.0, 1.0, nextfloat(1.0), nextfloat(1.0)]; move_knots = true)
5-element Vector{Float64}:
1.0
1.0000000000000002
1.0000000000000004
1.0000000000000007
1.0000000000000009
```
source
Interpolations.gradient_weightsFunction
w = gradient_weights(degree, δx)

Compute the weights for interpolation of the gradient at an offset δx from the "base" position. degree describes the interpolation scheme.

Example

julia> Interpolations.gradient_weights(Linear(), 0.2)
(-1.0, 1.0)

This defines the gradient of a linear interpolation at 3.2 as y[4] - y[3].

source
Interpolations.hessian_weightsFunction
w = hessian_weights(degree, δx)

Compute the weights for interpolation of the hessian at an offset δx from the "base" position. degree describes the interpolation scheme.

Example

julia> Interpolations.hessian_weights(Linear(), 0.2)
(0.0, 0.0)

Linear interpolation uses straight line segments, so the second derivative is zero.

source
Interpolations.inner_system_diagsFunction
dl, d, du = inner_system_diags{T,IT}(::Type{T}, n::Int, ::Type{IT})

Helper function to generate the prefiltering equation system: generates the diagonals for a n-by-n tridiagonal matrix with eltype T corresponding to the interpolation type IT.

dl, d, and du are intended to be used e.g. as in M = Tridiagonal(dl, d, du)

source
Interpolations.nextknotidxMethod
nextknotidx(iter::KnotIterator, x)

Returns the index of the first knot such that x < k or nothing if no such knot exists.

New boundary conditions should define:

nextknotidx(::Type{<:NewBoundaryCondition}, knots::Vector, x)

Where knots is iter.knots and NewBoundaryCondition is the new boundary conditions. This method is expected to handle values of x that are both inbounds or extrapolated.

source
Interpolations.prefiltering_systemFunction
M, b = prefiltering_system{T,TC,GT<:GridType,D<:Degree}m(::T, ::Type{TC}, n::Int, ::Type{D}, ::Type{GT})

Given element types (T, TC) and interpolation scheme (GT, D) as well the number of rows in the data input (n), compute the system used to prefilter spline coefficients. Boundary conditions determine the values on the first and last rows.

Some of these boundary conditions require that these rows have off-tridiagonal elements (e.g the [1,3] element of the matrix). To maintain the efficiency of solving tridiagonal systems, the Woodbury matrix identity is used to add additional elements off the main 3 diagonals.

The filtered coefficients are given by solving the equation system

M * c = v + b

where c are the sought coefficients, and v are the data points.

source
Interpolations.prefiltering_systemMethod

Cubic{Flat}, OnCell amounts to setting y_1'(x) = 0 at x = 1/2. Applying this condition yields

-9/8 cm + 11/8 c - 3/8 cp + 1/8 cpp = 0

or, equivalently,

-9 cm + 11 c -3 cp + 1 cpp = 0

(Note that we use y_1'(x) although it is strictly not valid in this domain; if we were to use y_0'(x) we would have to introduce new coefficients, so that would not close the system. Instead, we extend the outermost polynomial for an extra half-cell.)

source
Interpolations.prefiltering_systemMethod

Cubic{Line} OnCell amounts to setting y_1''(x) = 0 at x = 1/2. Applying this condition yields

3 cm -7 c + 5 cp -1 cpp = 0

(Note that we use y_1'(x) although it is strictly not valid in this domain; if we were to use y_0'(x) we would have to introduce new coefficients, so that would not close the system. Instead, we extend the outermost polynomial for an extra half-cell.)

source
Interpolations.prefiltering_systemMethod

Cubic{Free} OnGrid and Cubic{Free} OnCell amount to requiring an extra continuous derivative at the second-to-last cell boundary; this means y_1'''(2) = y_2'''(2), yielding

1 cm -3 c + 3 cp -1 cpp = 0
source
Interpolations.prefiltering_systemMethod

Cubic{Periodic} OnGrid closes the system by looking at the coefficients themselves as periodic, yielding

c0 = c(N+1)

where N is the number of data points.

source
Interpolations.prefiltering_systemMethod

Quadratic{Free} OnGrid and Quadratic{Free} OnCell amount to requiring an extra continuous derivative at the second-to-last cell boundary; this means that y_1''(3/2) = y_2''(3/2), yielding

1 cm -3 c + 3 cp - cpp = 0
source
Interpolations.prefiltering_systemMethod

Quadratic{Line} OnGrid and Quadratic{Line} OnCell amount to setting y_1''(x) = 0 at x=1 and x=1/2 respectively. Since y_i''(x) is independent of x for a quadratic b-spline, these both yield

1 cm -2 c + 1 cp = 0
source
Interpolations.prefiltering_systemMethod

Quadratic{Periodic} OnGrid and Quadratic{Periodic} OnCell close the system by looking at the coefficients themselves as periodic, yielding

c0 = c(N+1)

where N is the number of data points.

source
Interpolations.priorknotidxMethod
priorknotidx(iter::KnotIterator, x)

Returns the index of the last knot such that k < x or nothing ig not such knot exists.

New boundary conditions should define

priorknotidx(::Type{<:NewBoundaryCondition}, knots::Vector, x)

Where knots is iter.knots and NewBoundaryCondition is the new boundary condition. This method is expected to handle values of x that are both inbounds or extrapolated.

source
Interpolations.rescale_gradientFunction

rescale_gradient(r::AbstractRange)

Implements the chain rule dy/dx = dy/du * du/dx for use when calculating gradients with scaled interpolation objects.

source
Interpolations.show_rangedMethod
show_ranged(io, X, knots)

A replacement for the default array-show for types that may not have the canonical evaluation points. rngs is the tuple of knots along each axis.

source
Interpolations.value_weightsFunction
w = value_weights(degree, δx)

Compute the weights for interpolation of the value at an offset δx from the "base" position. degree describes the interpolation scheme.

Example

julia> Interpolations.value_weights(Linear(), 0.2)
(0.8, 0.2)

This corresponds to the fact that linear interpolation at x + 0.2 is 0.8*y[x] + 0.2*y[x+1].

source
Interpolations.weightedindexesMethod
weightedindexes(fs, itpflags, nodes, xs)

Compute WeightedIndex values for evaluation at the position xs.... fs is a function or tuple of functions indicating the types of index required, typically value_weights, gradient_weights, and/or hessian_weights. itpflags and nodes can be obtained from itpinfo(itp)....

See the "developer documentation" for further information.

source
Interpolations.BSplineInterpolationType
BSplineInterpolation{T,N,TCoefs,IT,Axs}

An interpolant-type for b-spline interpolation on a uniform grid with integer nodes. T indicates the element type for operations like collect(itp), and may also agree with the values obtained from itp(x, y, ...) at least for certain types of x and y. N is the dimensionality of the interpolant. The remaining type-parameters describe the types of fields:

  • the coefs field holds the interpolation coefficients. Depending on prefiltering, these may or may not be the same as the supplied array of interpolant values.
  • parentaxes holds the axes of the parent. Depending on prefiltering this may be "narrower" than the axes of coefs.
  • it holds the interpolation type, e.g., BSpline(Linear()) or (BSpline(Quadratic(OnCell()),BSpline(Linear())).

BSplineInterpolation objects are typically created with interpolate. However, for customized control you may also construct them with

BSplineInterpolation(TWeights, coefs, it, axs)

where T gets computed from the product of TWeights and eltype(coefs). (This is equivalent to indicating that you'll be evaluating at locations itp(x::TWeights, y::TWeights, ...).)

source
Interpolations.BoundaryConditionType
BoundaryCondition

An abstract type with one of the following values (see the help for each for details):

  • Throw(gt)
  • Flat(gt)
  • Line(gt)
  • Free(gt)
  • Periodic(gt)
  • Reflect(gt)
  • InPlace(gt)
  • InPlaceQ(gt)

where gt is the grid type, e.g., OnGrid() or OnCell().

source
Interpolations.BoundsCheckStyleType
BoundsCheckStyle(itp)

A trait to determine dispatch of bounds-checking for itp. Can return NeedsCheck(), in which case bounds-checking is performed, or CheckWillPass() in which case the check will return true.

source
Interpolations.GriddedInterpolationMethod
GriddedInterpolation(typeOfWeights::Type{<:Real},
                     knots::NTuple{N, AbstractUnitRange },
                     array::AbstractArray{TCoefs,N},
                     interpolationType::DimSpec{<:Gridded})

Construct a GriddedInterpolation for generic knots from an AbstractUnitRange.

AbstractUnitRanges are collected to an Array to not confuse bound calculations (See Issue #398)
source
Interpolations.GriddedInterpolationMethod
GriddedInterpolation(typeOfWeights::Type{<:Real},
                     knots::NTuple{N, Union{ AbstractVector{T}, Tuple } },
                     array::AbstractArray{Tel,N},
                     interpolationType::DimSpec{<:Gridded})

Construct a GriddedInterpolation for generic knots from an AbstractArray
source
Interpolations.KnotIteratorType
KnotIterator{T,ET}(k::AbstractArray{T}, bc::ET)

Defines an iterator over the knots in k based on the boundary conditions bc.

Fields

  • knots::Vector{T} The interpolated knots of the axis to iterate over
  • bc::ET The Boundary Condition for the axis
  • nknots::Int The number of interpolated knots (ie. length(knots))

ET is Union{BoundaryCondition,Tuple{BoundaryCondition,BoundaryCondition}}

Iterator Interface

The following methods defining Julia's iterator interface have been defined

IteratorSize(::Type{KnotIterator}) -> Will return one of the following

  • Base.IsInfinite if the iteration will produces an infinite sequence of knots
  • Base.HasLength if iteration will produce a finite sequence of knots
  • Base.SizeUnknown if we can't decided from only the type information

length and size -> Are defined if IteratorSize is HasLength, otherwise will raise a MethodError.

IteratorEltype will always return HasEltype, as we always track the data types of the knots

eltype will return the data type of the knots

iterate Defines iteration over the knots starting from the first one and moving in the forward direction along the axis.

Knots for Multi-dimensional Interpolants

Iteration over the knots of a multi-dimensional interpolant is done by wrapping multiple KnotIterator within Iterators.product.

Indexing

KnotIterator provides limited support for accessing knots via indexing

  • getindex is provided for KnotIterator but does not support Multidimensional interpolations (As wrapped by ProductIterator) or non-Int indexes.
  • A BoundsError will be raised if out of bounds and checkbounds has been implemented for KnotIterator
julia> using Interpolations;

julia> etp = LinearInterpolation([1.0, 1.2, 2.3, 3.0], rand(4); extrapolation_bc=Periodic());

julia> kiter = knots(etp);

julia> kiter[4]
3.0

julia> kiter[36]
24.3
source
Interpolations.KnotRangeType
KnotRange(iter::KnotIterator{T}, start, stop)

Defines an iterator over a range of knots such that start < k < stop.

Fields

  • iter::KnotIterator{T} Underlying KnotIterator providing the knots iterated
  • range::R Iterator defining the range of knot indices iterated. Where R <: Union{Iterators.Count, UnitRange}

Iterator Interface

The following methods defining the Julia's iterator interface have been defined

Base.IteratorSize -> Will return one of the following:

  • Base.HasLength if range is of finite length
  • Base.IsInfinite if range is of infinite length
  • Base.SizeUnknown if the type of range is unspecified

Base.IteratorEltype -> Returns Base.EltypeUnknown if type parameter not provided, otherwise Base.HasEltype

length and size -> Returns the number of knots to be iterated if IteratorSize !== IsInfinite, otherwise will raise MethodError

Multidimensional Interpolants

Iteration over the knots of a multi-dimensional interpolant is done by wrapping multiple KnotRange iterators within Iterators.product.

source
Interpolations.WeightedIndexType
wi = WeightedIndex(indexes, weights)

Construct a weighted index wi, which can be thought of as a generalization of an ordinary array index to the context of interpolation. For an ordinary vector a, a[i] extracts the element at index i. When interpolating, one is typically interested in a range of indexes and the output is some weighted combination of array values at these indexes. For example, for linear interpolation at a location x between integers i and i+1, we have

ret = (1-f)*a[i] + f*a[i+1]

where f = x-i lies between 0 and 1. This can be represented as a[wi], where

wi = WeightedIndex(i:i+1, (1-f, f))

i.e.,

ret = sum(a[indexes] .* weights)

Linear interpolation thus constructs weighted indices using a 2-tuple for weights and a length-2 indexes range. Higher-order interpolation would involve more positions and weights (e.g., 3-tuples for quadratic interpolation, 4-tuples for cubic).

In multiple dimensions, separable interpolation schemes are implemented in terms of multiple weighted indices, accessing A[wi1, wi2, ...] where each wi is the WeightedIndex along the corresponding dimension.

For value interpolation, weights will typically sum to 1. However, for gradient and Hessian computation this will not necessarily be true. For example, the gradient of one-dimensional linear interpolation can be represented as

gwi = WeightedIndex(i:i+1, (-1, 1))
g1 = a[gwi]

For a three-dimensional array A, one might compute ∂A/∂x₂ (the second component of the gradient) as A[wi1, gwi2, wi3], where wi1 and wi3 are "value" weights and gwi2 "gradient" weights.

indexes may be supplied as a range or as a tuple of the same length as weights. The latter is applicable, e.g., for periodic boundary conditions.

source