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 (unconstrained) - used for y-values, coefficients, return values. Any type supporting the 5 core operations (+, -, TgTv, TvTg, Int*Tv).

Design Invariant

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

Subtypes

• AbstractInterpolant1D{Tg, Tv}: 1D interpolants (shared callable protocol)
• AbstractSeriesInterpolant{Tg, Tv}: Multi-series interpolants
• AbstractInterpolantND{Tg, Tv, N}: N-dimensional interpolants

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 (unconstrained)

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 (unconstrained)
• 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)

AbstractAdjointND{Tg<:AbstractFloat, N} <: AbstractAdjoint{Tg}

Abstract supertype for N-dimensional adjoint operators.

Subtypes automatically inherit shared callable dispatch from nd_adjoint_protocol.jl (allocating, in-place, scalar, tuple) by implementing the required interface:

_n_queries(adj)::Int                          — number of baked query points
_grid_size(adj)::NTuple{N,Int}                — internal grid size
_adjoint_bcs(adj)                             — boundary conditions tuple
_adjoint_nd_apply!(f_bar, adj, y_bar, ops)    — scatter y_bar into f_bar (accumulate)

Subtypes

• CubicAdjointND{Tg, N, ...}: Adjoint of cubic spline interpolation (ND)
• LinearAdjointND{Tg, N, ...}: Adjoint of linear interpolation (ND)

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)
• AutoCoeffs: Automatic selection based on context

PreCompute <: AbstractCoeffStrategy

Precompute all partial derivatives at construction time.

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

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)

OnTheFly <: AbstractCoeffStrategy

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

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)

AutoCoeffs <: AbstractCoeffStrategy

Automatic coefficient strategy selection based on context.

Rules

• 1D scalar query → OnTheFly() (O(1) local slopes beats O(n) bulk)
• 1D vector query → length(xq) > length(x) ? PreCompute() : OnTheFly()
• 1D interpolant → PreCompute() (amortize for reuse)
• ND trivial methods (Linear/Constant) → PreCompute() (no slopes)
• ND N ≥ 3 → OnTheFly() (2^N memory savings)
• ND any local Hermite method → OnTheFly()

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 (unconstrained)
• 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)

HeteroInterpolantND{Tg, Tv, N, G, S, M, E, P, D} <: AbstractInterpolantND{Tg, Tv, N}

N-dimensional interpolant with per-axis method specification.

Supports heterogeneous methods across dimensions (e.g., cubic on axis 1, linear on axis 2) with two evaluation strategies:

• OnTheFly(): Sequential 1D interpolation per query (O(n) per query, zero build cost)
• PreCompute(): Precomputed partial derivatives with local kernel eval (O(1) per query)

Type Parameters

• Tg: Grid/coordinate type (Float32 or Float64)
• Tv: Value type (unconstrained)
• N: Number of dimensions
• G: Tuple type for grids
• S: Tuple type for spacings
• M: Tuple type for per-axis methods (heterogeneous)
• E: Tuple type for extrapolation modes
• P: Tuple type for search policies
• D: Data storage type — Array{Tv,N} (OnTheFly) or _HeteroPartials{Tv,N,N+1} (PreCompute)

Example

x, y = range(0, 1, 50), range(0, 1, 30)
data = [sin(xi) * cos(yj) for xi in x, yj in y]

# On-the-fly (default)
itp = interp((x, y), data; method=(CubicInterp(), LinearInterp()))

# Precomputed (fast eval)
itp = interp((x, y), data; method=(CubicInterp(), LinearInterp()), coeffs=PreCompute())
itp((0.5, 0.3))

AbstractInterpMethod

Abstract supertype for per-axis interpolation method specification. Used by HeteroInterpolantND to specify interpolation method per dimension.

Subtypes

• CubicInterp: Cubic spline interpolation (C² continuous)
• LinearInterp: Linear interpolation (C⁰ continuous)
• QuadraticInterp: Quadratic spline interpolation (C¹ continuous)
• ConstantInterp: Constant (nearest-neighbor) interpolation

Example

methods = (CubicInterp(), LinearInterp())  # cubic on axis 1, linear on axis 2
itp = interp((x, y), data; method=methods)

PchipInterp <: AbstractInterpMethod

PCHIP (monotone-preserving) interpolation method for one axis. Slopes computed via Fritsch-Carlson algorithm. Requires ≥2 grid points.

CardinalInterp{T} <: AbstractInterpMethod

Cardinal spline interpolation method for one axis. tension=0 gives Catmull-Rom. Requires ≥2 grid points.

Arguments

• tension: Tension parameter (0 = CatmullRom, 1 = zero slopes)

AkimaInterp <: AbstractInterpMethod

Akima (1970) outlier-robust interpolation method for one axis. Requires ≥2 grid points.

CubicHermiteInterp <: AbstractInterpMethod

