Type Reference

AbstractInterpolant{Tg<:AbstractFloat, Tv}

Abstract supertype for all interpolant objects.

Type Parameters

• Tg: Grid/coordinate type (Float32 or Float64) - used for x-coordinates, spacing, search
• Tv: Value type - used for y-values, coefficients, return values. Can be Tg (real-valued), Complex{Tg} (complex-valued), or other Number types.

Design Invariant

• Grid operations (search, spacing) always use Tg
• Value operations (kernel, coefficients) use Tv
• Evaluation returns type based on promotetype(Tv, querytype)

Subtypes

• AbstractSeriesInterpolant{Tg, Tv}: Multi-series interpolants
• LinearInterpolant{Tg, Tv}: Piecewise linear interpolation
• ConstantInterpolant{Tg, Tv}: Piecewise constant (step) interpolation
• QuadraticInterpolant{Tg, Tv}: C1 piecewise quadratic spline
• CubicInterpolant{Tg, Tv}: C2 zero-curvature/zero-slope/periodic cubic spline

Note

This is a pure type hierarchy - no methods are defined on AbstractInterpolant itself. All functionality is implemented in concrete subtypes.

AbstractSeriesInterpolant{Tg<:AbstractFloat, Tv}

Abstract supertype for multi-series interpolant objects. Series interpolants handle multiple y-series sharing the same x-grid.

Type Parameters

• Tg: Grid/coordinate type (Float32 or Float64)
• Tv: Value type - can be real or complex

Subtypes

• LinearSeriesInterpolant{Tg, Tv}: Multiple linear interpolants sharing x-grid
• ConstantSeriesInterpolant{Tg, Tv}: Multiple constant interpolants sharing x-grid
• QuadraticSeriesInterpolant{Tg, Tv}: Multiple quadratic interpolants sharing x-grid
• CubicSeriesInterpolant{Tg, Tv}: Multiple cubic splines sharing x-grid

Key Features

• Anchor optimization: compute interval once, evaluate all series
• Matrix storage: unified storage for optimal SIMD vectorization
• Zero-allocation batch evaluation with pre-built anchors

Usage

x = collect(range(0.0, 1.0, 101))
y1, y2, y3 = sin.(2π .* x), cos.(2π .* x), exp.(-x)

sitp = cubic_interp(x, [y1, y2, y3])  # Creates CubicSeriesInterpolant

vals = sitp(0.5)            # Returns Vector of 3 values
sitp(output, 0.5)           # In-place evaluation

Note

This is a pure type hierarchy - no methods are defined on AbstractSeriesInterpolant itself. All functionality is implemented in concrete subtypes.

AbstractInterpolantND{Tg<:AbstractFloat, Tv, N}

Abstract supertype for N-dimensional interpolant objects.

Type Parameters

• Tg: Grid/coordinate type (Float32 or Float64)
• Tv: Value type - can be real, complex, or other Number types
• N: Number of dimensions

API Differences from 1D Interpolants

• Evaluation: itp(x::NTuple{N}) or itp(x::AbstractVector) instead of itp(x::Real)
• Derivatives: Use deriv keyword (e.g., itp(x; deriv=(1,0))) or deriv_view(itp, (1,0))
• Vector Calculus: Supports gradient, hessian, laplacian

Subtypes

• CubicInterpolantND{Tg, Tv, N}: N-dimensional cubic Hermite interpolation

Example

x, y = range(0, 1, 50), range(0, 1, 50)
data = [sin(xi) * cos(yj) for xi in x, yj in y]
itp = cubic_interp((x, y), data)  # Returns CubicInterpolantND{..., 2}

itp((0.5, 0.5))                    # Evaluate
itp((0.5, 0.5); deriv=(1, 0))      # ∂f/∂x
gradient(itp, (0.5, 0.5))          # (∂f/∂x, ∂f/∂y)

AbstractCoeffStrategy

Abstract supertype for coefficient computation strategies in ND interpolation.

Implemented Strategies

• PreCompute: Precompute all partial derivatives at construction (O(1) query)
• OnTheFly: Compute coefficients lazily at query time (O(n) query)

PreCompute <: AbstractCoeffStrategy

Precompute all partial derivatives at construction time.

For N-dimensional interpolation, stores 2^N partial derivatives per grid point:

• 2D: 4 partials (f, ∂f/∂x, ∂f/∂y, ∂²f/∂x∂y)
• 3D: 8 partials (f, ∂f/∂x, ∂f/∂y, ∂f/∂z, ∂²f/∂x∂y, ∂²f/∂x∂z, ∂²f/∂y∂z, ∂³f/∂x∂y∂z)

Trade-offs

• Memory: O(2^N × n^N) - higher than OnTheFly
• Construction: O(N × n^N) - expensive (N passes of 1D spline solving)
• Query: O(1) - ultra-fast (just polynomial evaluation)

Best for

• Static data with many queries
• Real-time applications requiring fast evaluation
• Scenarios where construction cost can be amortized

OnTheFly <: AbstractCoeffStrategy

Compute coefficients lazily at query time using tensor-product (separable) approach.

Stores only the original data; computes 1D spline solutions per-axis at each query.

Trade-offs

• Memory: O(n^N) - minimal (only original data)
• Construction: O(1) - instant (just store data reference)
• Query: O(n) per axis - slower (must solve 1D systems)

Best for

• Few queries on large datasets
• Memory-constrained environments
• Data that changes frequently

CubicInterpolantND{Tg, Tv, N, NP1, G, S, B, E, P}

Generic N-dimensional cubic Hermite interpolant with precomputed partial derivatives.

Stores function values AND all partial derivatives at grid nodes, enabling ultra-fast O(1) evaluation via tensor-product Hermite polynomials.

