Interpolations.jl

Interpolations.jl provides high-performance interpolation schemes for Julia.

Supported features include:

• B-splines
• Irregular grids
• Lanczos interpolation

The package is designed for ease-of-use and broad algorithmic support, with an API intended to evolve as more methods are added.

A WeightedIndex is an abstract type used to perform interpolation by combining an index and a weight. There are two concrete subtypes:

• WeightedAdjIndex: Used for indices addressing adjacent points (where the index increments by 1). This is used for schemes like Linear interpolation or when prefiltering provides sufficient padding.
• WeightedArbIndex: Stores both the weight and the index for each accessed grid point. This is used when access patterns are non-adjacent, such as with periodic boundary conditions.

Index transformation extrapolation handles coordinates outside the domain by mapping them back into the domain using a specific mathematical rule. Once the coordinate is transformed to an in-bounds index, the standard interpolation process proceeds.

Common index transformation strategies include:

• ExtrapPeriodic: Uses modulo calculations to wrap coordinates back into the domain (useful for periodic/cyclic data).
• ExtrapConstant: Clamps out-of-bounds coordinates to the nearest available in-bounds data point.

Interpolant construction in Interpolations.jl involves creating an object that stores the coefficient array and the interpolation scheme.

• Simple Schemes: For NoInterp, Constant, or Linear interpolation, construction simply records the array and settings. Values are computed directly from on-grid values.
• Higher-Order Schemes: For Quadratic and higher orders, the array of values must be prefiltered. This process 'inverts' the interpolation computation to approximately reconstruct original values at on-grid points, typically by solving a nearly-tridiagonal system of equations. This prefiltering is implemented independently along each axis using the AxisAlgorithms package.

When using directional boundary conditions (e.g., (Throw(), Periodic())), the iterator's finiteness is determined by the forward boundary condition.

• If the forward condition is infinite (e.g., Periodic), knots(etp) produces an infinite sequence. You can check this using Base.IteratorSize(kiter) == Base.IsInfinite().
• If the forward condition is finite (e.g., Throw), the sequence is finite, and Base.IteratorSize(kiter) == Base.HasLength().

x = [1.0, 1.2, 1.3, 2.0];

# Unbounded (Infinite) sequence
etp = linear_interpolation(x, x.^2, extrapolation_bc=((Throw(), Periodic()),));
kiter = knots(etp);
Base.IteratorSize(kiter) # Base.IsInfinite()

# Bounded (Finite) sequence
etp = linear_interpolation(x, x.^2, extrapolation_bc=((Periodic(), Throw()),));
kiter = knots(etp);
length(kiter) # 4
Base.IteratorSize(kiter) # Base.HasLength()

Interpolant usage (evaluating itp(x, y...) or Interpolations.gradient(itp, x, y...)) consists of two steps:

1. Weight Computation: Calculating weights based on the interpolation scheme and the location x, y.... These weights are independent of the coefficient array, allowing for reuse.
2. Interpolation: Performing the actual interpolation using the computed weights.

Internally, this is achieved by indexing the coefficient array with WeightedIndex objects. For example, itp(x...) is equivalent to InterpGetindex(itp.coefs)[wis...], where wis is a tuple of WeightedIndex indices.

The interpolation degree determines the continuity properties of the resulting interpolant.

• Linear interpolation: A piecewise linear function that is continuous but has discontinuous derivatives (gradients).
• Higher degrees: Higher degrees provide smoother transitions (higher-order continuity).

In Interpolations.jl, the degree is referred to by names like linear, quadratic, etc., to avoid confusion with B-spline order terminology.

For interpolation degrees of quadratic and higher, the B-spline support is too large to solve the equations at the edges of the data set without additional information. Interpolations.jl uses boundary conditions to close these equation systems.

Common boundary condition strategies include:

• Flat boundary conditions: Assuming the derivative is 0 at the edges.
• Extension: Extending the second-to-outermost value all the way to the edge.

Boundary conditions can be applied at different locations:

• OnGrid(): Applied at the edge grid point.
• OnCell(): Applied at the halfway mark to the first beyond-the-edge index.

Note: Constant and linear interpolations do not require boundary conditions and therefore do not support them.

Standard interpolation methods can produce non-monotonic results (overshooting) even when the input data is monotonic. Use a monotonic interpolation algorithm to ensure the interpolating function preserves the monotonicity of the data.

Available Algorithms

# Example: Creating a monotonic CDF
y = sort(randn(100))
x = range(0, 1, length=100)

itp_cdf = extrapolate(
    interpolate(y, x, SteffenMonotonicInterpolation()),
    Flat()
)

Interpolations.jl is highly optimized for regular grids with uniform spacing.

• Highest Performance: Use knots that are an AbstractUnitRange (e.g., 2:5 or Base.OneTo(9)). If no knots are specified, the package defaults to a UnitRange starting at 1.
• Scaling Penalty: If knots are not unit-spaced or do not start at 1, you must use the scale function. This provides flexibility but introduces a performance penalty.
• Irregular Grids: For irregularly spaced knots, the ranges between knots must be scaled (similar to Gridded interpolation types).