Cubic Hermite interpolation with user-supplied slopes. Type-only in Phase 1 — ND support requires separate slope design.

CubicInterp(; bc::AbstractBC = CubicFit())

Cubic spline interpolation method for one axis. Requires ≥4 grid points (≥3 for PeriodicBC).

Arguments

• bc: Boundary condition — CubicFit(), ZeroCurvBC(), PeriodicBC(), etc.

LinearInterp()

Linear interpolation method for one axis. Requires ≥2 grid points.

QuadraticInterp(; bc::AbstractBC = Left(QuadraticFit()))

Quadratic spline interpolation method for one axis. Requires ≥3 grid points.

Arguments

• bc: Boundary condition — Left(QuadraticFit()), Right(QuadraticFit()), MinCurvFit, etc. Validated by the 1D quadratic_interp constructor.

ConstantInterp(; side::AbstractSide = NearestSide())

Constant (nearest-neighbor) interpolation method for one axis. Requires ≥2 grid points.

Arguments

• side: Side selection — NearestSide(), LeftSide(), RightSide().

NoInterp()

No-interpolation method for one axis. Declares that this axis is discrete and will be queried via GridIdx at evaluation time.

At build time, no differentiation or precomputation is performed for NoInterp axes. At eval time, the axis is sliced at the given GridIdx index — no search, no kernel.

Examples

# 3D data: interpolate x,y; select z by index
itp = interp((x, y, z), data; method=(CubicInterp(), LinearInterp(), NoInterp()))
itp((0.5, 0.3, GridIdx(10)))    # cubic×linear on the z=10 slice
itp((0.5, 0.3, GridIdx(1)))     # same interpolant, different slice

GridIdx(k::Integer)

Grid-index query coordinate. Wraps an integer index for direct grid-point lookup.

After resolution (internal), carries both the index and the grid coordinate value. Search functions short-circuit when they see a GridIdx — zero search cost.

GridIdx <: Real: it flows through Tuple{Vararg{Real, N}} dispatch transparently. Before resolution, val = NaN — a poison sentinel that propagates visibly if _resolve_grididx is ever skipped. After resolution, val holds the grid coordinate and arithmetic auto-promotes via promote_rule (stripping the GridIdx wrapper).

Examples

# One-shot: slice axis 2 at index 5, interpolate axis 1
interp((x, y), data, (0.5, GridIdx(5)); method=(CubicInterp(), NoInterp()))

# Interpolant: query-time slicing
itp = interp((x, y), data; method=(CubicInterp(), NoInterp()))
itp((0.5, GridIdx(5)))

# GridIdx works on ANY axis (not just NoInterp):
itp_hetero = interp((x, y), data; method=(CubicInterp(), LinearInterp()))
itp_hetero((0.5, GridIdx(10)))   # search short-circuited on axis 2

interp(grids, data; method, coeffs=PreCompute(), extrap=NoExtrap(), search=AutoSearch())

Unified N-dimensional interpolation constructor with per-axis method specification.

Automatically dispatches to the optimal implementation:

• Homogeneous (all axes same method): delegates to existing optimized types (CubicInterpolantND, LinearInterpolantND, etc.) with full feature support (adjoint, oneshot, AD)
• Heterogeneous (mixed methods): creates HeteroInterpolantND

Arguments

• grids::NTuple{N, AbstractVector}: Grid vectors per dimension
• data::AbstractArray{<:Any, N}: N-dimensional data array

Keyword Arguments

• method: Interpolation method(s) (required)
  • Single AbstractInterpMethod: broadcast to all axes (e.g., method=CubicInterp())
  • Tuple{Vararg{AbstractInterpMethod, N}}: per-axis (e.g., method=(CubicInterp(), LinearInterp()))
• coeffs=AutoCoeffs(): Coefficient strategy (default: automatic selection)
  • AutoCoeffs(): Automatic — picks best strategy based on methods and dimensionality
  • PreCompute(): Precompute partial derivatives (O(1) eval)
  • OnTheFly(): Build 1D per query (zero build cost, O(n) eval)
• extrap=NoExtrap(): Extrapolation mode(s) — single or per-axis tuple
• search=AutoSearch(): Search policy(ies) — single or per-axis tuple

Examples

x, y = range(0, 1, 50), range(0, 1, 30)
data = [sin(xi) * cos(yj) for xi in x, yj in y]

# Single method → all axes (dispatches to CubicInterpolantND)
itp = interp((x, y), data; method=CubicInterp())

# Heterogeneous → HeteroInterpolantND
itp = interp((x, y), data; method=(CubicInterp(), LinearInterp()))

# Per-axis BCs
itp = interp((x, y), data; method=(CubicInterp(CubicFit()), CubicInterp(ZeroCurvBC())))

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