Type Parameters

• Tg: Grid/coordinate type (Float32 or Float64)
• Tv: Value type (can be Tg, Complex{Tg}, or other Number)
• N: Number of dimensions
• NP1: N + 1 (partials array dimensionality)
• G: Tuple type for grids, NTuple{N, <:AbstractVector{Tg}}
• S: Tuple type for spacings, NTuple{N, <:AbstractGridSpacing{Tg}}
• B: Tuple type for boundary conditions, NTuple{N, <:AbstractBC}
• E: Tuple type for extrapolation modes, Tuple{Vararg{AbstractExtrap, N}}
• P: Tuple type for search policies, NTuple{N, <:AbstractSearchPolicy}

Fields

• grids: N-tuple of grid vectors for each dimension
• spacings: N-tuple of grid spacing info (for O(1) h lookup)
• nodal_derivs: NodalDerivativesND containing partial derivatives at grid nodes
• bcs: N-tuple of boundary conditions used for construction
• extraps: N-tuple of extrapolation modes
• searches: N-tuple of search policies

Performance

• Construction: O(2^N × n^N) - computes all partial derivatives
• Query: O(1) - tensor-product Hermite polynomial evaluation
• Memory: 2^N × n^N values

Thread-Safety

Immutable after construction; safe for concurrent read access.

Example

x = range(0.0, 2π, 50)
y = range(0.0, π, 30)
z = range(0.0, 1.0, 20)
data = [sin(xi) * cos(yj) * zk for xi in x, yj in y, zk in z]

itp = cubic_interp((x, y, z), data)  # Returns CubicInterpolantND{..., 3, 4, ...}
itp((1.0, 0.5, 0.3))                  # Evaluate at (1.0, 0.5, 0.3)

grid_type(::AbstractInterpolant{Tg, Tv}) -> Type{Tg}

Get the grid/coordinate type of an interpolant.

value_type(::AbstractInterpolant{Tg, Tv}) -> Type{Tv}

Get the value type of an interpolant.

eval_type(::AbstractInterpolant{Tg, Tv}, ::Type{Tq}) -> Type

Compute the output type when evaluating an interpolant with query type Tq. This is promote_type(Tv, Tq), accounting for value type (real or complex) and query type (standard float or ForwardDiff.Dual).

Examples

itp = linear_interp([0.0, 1.0], [1.0, 2.0])
eval_type(itp, Float64)  # Float64

itp_c = linear_interp([0.0, 1.0], [1.0+0im, 2.0+0im])
eval_type(itp_c, Float64)  # ComplexF64

AbstractEvalOp

Abstract type for evaluation operations (value, derivatives). Used for compile-time dispatch to select appropriate kernel.

Subtypes:

• DerivOp{0} (alias EvalValue): Evaluate function value f(x)
• DerivOp{1} (alias EvalDeriv1): Evaluate first derivative f'(x)
• DerivOp{2} (alias EvalDeriv2): Evaluate second derivative f''(x)
• DerivOp{3} (alias EvalDeriv3): Evaluate third derivative f'''(x)

DerivOp{N} <: AbstractEvalOp

Parametric singleton for derivative order dispatch. N is the derivative order (0 = value, 1 = first derivative, etc.).

Construction

• DerivOp{0}(), DerivOp{1}() — direct parametric construction
• DerivOp(n::Int) — convenience: DerivOp(1) → DerivOp{1}()
• DerivOp(n1, n2, ...) — ND: DerivOp(1, 0) → (DerivOp{1}(), DerivOp{0}())

Backward-compatible aliases

• EvalValue = DerivOp{0}
• EvalDeriv1 = DerivOp{1}
• EvalDeriv2 = DerivOp{2}
• EvalDeriv3 = DerivOp{3}

Examples

# 1D scalar evaluation
itp(x; deriv=DerivOp(1))        # first derivative
itp(x; deriv=EvalDeriv2())      # second derivative (using alias)

# ND mixed partials
itp(q; deriv=DerivOp(1, 0))     # ∂f/∂x
itp(q; deriv=DerivOp(0, 2))     # ∂²f/∂y²

deriv_order(::DerivOp{N}) -> Int

Extract derivative order from a DerivOp instance.

DerivOp{N} <: AbstractEvalOp

Parametric singleton for derivative order dispatch. N is the derivative order (0 = value, 1 = first derivative, etc.).

Construction

• DerivOp{0}(), DerivOp{1}() — direct parametric construction
• DerivOp(n::Int) — convenience: DerivOp(1) → DerivOp{1}()
• DerivOp(n1, n2, ...) — ND: DerivOp(1, 0) → (DerivOp{1}(), DerivOp{0}())

Backward-compatible aliases

• EvalValue = DerivOp{0}
• EvalDeriv1 = DerivOp{1}
• EvalDeriv2 = DerivOp{2}
• EvalDeriv3 = DerivOp{3}

Examples

# 1D scalar evaluation
itp(x; deriv=DerivOp(1))        # first derivative
itp(x; deriv=EvalDeriv2())      # second derivative (using alias)

# ND mixed partials
itp(q; deriv=DerivOp(1, 0))     # ∂f/∂x
itp(q; deriv=DerivOp(0, 2))     # ∂²f/∂y²

DerivOp{N} <: AbstractEvalOp

Parametric singleton for derivative order dispatch. N is the derivative order (0 = value, 1 = first derivative, etc.).

Construction

• DerivOp{0}(), DerivOp{1}() — direct parametric construction
• DerivOp(n::Int) — convenience: DerivOp(1) → DerivOp{1}()
• DerivOp(n1, n2, ...) — ND: DerivOp(1, 0) → (DerivOp{1}(), DerivOp{0}())

Backward-compatible aliases

• EvalValue = DerivOp{0}
• EvalDeriv1 = DerivOp{1}
• EvalDeriv2 = DerivOp{2}
• EvalDeriv3 = DerivOp{3}

Examples

# 1D scalar evaluation
itp(x; deriv=DerivOp(1))        # first derivative
itp(x; deriv=EvalDeriv2())      # second derivative (using alias)

# ND mixed partials
itp(q; deriv=DerivOp(1, 0))     # ∂f/∂x
itp(q; deriv=DerivOp(0, 2))     # ∂²f/∂y²

DerivOp{N} <: AbstractEvalOp

Parametric singleton for derivative order dispatch. N is the derivative order (0 = value, 1 = first derivative, etc.).

Construction

• DerivOp{0}(), DerivOp{1}() — direct parametric construction
• DerivOp(n::Int) — convenience: DerivOp(1) → DerivOp{1}()
• DerivOp(n1, n2, ...) — ND: DerivOp(1, 0) → (DerivOp{1}(), DerivOp{0}())

Backward-compatible aliases

• EvalValue = DerivOp{0}
• EvalDeriv1 = DerivOp{1}
• EvalDeriv2 = DerivOp{2}
• EvalDeriv3 = DerivOp{3}

Examples

# 1D scalar evaluation
itp(x; deriv=DerivOp(1))        # first derivative
itp(x; deriv=EvalDeriv2())      # second derivative (using alias)

# ND mixed partials
itp(q; deriv=DerivOp(1, 0))     # ∂f/∂x
itp(q; deriv=DerivOp(0, 2))     # ∂²f/∂y²

DerivOp{N} <: AbstractEvalOp

Parametric singleton for derivative order dispatch. N is the derivative order (0 = value, 1 = first derivative, etc.).

Construction

• DerivOp{0}(), DerivOp{1}() — direct parametric construction
• DerivOp(n::Int) — convenience: DerivOp(1) → DerivOp{1}()
• DerivOp(n1, n2, ...) — ND: DerivOp(1, 0) → (DerivOp{1}(), DerivOp{0}())

Backward-compatible aliases

• EvalValue = DerivOp{0}
• EvalDeriv1 = DerivOp{1}
• EvalDeriv2 = DerivOp{2}
• EvalDeriv3 = DerivOp{3}

Examples

# 1D scalar evaluation
itp(x; deriv=DerivOp(1))        # first derivative
itp(x; deriv=EvalDeriv2())      # second derivative (using alias)

# ND mixed partials
itp(q; deriv=DerivOp(1, 0))     # ∂f/∂x
itp(q; deriv=DerivOp(0, 2))     # ∂²f/∂y²

coeffs(itp, x; search=..., hint=...) → CellPoly

Extract the local polynomial coefficients of itp at query point x. Returns a CellPoly with ascending-power coefficients for evalpoly.

The returned CellPoly is callable: cell(x) evaluates the polynomial at x. Alternatively, use evalpoly(x, cell) for the same result.

Cubic

For CubicInterpolant, returns CellPoly{4} with coefficients (d, c, b, a) representing S(u) = d + c·u + b·u² + a·u³ where u = x - xL.

Quadratic

For QuadraticInterpolant, returns CellPoly{3} with coefficients (y₀, d, a) representing S(u) = y₀ + d·u + a·u².

Linear

For LinearInterpolant, returns CellPoly{2} with coefficients (y₀, slope) representing S(u) = y₀ + slope·u.

Constant

For ConstantInterpolant, returns CellPoly{1} with coefficient (y₀,).

Example

itp = cubic_interp([0.0, 1.0, 2.0, 3.0], [0.0, 1.0, 4.0, 9.0])
cell = coeffs(itp, 1.5)
cell(1.5)            # same as itp(1.5)
evalpoly(1.5, cell)  # same result
cell.p               # (d, c, b, a) in ascending power order

CellPoly{N, Tv, Tg}

Lightweight immutable polynomial on a spline cell. Coefficients in ascending power order for evalpoly compatibility.

Type Parameters

• N: Number of coefficients (= polynomial degree + 1)
• Tv: Value type (Float64, ComplexF64, etc.)
• Tg: Grid coordinate type (Float64, Float32)

Fields

• p::NTuple{N, Tv}: Polynomial coefficients in ascending power order
• xL::Tg: Left endpoint of cell
• xR::Tg: Right endpoint of cell

Polynomial form

S(x) = p₁ + p₂·u + p₃·u² + ... + pₙ·uⁿ⁻¹    where u = x - xL

Usage

cell = coeffs(itp, 1.5)
cell(1.5)                   # callable: evaluates polynomial at global x
evalpoly(1.5, cell)          # overloaded: auto-computes u = x - xL
cell.p                       # raw coefficients for custom operations
cell.xL, cell.xR             # cell boundaries

gradient(itp::AbstractInterpolantND, query)

Compute the gradient (vector of partial derivatives) at query.

Returns an NTuple{N} of partial derivatives (∂f/∂x₁, ∂f/∂x₂, ..., ∂f/∂xₙ).

Performance

~9x faster than ForwardDiff.gradient by using analytical derivatives. Uses locate-once optimization: interval search performed only once per query point.

Examples

itp = cubic_interp((x, y), data)
gradient(itp, (0.5, 0.5))    # → (∂f/∂x, ∂f/∂y)
gradient(itp, [0.5, 0.5])    # Vector input also supported

See also: gradient!, hessian, laplacian

gradient!(G, itp::AbstractInterpolantND, query)

Compute the gradient in-place, writing partial derivatives into G.

Zero-allocation version of gradient for use in optimization loops.

Examples

itp = cubic_interp((x, y), data)
G = zeros(2)
gradient!(G, itp, (0.5, 0.5))    # G .= (∂f/∂x, ∂f/∂y)
gradient!(G, itp, [0.5, 0.5])    # G .= (∂f/∂x, ∂f/∂y)

# Optim.jl compatible:
grad!(G, x) = gradient!(G, itp, x)
result = optimize(f, grad!, x0, LBFGS())

See also: gradient, hessian!

hessian(itp::AbstractInterpolantND, query)

Compute the Hessian matrix (matrix of second partial derivatives) at query.

Returns an N×N matrix where H[i,j] = ∂²f/∂xᵢ∂xⱼ.

Performance

~9x faster than ForwardDiff.hessian by using analytical derivatives. Exploits symmetry: computes only N(N+1)/2 unique elements. Uses locate-once optimization: interval search performed only once per query point.

Examples

itp = cubic_interp((x, y), data)
H = hessian(itp, (0.5, 0.5))
# H = [∂²f/∂x²    ∂²f/∂x∂y]
#     [∂²f/∂x∂y   ∂²f/∂y² ]

See also: gradient, hessian!, laplacian

hessian!(H, itp::AbstractInterpolantND, query)

Compute the Hessian matrix in-place, writing second partial derivatives into H.

Zero-allocation version of hessian for use in optimization loops. Exploits symmetry: computes only N(N+1)/2 unique elements.

Examples

itp = cubic_interp((x, y), data)
H = zeros(2, 2)
hessian!(H, itp, (0.5, 0.5))

# Optim.jl compatible:
hess!(H, x) = hessian!(H, itp, x)
result = optimize(f, grad!, hess!, x0, NewtonTrustRegion())

See also: hessian, gradient!

laplacian(itp::AbstractInterpolantND, query)

Compute the Laplacian (sum of second partial derivatives) at query.

Returns a scalar: ∇²f = ∂²f/∂x₁² + ∂²f/∂x₂² + ... + ∂²f/∂xₙ²

Performance

Faster than computing full Hessian when only the trace is needed. Uses locate-once optimization: interval search performed only once per query point.

Examples

itp = cubic_interp((x, y), data)
∇²f = laplacian(itp, (0.5, 0.5))  # ∂²f/∂x² + ∂²f/∂y²

# Equivalent to (but faster than):
# tr(hessian(itp, (0.5, 0.5)))

Applications

• Heat equation: ∂T/∂t = α∇²T
• Wave equation: ∂²u/∂t² = c²∇²u
• Poisson equation: ∇²φ = ρ

See also: gradient, hessian

AbstractBC

Abstract base type for all boundary condition specifications. Type-free design: no type parameter on the abstract type. Concrete subtypes with values (Deriv1, Deriv2, Deriv3) carry their own Tv parameter.

Subtypes

• ZeroCurvBC: Zero-Curvature BC (zero curvature at both ends) - singleton
• ZeroSlopeBC: Zero-Slope BC (zero slope at both ends) - singleton
• PeriodicBC{E,P}: Periodic boundary condition (inclusive/exclusive endpoint)
• PointBC: Single-point derivative conditions (abstract)
• BCPair{L,R}: Pair of left/right boundary conditions
• Left{B}: Endpoint wrapper for BC at left (x[1]) - used by quadratic splines
• Right{B}: Endpoint wrapper for BC at right (x[end]) - used by quadratic splines

PointBC <: AbstractBC

Abstract type for single-point boundary conditions. Represents a derivative condition at one endpoint. Type-free design: concrete subtypes carry their own Tv parameter.

Subtypes

• Deriv1{Tv}: First derivative (slope) BC
• Deriv2{Tv}: Second derivative (curvature) BC
• Deriv3{Tv}: Third derivative BC
• PolyFit{D}: Polynomial fit (lazy, no Tv until materialization)

Deriv1{Tv} <: PointBC

First derivative (slope) boundary condition: S'(endpoint) = val

The type parameter Tv is the value type, supporting both Real and Complex.

Example

Deriv1(0.5)           # Slope of 0.5 at endpoint (Float64)
Deriv1(0)             # Zero slope (horizontal tangent)
Deriv1(1.0+2.0im)     # Complex slope (ComplexF64)

Deriv2{Tv} <: PointBC

Second derivative (curvature) boundary condition: S''(endpoint) = val

The type parameter Tv is the value type, supporting both Real and Complex.

Example

Deriv2(0)             # Zero curvature (same as ZeroCurvBC)
Deriv2(1.5)           # Specified curvature at endpoint
Deriv2(0.0+0.0im)     # Complex curvature

Deriv3{Tv} <: PointBC

Third derivative boundary condition: S'''(endpoint) = val