interp(grids, data, query; method, coeffs=AutoCoeffs(), deriv=EvalValue(), extrap=NoExtrap(), search=AutoSearch(), hint=nothing)

One-shot N-dimensional interpolation at a single point. Zero-allocation after warmup. No interpolant object is created.

With coeffs=AutoCoeffs() (default), scalar queries use OnTheFly() strategy (2^N× less work than PreCompute() for a single point).

Examples

x, y = range(0, 1, 50), range(0, 1, 30)
data = [sin(xi) * cos(yj) for xi in x, yj in y]

# One-shot (no interpolant created)
val = interp((x, y), data, (0.5, 0.3); method=(CubicInterp(), LinearInterp()))

# With derivative
dfdx = interp((x, y), data, (0.5, 0.3);
    method=(CubicInterp(), LinearInterp()), deriv=(DerivOp(1), DerivOp(0)))

interp(grids, data, queries; method, coeffs=AutoCoeffs(), kwargs...)

Allocating one-shot N-dimensional interpolation at multiple points. Returns a Vector of interpolated values.

interp!(output, grids, data, queries; method, coeffs=AutoCoeffs(), kwargs...)

In-place one-shot N-dimensional interpolation at multiple points. Builds partials once, evaluates at all query points.

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

Get the grid/coordinate type of an interpolant.

grid_type(::AbstractAdjoint{Tg}) -> Type{Tg}

Get the grid/coordinate type of an adjoint operator.

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 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} for N ≥ 4: Returns zero (beyond polynomial degree of supported interpolants)

DerivOp{N} <: AbstractEvalOp

Parametric singleton for derivative order dispatch. N is the derivative order (0 = value, 1 = first derivative, etc.). Any non-negative integer is valid; orders beyond the interpolant's polynomial degree (e.g., DerivOp{4} for cubic) return zero automatically.

Construction

• DerivOp{0}(), DerivOp{1}() — direct parametric construction
• DerivOp(n::Int) — convenience: DerivOp(1) → DerivOp{1}() (any n ≥ 0)
• 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
deriv_order(::Type{DerivOp{N}}) -> Int

Extract derivative order from a DerivOp instance or type. The type-level method is used inside @generated functions where only types are available.

DerivOp{N} <: AbstractEvalOp

Parametric singleton for derivative order dispatch. N is the derivative order (0 = value, 1 = first derivative, etc.). Any non-negative integer is valid; orders beyond the interpolant's polynomial degree (e.g., DerivOp{4} for cubic) return zero automatically.

Construction

• DerivOp{0}(), DerivOp{1}() — direct parametric construction
• DerivOp(n::Int) — convenience: DerivOp(1) → DerivOp{1}() (any n ≥ 0)
• 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.). Any non-negative integer is valid; orders beyond the interpolant's polynomial degree (e.g., DerivOp{4} for cubic) return zero automatically.

Construction

• DerivOp{0}(), DerivOp{1}() — direct parametric construction
• DerivOp(n::Int) — convenience: DerivOp(1) → DerivOp{1}() (any n ≥ 0)
• 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.). Any non-negative integer is valid; orders beyond the interpolant's polynomial degree (e.g., DerivOp{4} for cubic) return zero automatically.

Construction

• DerivOp{0}(), DerivOp{1}() — direct parametric construction
• DerivOp(n::Int) — convenience: DerivOp(1) → DerivOp{1}() (any n ≥ 0)
• 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.). Any non-negative integer is valid; orders beyond the interpolant's polynomial degree (e.g., DerivOp{4} for cubic) return zero automatically.

Construction

• DerivOp{0}(), DerivOp{1}() — direct parametric construction
• DerivOp(n::Int) — convenience: DerivOp(1) → DerivOp{1}() (any n ≥ 0)
• 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 (unconstrained)
• 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

nodal_partials(itp, order::NTuple{N, Integer}) -> AbstractArray{Tv, N}

Return a zero-copy view of precomputed partial derivatives at grid nodes.

The order tuple specifies which axes are differentiated: 0 means no differentiation, 1 means first derivative along that axis.

Examples

itp = interp((x, y), data; method=CubicInterp())

nodal_partials(itp, (0, 0))  # f(xᵢ, yⱼ)      — function values
nodal_partials(itp, (1, 0))  # ∂f/∂x(xᵢ, yⱼ)  — x-derivative at nodes
nodal_partials(itp, (0, 1))  # ∂f/∂y(xᵢ, yⱼ)  — y-derivative at nodes
nodal_partials(itp, (1, 1))  # ∂²f/∂x∂y        — mixed partial at nodes

For heterogeneous interpolants, only axes using CubicInterp or QuadraticInterp support order[d] = 1. Requesting a derivative on a LinearInterp/ConstantInterp/NoInterp axis throws ArgumentError.

Supported types

• CubicInterpolantND — all 2^N derivative combinations
• QuadraticInterpolantND — all 2^N combinations (mixed partials are zero)
• HeteroInterpolantND with PreCompute() — per-axis validation

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!, value_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, value_gradient, hessian!

value_gradient(itp::AbstractInterpolantND, query)

Compute the value and gradient simultaneously at query.

Returns (value, gradient) where gradient is an NTuple{N} of partial derivatives.

Performance

Faster than calling itp(query) and gradient(itp, query) separately: interval search is performed only once per query point.

Examples

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

# Optim.jl fg! pattern:
function fg!(F, G, x)
    if G !== nothing && F !== nothing
        val, grad = value_gradient(itp, Tuple(x))
        G .= grad
        return val
    elseif G !== nothing
        gradient!(G, itp, Tuple(x))
        return nothing
    else
        return itp(Tuple(x))
    end
end
result = optimize(Optim.only_fg!(fg!), x0, LBFGS())

See also: gradient, gradient!

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 (unconstrained — supports any type with +, -, scalar *).

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)
Deriv1(MyStruct(0.0))  # Custom type (duck typing)

Deriv2{Tv} <: PointBC

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

The type parameter Tv is the value type (unconstrained — supports any type with +, -, scalar *).

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 (unconstrained — supports any type with +, -, scalar *).

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, C} <: 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)
• C::Bool: whether to validate y[1] ≈ y[end] at construction time (default true)

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))

# Skip endpoint check (e.g., scaled data where noise > 8eps)
itp = cubic_interp(x, 1e6 .* sin.(x); bc=PeriodicBC(check=false))

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
• ClampExtrap: Clamp to nearest boundary value
• FillExtrap: Return a user-specified constant for out-of-domain queries
• ExtendExtrap: Extend interpolation polynomial beyond domain
• WrapExtrap: Wrap queries into domain (periodic)
• InBounds: Caller guarantees in-domain; skip domain checks (advanced/internal)

See also: ConstExtrap() (deprecated factory, forwards to ClampExtrap()).

NoExtrap <: AbstractExtrap

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

Example

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

ClampExtrap <: AbstractExtrap

Constant extrapolation — clamps to nearest boundary value for out-of-domain queries.

Returns y[1] for queries below domain, y[end] for queries above domain. In ND, derivatives along OOB axes are zero; orthogonal derivatives are computed at the clamped boundary point.

Example

itp = cubic_interp(x, y; extrap=ClampExtrap())  # clamp to boundary values
itp = cubic_interp(x, y; extrap=ClampExtrap())    # equivalent

FillExtrap{T} <: AbstractExtrap

Fill extrapolation — returns a user-specified constant value for out-of-domain queries. All derivatives are zero when out-of-domain (constant function).

Standard numerics auto-promote to float; duck types (SVector, etc.) are stored as-is. Follows the same pattern as Deriv1{Tv}.

Examples

itp = cubic_interp(x, y; extrap=FillExtrap(NaN))      # NaN outside domain
itp = cubic_interp(x, y; extrap=FillExtrap(0.0))       # zero outside domain
itp = cubic_interp(x, y; extrap=FillExtrap(; fill_value=0))  # kwarg form

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())

InBounds <: AbstractExtrap

Caller guarantees all queries are within the interpolation domain. Skips domain validation for maximum performance.

Used internally by vector loops after batch _check_domain validation (NoExtrap → InBounds conversion), and available to advanced users who have pre-validated their query points.

Example

# Skip domain check when you know queries are in-domain
linear_interp(x, y, xq; extrap=InBounds())

ConstExtrap()

Deprecated factory function. Use ClampExtrap() directly instead.

Examples

ConstExtrap()  # → ClampExtrap() (with deprecation warning)

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

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

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

Series(y1, y2, ...)                     # varargs
Series([y1, y2, ...])                   # vector of vectors
Series(Y::AbstractMatrix)               # matrix (columns = series)

Input wrapper for multi-series interpolation. Wraps series data to disambiguate from vector-valued interpolation inputs.

The wrapped data is not copied or promoted — type promotion is handled by the interpolation constructor.

Examples

# One-shot evaluation (search once, kernel per y)
linear_interp(x, Series(y_sin, y_cos), 0.5)  # → [sin(0.5), cos(0.5)]

# Interpolant construction
itp = linear_interp(x, Series(y_sin, y_cos))
itp(0.5)  # → [sin(0.5), cos(0.5)]

Number of series in a Series wrapper.

Number of series in the interpolant.

Number of series in the interpolant.

Number of series in the interpolant.

Number of series in the interpolant.

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 (unconstrained)
• 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 (unconstrained)
• 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 (unconstrained)
• 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
• spacing::S: Precomputed grid spacing (ScalarSpacing for Range, VectorSpacing for Vector)
• _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 grid spacing (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 (unconstrained)
• 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.