For cubic splines, the third derivative is constant within each interval: S'''(x) = (z[i+1] - z[i]) / h[i]. This BC specifies the third derivative value at the first (or last) interval.

The type parameter Tv is the value type, supporting both Real and Complex.

Example

Deriv3(0)             # Zero third derivative at endpoint
Deriv3(1.0)           # Specified third derivative
Deriv3(0.0+0.0im)     # Complex third derivative

BCPair{L<:PointBC, R<:PointBC} <: AbstractBC

Container for left and right boundary conditions with type parameters for zero-overhead dispatch. The BC types are encoded in the type parameters, enabling compile-time specialization.

Type-free design: no Tv parameter on BCPair itself. The value type is determined by the concrete L and R types (e.g., BCPair{Deriv1{Float64}, Deriv2{Float64}}).

Type Parameters

• L: Left boundary condition type (Deriv1{Tv}, Deriv2{Tv}, PolyFit{D}, etc.)
• R: Right boundary condition type (Deriv1{Tv}, Deriv2{Tv}, PolyFit{D}, etc.)

Example

bc = BCPair(Deriv1(0.5), Deriv2(0))           # Left: slope=0.5, Right: zero curvature
bc = BCPair(Deriv1(1.0+0.0im), Deriv2(0.0im)) # Complex BC
bc = BCPair(CubicFit(), Deriv2(0.0))          # Mixed: lazy left, concrete right

ZeroCurvBC <: AbstractBC

Zero-curvature boundary condition: S''(endpoints) = 0 (zero curvature at both ends). Equivalent to BCPair(Deriv2(0), Deriv2(0)).

This is a singleton type (structure-only). Normalized to BCPair with Deriv2(zero(Tv)) at both ends based on the value type during interpolant construction.

Example

itp = cubic_interp(x, y; bc=ZeroCurvBC())  # Explicit zero-curvature BC

ZeroSlopeBC <: AbstractBC

Zero-slope boundary condition: S'(endpoints) = 0 (zero slope at both ends). Equivalent to BCPair(Deriv1(0), Deriv1(0)).

This is a singleton type (structure-only). Normalized to BCPair with Deriv1(zero(Tv)) at both ends based on the value type during interpolant construction.

Example

itp = cubic_interp(x, y; bc=ZeroSlopeBC())

PeriodicBC{E, P} <: AbstractBC

Periodic boundary condition: S(x0) = S(xn), S'(x0) = S'(xn), S''(x0) = S''(xn)

Internally, periodic BC uses Sherman-Morrison solver with PeriodicData{T} for the cache.

Type Parameters

• E::Symbol: :inclusive or :exclusive (compile-time endpoint convention)
• P: Nothing (inclusive or auto-infer) or <:AbstractFloat (explicit period)

Endpoint Conventions

• Inclusive (endpoint=:inclusive, default): y[1] ≈ y[end] required (standard convention)
• Exclusive (endpoint=:exclusive): y[end] is the last unique data point; the period boundary is handled internally. For AbstractRange grids, the period is inferred from step(x) * length(x). For non-uniform grids, period must be provided explicitly.

Examples

# Inclusive (default, backward compatible)
cache = CubicSplineCache(x; bc=PeriodicBC())

# Exclusive with explicit period
itp = cubic_interp(x, y; bc=PeriodicBC(endpoint=:exclusive, period=2π))

# Exclusive with Range grid (period auto-inferred)
x = range(0, step=2π/64, length=64)
itp = cubic_interp(x, sin.(x); bc=PeriodicBC(endpoint=:exclusive))

PolyFit{D} <: PointBC

Generic polynomial fitting boundary condition (lazy).

Fits a degree-D polynomial through (D+1) points at the endpoint and evaluates its derivative. This automatically estimates the first derivative from data without requiring user-specified values.

This is a lazy type with no value parameter. The actual derivative value is computed during materialize_bc() based on the data's value type (Tv).

Type Parameters

• D::Int: Polynomial degree (1=linear, 2=quadratic, 3=cubic, ...)

Relationships

Points needed = D + 1
Accuracy = O(h^D) for smooth functions

Aliases (Recommended for common cases)

• LinearFit = PolyFit{1} → 2 points, O(h)
• QuadraticFit = PolyFit{2} → 3 points, O(h²)
• CubicFit = PolyFit{3} → 4 points, O(h³)

Mathematical Background

All polynomial fitting methods are mathematically equivalent to finite difference formulas of the same order. For example, CubicFit (4-point) gives identical coefficients to 4-point one-sided finite difference:

• Left: f'(x₁) ≈ (-11f₁ + 18f₂ - 9f₃ + 2f₄) / (6h)
• Right: f'(xₙ) ≈ (-2fₙ₋₃ + 9fₙ₋₂ - 18fₙ₋₁ + 11fₙ) / (6h)

Key Property

Polynomial Reproduction: PolyFit{D} exactly reproduces polynomials up to degree D.

Example

# Recommended: use named aliases
itp = cubic_interp(x, y; bc=CubicFit())
itp = quadratic_interp(x, y; bc=Left(QuadraticFit()))

# Generic form (equivalent)
itp = cubic_interp(x, y; bc=PolyFit{3}())

See also: LinearFit, QuadraticFit, CubicFit, Deriv1

LinearFit = PolyFit{1}

2-point linear fit boundary condition. Estimates derivative using forward/backward difference: f'(x) ≈ (f₂ - f₁) / h.

Accuracy: O(h) - first order. Points needed: 2

Example

itp = cubic_interp(x, y; bc=LinearFit())

QuadraticFit = PolyFit{2}

3-point quadratic fit boundary condition. Fits a parabola through 3 points and evaluates its derivative at the endpoint.

Accuracy: O(h²) - second order. Points needed: 3

For uniform grids: f'(x₁) ≈ (-3f₁ + 4f₂ - f₃) / (2h)

Key Property

Polynomial Reproduction: Exactly reproduces quadratic polynomials.

Example

# Default BC for quadratic splines
itp = quadratic_interp(x, y; bc=Left(QuadraticFit()))

CubicFit = PolyFit{3}

4-point cubic fit boundary condition. Fits a cubic polynomial through 4 points and evaluates its derivative at the endpoint.

Accuracy: O(h³) - third order. Points needed: 4

For uniform grids:

• Left: f'(x₁) ≈ (-11f₁ + 18f₂ - 9f₃ + 2f₄) / (6h)
• Right: f'(xₙ) ≈ (-2fₙ₋₃ + 9fₙ₋₂ - 18fₙ₋₁ + 11fₙ) / (6h)

Key Property

Polynomial Reproduction: Exactly reproduces cubic polynomials.

Example

# Recommended BC for cubic splines when derivative is unknown
itp = cubic_interp(x, y; bc=CubicFit())

# Mixed with other BCs
itp = cubic_interp(x, y; bc=BCPair(CubicFit(), Deriv2(0)))

Left{B<:PointBC} <: AbstractBC

Wrapper indicating BC is applied at left endpoint (x[1]). Used for quadratic splines where only one endpoint BC is specified.

The type parameter B is the inner PointBC type, which carries the value type Tv for concrete types (Deriv1{Tv}, etc.) or just the degree D for lazy types (PolyFit{D}).

Example

bc = Left(Deriv1(0.5))        # slope = 0.5 at left endpoint
bc = Left(Deriv2(0.0))        # curvature = 0 at left endpoint
bc = Left(Deriv1(1.0+2.0im))  # Complex slope
bc = Left(QuadraticFit())     # Lazy polynomial fit

Right{B<:PointBC} <: AbstractBC

Wrapper indicating BC is applied at right endpoint (x[end]). Used for quadratic splines where only one endpoint BC is specified.

The type parameter B is the inner PointBC type, which carries the value type Tv for concrete types (Deriv1{Tv}, etc.) or just the degree D for lazy types (PolyFit{D}).

Example

bc = Right(Deriv1(2.0))        # slope = 2.0 at right endpoint
bc = Right(Deriv2(0.0))        # curvature = 0 at right endpoint
bc = Right(Deriv1(1.0+2.0im))  # Complex slope
bc = Right(QuadraticFit())     # Lazy polynomial fit

AbstractExtrap

Abstract type for typed extrapolation mode specification. Use concrete subtypes at the API boundary for compile-time dispatch.

Concrete subtypes

• NoExtrap: Throw DomainError for out-of-domain queries
• ConstExtrap: Clamp to nearest boundary value
• ExtendExtrap: Extend interpolation polynomial beyond domain
• WrapExtrap: Wrap queries into domain (periodic)

NoExtrap <: AbstractExtrap

No extrapolation — throws DomainError for out-of-domain queries.

Example

itp = cubic_interp((x, y), data; extrap=NoExtrap())

ConstExtrap <: AbstractExtrap

Constant extrapolation — clamps queries to the nearest boundary value.

Example

itp = cubic_interp((x, y), data; extrap=ConstExtrap())

ExtendExtrap <: AbstractExtrap

Extension extrapolation — extends the interpolation polynomial beyond the domain.

Example

itp = cubic_interp((x, y), data; extrap=ExtendExtrap())

WrapExtrap <: AbstractExtrap

Wrap extrapolation — wraps queries into the domain using modular arithmetic. For periodic data.

Example

itp = cubic_interp((x, y), data; extrap=WrapExtrap())

AbstractSearchPolicy

Abstract supertype for user-facing search policy selection. Concrete subtypes encode the search algorithm choice.

See also: BinarySearch, LinearBinarySearch

AutoSearch <: AbstractSearchPolicy

Adaptive search policy that resolves at call time based on query type:

• Scalar queries (Real, Tuple{Vararg{Real}}) → BinarySearch() — no hint locality to exploit
• Vector queries (AbstractVector, Tuple{Vararg{AbstractVector}}) → LinearBinarySearch() — hint continuity benefits sorted sequences
• Broadcast (itp.(xs)) → BinarySearch() per element — fresh searcher each call

This is the default search policy for all interpolants. For known query patterns, specify BinarySearch() or LinearBinarySearch() explicitly for optimal performance.

Example

itp = linear_interp(x, y)              # stores AutoSearch()
itp(0.5)                               # → BinarySearch() internally (scalar)
itp([0.1, 0.5, 0.9])                   # → LinearBinarySearch() internally (vector)
itp(0.5; search=LinearBinarySearch())        # override: use LinearBinarySearch explicitly

See also: BinarySearch, LinearBinarySearch

BinarySearch <: AbstractSearchPolicy

Binary search algorithm. O(log n) per query for vectors, O(1) for ranges. Stateless and thread-safe. This is the default search policy.

Hint Behavior

When a hint argument is provided with BinarySearch(), the search automatically upgrades to LinearBinarySearch() (default window) to utilize the hint for locality. Without a hint, pure binary search is used.

Example

val = linear_interp(x, y, 0.5; search=BinarySearch())  # explicit binary search
val = linear_interp(x, y, 0.5)                    # default: AutoSearch()

# With hint: auto-upgrades to LinearBinarySearch() (default window)
hint = Ref(1)
val = itp(0.5; search=BinarySearch(), hint=hint)  # uses LinearBinarySearch() internally

See also: LinearBinarySearch

LinearSearch <: AbstractSearchPolicy

Pure linear search for strictly monotonic query sequences. Maximum speed with no binary fallback and minimal bounds checking.

Safety Contract (User Responsibility)

Preconditions that MUST be satisfied:

1. Queries must be monotonically ordered (all increasing OR all decreasing)
2. All queries must satisfy x[1] <= xi <= x[end] (within interpolation domain)

If violated: Undefined behavior (infinite loop or out-of-bounds access)

Performance Characteristics

• Best case: O(1) per query (consecutive intervals)
• Amortized: O(1) for properly monotonic sequences
• Worst case: O(n) if sequence is not monotonic

When to Use

• ODE integration with strictly monotonic time stepping
• Streaming data evaluation with guaranteed ordering
• Performance-critical loops where ALL inputs are controlled

When NOT to Use

• Random access patterns (use BinarySearch)
• Queries that might exceed domain bounds (use LinearBinarySearch)
• Untrusted input data (use LinearBinarySearch)

Example

# ODE-style monotonic evaluation
hint = Ref(1)
for t in t_values  # strictly increasing
    y = itp(t; search=LinearSearch(), hint=hint)
end

See also: LinearBinarySearch, BinarySearch

LinearBinarySearch{MAX} <: AbstractSearchPolicy

Bounded linear search within a window of MAX positions from hint, then binary fallback. Optimal for sorted/monotonic query sequences where consecutive queries tend to fall in adjacent or nearby intervals.

Type Parameter

• MAX::Int: Size of the linear search window before falling back to binary search. This is a compile-time constant encoded in the type parameter.

Performance Characteristics

• Best case: O(1) when query is within MAX positions of the hint
• Worst case: O(log n) binary search fallback
• Ideal for: Sorted queries, time-series data, streaming evaluations

Construction

Use the factory function (recommended) to construct with a curated set of values:

LinearBinarySearch()                   # default MAX=8
LinearBinarySearch(linear_window=0)    # hint check only, no walk
LinearBinarySearch(linear_window=4)    # custom MAX=4

Or construct the parametric type directly (advanced):

LinearBinarySearch{8}()            # explicit type parameter

Example

sorted_queries = sort(rand(1000))
vals = linear_interp(x, y, sorted_queries; search=LinearBinarySearch(linear_window=8))

See also: BinarySearch

ConstantSeriesInterpolant{Tg, Tv, P, X}

Multi-series constant (step) interpolant with unified matrix storage and SIMD optimization. Shares a single x-grid across N y-series for efficient batch evaluation.

Type Parameters

• Tg: Grid type (Float32 or Float64)
• Tv: Value type (Tg for real, Complex{Tg} for complex)
• P: Search policy type
• X: Grid container type (Vector or Range)

Fields

• x::X: Shared x-grid (Vector or Range)
• y::Matrix{Tv}: Function values (npoints × nseries) series-contiguous
• _transpose::LazyTranspose{Tv}: Lazy point-contiguous layout for scalar SIMD
• extrap::AbstractExtrap: Extrapolation mode
• side::SD: Side selection (NearestSide(), LeftSide(), RightSide())

Memory Layout

Primary storage is series-contiguous (npoints × nseries):

• y[i, k] = value of series k at grid point i
• y[:, k] is contiguous → optimal for vector queries

Lazy transpose (nseries × npoints) for scalar queries:

• y_point[:, i] is contiguous → optimal for SIMD scalar queries

Usage

x = collect(range(0.0, 1.0, 101))
y1, y2, y3 = sin.(2π .* x), cos.(2π .* x), exp.(-x)

sitp = constant_interp(x, [y1, y2, y3])

# Scalar evaluation
vals = sitp(0.5)            # Returns Vector{Float64} of length 3
sitp(output, 0.5)           # In-place

# Vector evaluation
vals = sitp([0.1, 0.5, 0.9])    # Returns Vector of Vectors
sitp([out1, out2, out3], xq)    # In-place (zero allocation)

# Complex values are also supported
y_complex = [exp.(2im * π * x), (1.0+2.0im) .* x]
sitp_complex = constant_interp(x, y_complex)

Implementation Note: mutable struct with const fields

This type uses mutable struct with all const fields (Julia 1.8+) instead of plain struct for performance reasons. See CubicSeriesInterpolant for details.

LinearSeriesInterpolant{Tg, Tv, P, X}

Multi-series linear interpolant with unified matrix storage and SIMD optimization. Shares a single x-grid across N y-series for efficient batch evaluation.

Type Parameters

• Tg: Grid type (Float32 or Float64)
• Tv: Value type (Tg for real, Complex{Tg} for complex)
• P: Search policy type
• X: Grid container type (Vector or Range)

Fields

• x::X: Shared x-grid (Vector or Range)
• y::Matrix{Tv}: Function values (npoints × nseries) series-contiguous
• _transpose::LazyTranspose{Tv}: Lazy point-contiguous layout for scalar SIMD
• extrap::E: Extrapolation mode (compile-time specialized via type parameter)

Memory Layout

Primary storage is series-contiguous (npoints × nseries):

• y[i, k] = value of series k at grid point i
• y[:, k] is contiguous → optimal for vector queries

Lazy transpose (nseries × npoints) for scalar queries:

• y_point[:, i] is contiguous → optimal for SIMD scalar queries

Usage

x = collect(range(0.0, 1.0, 101))
y1, y2, y3 = sin.(2π .* x), cos.(2π .* x), exp.(-x)

sitp = linear_interp(x, [y1, y2, y3])

# Scalar evaluation
vals = sitp(0.5)            # Returns Vector{Float64} of length 3
sitp(output, 0.5)           # In-place

# Vector evaluation
vals = sitp([0.1, 0.5, 0.9])    # Returns Vector of Vectors
sitp([out1, out2, out3], xq)    # In-place (zero allocation)

# Complex values are also supported
y_complex = [exp.(2im * π * x), (1.0+2.0im) .* x]
sitp_complex = linear_interp(x, y_complex)

Performance

• Vector queries use series-contiguous layout directly
• Scalar queries trigger lazy transpose on first call

Implementation Note: mutable struct with const fields

This type uses mutable struct with all const fields (Julia 1.8+) instead of plain struct for performance reasons. See CubicSeriesInterpolant for details.

QuadraticSeriesInterpolant{Tg, Tv, E, P, X}

Multi-series quadratic spline interpolant with unified matrix storage and SIMD optimization. Shares a single x-grid across N y-series for efficient batch evaluation.

Type Parameters

• Tg: Grid type (Float32 or Float64)
• Tv: Value type (Tg for real, Complex{Tg} for complex)
• E: Extrapolation mode type (compile-time specialized)
• P: Search policy type
• X: Grid container type (Vector or Range)

Fields

• x::X: Grid points (sorted, Vector or Range)
• y::Matrix{Tv}: Function values (npoints × nseries) series-contiguous
• a::Matrix{Tv}: Quadratic coefficients (npoints × nseries) series-contiguous
• d::Matrix{Tv}: Slope coefficients (npoints × nseries) series-contiguous
• h::Vector{Tg}: Grid spacing (shared across all series, always real)
• _transpose::LazyTransposeTriple{Tv}: Lazy point-contiguous layout for SIMD
• extrap::E: Extrapolation mode (compile-time specialized via type parameter)

Memory Layout

Primary storage is series-contiguous (npoints × nseries):

• y[i, k] = value of series k at grid point i
• y[:, k] is contiguous → optimal for vector queries

Lazy transpose (nseries × npoints) for scalar queries:

• y_point[:, i] is contiguous → optimal for SIMD scalar queries

Usage

x = collect(range(0.0, 1.0, 101))
y1, y2, y3 = sin.(2π .* x), cos.(2π .* x), exp.(-x)

sitp = quadratic_interp(x, [y1, y2, y3])

# Scalar evaluation
vals = sitp(0.5)            # Returns Vector{Float64} of length 3
sitp(output, 0.5)           # In-place

# Vector evaluation
vals = sitp([0.1, 0.5, 0.9])    # Returns Vector of Vectors
sitp([out1, out2, out3], xq)    # In-place (zero allocation)

# Derivatives
d1 = sitp(0.5; deriv=DerivOp(1))     # First derivatives
d2 = sitp(0.5; deriv=DerivOp(2))     # Second derivatives

# Complex values are also supported
y_complex = [exp.(2im * π * x), (1.0+2.0im) .* x]
sitp_complex = quadratic_interp(x, y_complex)

Performance

• Vector queries use series-contiguous layout directly
• Scalar queries trigger lazy transpose on first call
• All series share same h[] array (O(1) memory overhead)

Implementation Note: mutable struct with const fields

This type uses mutable struct with all const fields (Julia 1.8+) instead of plain struct for performance reasons. See CubicSeriesInterpolant for details.

CubicSeriesInterpolant{Tg, Tv, C, B}

Multi-series cubic spline interpolant with unified matrix storage and SIMD optimization. Shares a single x-grid across N y-series for efficient batch evaluation.

Type Parameters

• Tg: Grid coordinate type (Float32 or Float64) - always real
• Tv: Value type (real or Complex{Tg})
• C: Cache type (CubicSplineCache{Tg})
• B: Boundary condition config type (BCPair or PeriodicData)

Fields

• cache::C: Shared CubicSplineCache with LU factorization
• bc_for_solve::B: BC configuration for solving systems
• y::Matrix{Tv}: Function values (npoints × nseries) series-contiguous
• z::Matrix{Tv}: Second derivatives (npoints × nseries) series-contiguous
• _point_snapshot: Atomic field for lazy point-contiguous layout
• extrap::E: Extrapolation mode (compile-time specialized via type parameter)

Memory Layout

Primary storage is series-contiguous (npoints × nseries):

• y[i, k] = value of series k at grid point i
• y[:, k] is contiguous → optimal for vector queries

Lazy transpose (nseries × npoints) for scalar queries:

• y_point[:, i] is contiguous → optimal for SIMD scalar queries

Usage

x = collect(range(0.0, 1.0, 101))
y1, y2, y3 = sin.(2π .* x), cos.(2π .* x), exp.(-x)

sitp = cubic_interp(x, [y1, y2, y3])

# Scalar evaluation
vals = sitp(0.5)            # Returns Vector{Float64} of length 3
sitp(output, 0.5)           # In-place

# Vector evaluation
vals = sitp([0.1, 0.5, 0.9])    # Returns Vector of Vectors
sitp([out1, out2, out3], xq)    # In-place (zero allocation)

# Complex values
y_complex = exp.(2im .* π .* x)
sitp_c = cubic_interp(x, [y_complex])  # CubicSeriesInterpolant{Float64, ComplexF64, ...}

Performance

• Vector queries use series-contiguous layout directly
• Scalar queries trigger lazy transpose on first call
• All series share same cache (O(1) memory overhead)

Implementation Note: mutable struct with const fields

This type uses mutable struct with all const fields (Julia 1.8+) instead of plain struct for performance reasons. The const annotation ensures fields cannot be reassigned while allowing heap allocation. This pattern provides:

• Stable memory addresses for the compiler to optimize field access
• Better inlining of field accesses compared to plain immutable structs
• Compatibility with mutable inner types (LazyTransposePair with @atomic)

Benchmarks show ~15% regression when using plain struct instead.

precompute_transpose!(lt::LazyTranspose, src::Matrix) -> LazyTranspose

Pre-compute transpose before hot loops. Returns the holder for chaining.

precompute_transpose!(ltp::LazyTransposePair, y::Matrix, z::Matrix) -> LazyTransposePair

Pre-compute transpose pair before hot loops. Returns the holder for chaining.

precompute_transpose!(ltt::LazyTransposeTriple, y::Matrix, a::Matrix, d::Matrix) -> LazyTransposeTriple

Pre-compute transpose triple before hot loops. Returns the holder for chaining.

precompute_transpose!(sitp::LinearSeriesInterpolant) -> sitp

Pre-allocate point-contiguous matrix for scalar queries. Call before hot loops to avoid first-call latency.

precompute_transpose!(sitp::ConstantSeriesInterpolant) -> sitp

Pre-allocate point-contiguous matrix for scalar queries. Call before hot loops to avoid first-call latency.

precompute_transpose!(sitp::CubicSeriesInterpolant) -> sitp

Pre-allocate point-contiguous matrices for scalar queries. Call before hot loops to avoid first-call latency.