Cubic Spline API

cubic_interp(cache, y, x_query; extrap=:none, deriv=0, search=Binary()) -> Vector{T}

Allocating version of cubic spline interpolation using cached LU factorization.

Arguments

• deriv::Int=0: Derivative order (0=value, 1=first derivative, 2=second derivative)
• search::AbstractSearchPolicy=Binary(): Search algorithm for interval finding

Example

cache = CubicSplineCache(collect(range(0.0, 1.0, 51)))
y = sin.(cache.x)
result = cubic_interp(cache, y, [0.25, 0.5, 0.75])
derivs = cubic_interp(cache, y, [0.25, 0.5, 0.75]; deriv=1)  # First derivative

# Optimized for sorted queries
sorted_queries = sort(rand(1000))
vals = cubic_interp(cache, y, sorted_queries; search=LinearBinary(linear_window=8))

cubic_interp(x, y, x_query; bc=NaturalBC(), extrap=:none, autocache=true, deriv=0, search=Binary()) -> Vector{T}

Cubic spline interpolation with optional automatic caching.

Arguments

• deriv::Int=0: Derivative order (0=value, 1=first derivative, 2=second derivative)
• search::AbstractSearchPolicy=Binary(): Search algorithm for interval finding

Extrapolation Modes

• :none (default): Throws DomainError if query point is outside domain
• :constant: Returns boundary values outside domain (0 for derivatives)
• :extension: Extends boundary polynomial outside domain
• :wrap: Wraps coordinates to domain (for sawtooth/triangle patterns)
• For bc=PeriodicBC(): extrapolation is ignored (coordinates are always wrapped)

Example

result = cubic_interp(x, y, x_query)                     # Auto-cached (default)
derivs = cubic_interp(x, y, x_query; deriv=1)            # First derivative
result = cubic_interp(x, y, x_query; extrap=:extension)  # Extend beyond domain

# Optimized for sorted queries
sorted_queries = sort(rand(1000))
vals = cubic_interp(x, y, sorted_queries; search=LinearBinary(linear_window=8))

cubic_interp(cache, y; extrap=:none, search=Binary()) -> CubicInterpolant

Create a callable interpolant from a pre-built cache.

Note on BC Values

Uses cache.bc_config for boundary condition values. This is correct when:

• Cache was built via CubicSplineCache(x; bc=...) with actual BC values
• BC is NaturalBC/ClampedBC/PeriodicBC (values are always zero)

Warning: Caches from get_cubic_cache contain placeholder zeros in bc_config. For non-zero BC values, use the full API: cubic_interp(x, y; bc=Deriv1(val)).

Tg = grid type, Tv = value type (can be Complex)

Thread-Safety

Pool-allocated tmp_z is copied by the CubicInterpolant constructor, so the pool memory can be safely reused after this function returns.

cubic_interp(x, ys::AbstractVector{<:AbstractVector}; bc=NaturalBC(), extrap=:none, autocache=true, precompute_transpose=false)

Create a multi-Y cubic spline interpolant for multiple y-data series sharing the same x-grid.

Arguments

• x::AbstractVector: x-coordinates (sorted, length ≥ 2)
• ys: Vector of y-value vectors (all same length as x)
• bc: Boundary condition (NaturalBC, ClampedBC, PeriodicBC, or Vector of BC for per-series)
• extrap: Extrapolation mode (:none, :constant, :extension, :wrap)
• autocache: If true, reuse cached LU factorization (default: true)
• precompute_transpose: If true, build point-contiguous layout immediately

Returns

CubicSeriesInterpolant object with matrix storage.

Example

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

# Uniform BC for all series
sitp = cubic_interp(x, [y1, y2, y3])
vals = sitp(0.5)  # [sin(π), cos(π), exp(-0.5)]

# Per-series BC (each series can have different BC)
sitp = cubic_interp(x, [y1, y2, y3]; bc=[
    NaturalBC(),
    BCPair(Deriv1(2.0), Deriv1(0.0)),
    BCPair(Deriv2(0.0), Deriv3(5.0)),
])

cubic_interp(x, Y::AbstractMatrix; bc=NaturalBC(), extrap=:none, autocache=true, precompute_transpose=false)

Create a multi-Y cubic spline interpolant from a matrix where each column is a y-series.

Arguments

• x::AbstractVector: x-coordinates (length n)
• Y::AbstractMatrix: n×m matrix, each column is a y-series
• bc, extrap, autocache, precompute_transpose: Same as vector form

Example

x = collect(range(0.0, 1.0, 101))
Y = hcat(sin.(2π .* x), cos.(2π .* x))  # 101×2 matrix

sitp = cubic_interp(x, Y)

cubic_interp(grids::NTuple{N,AbstractVector}, data::AbstractArray{<:Any,N}; kwargs...)

Create an N-dimensional cubic Hermite interpolant from grid vectors and data array.

Arguments

• grids::NTuple{N,AbstractVector}: Tuple of grid vectors for each dimension
• data::AbstractArray{<:Any,N}: Function values at grid points

Keywords

• bc=NaturalBC(): Boundary condition(s). Can be:
  • Single AbstractBC: Applied to all axes
  • NTuple{N,AbstractBC}: Per-axis BCs
• extrap=:none: Extrapolation mode(s). Can be:
  • Single Symbol: Applied to all axes (:none, :constant, :wrap)
  • NTuple{N,Symbol}: Per-axis modes
• search=Binary(): Search policy(s). Can be:
  • Single AbstractSearchPolicy: Applied to all axes
  • NTuple{N,AbstractSearchPolicy}: Per-axis policies
• coeffs=PreCompute(): Coefficient computation strategy

Returns

• CubicInterpolantND{Tg, Tv, N, ...}: Callable interpolant object

Type Inference

• Grid type Tg: Promoted from all grid element types (always AbstractFloat)
• Value type Tv: Element type of data (can be real, complex, or AD types)

Examples

# 3D interpolation
x = range(0.0, 2π, 20)
y = range(0.0, π, 15)
z = range(0.0, 1.0, 10)
data = [sin(xi) * cos(yj) * zk for xi in x, yj in y, zk in z]
itp = cubic_interp((x, y, z), data)
itp((1.0, 0.5, 0.3))  # Evaluate at (1.0, 0.5, 0.3)

# With per-axis options
itp = cubic_interp((x, y, z), data;
    bc=(NaturalBC(), PeriodicBC(), NaturalBC()),
    extrap=(:none, :wrap, :constant))

# Complex-valued data
data_c = [sin(xi) * cos(yj) * zk + im * cos(xi) for xi in x, yj in y, zk in z]
itp_c = cubic_interp((x, y, z), data_c)

cubic_interp(grids, data, query; deriv=0, kwargs...)

One-shot ND cubic interpolation at a single point.

Keywords

• deriv: Int (0-3) or Val((d1,d2,...)) for mixed partials

cubic_interp(grids, data, queries; deriv=0, kwargs...)

One-shot ND cubic interpolation at multiple points (batch).

Keywords

• deriv: Int (0-3) or Val((d1,d2,...)) for mixed partials

cubic_interp!(output, cache, y, x_query; extrap=:none, deriv=0, search=Binary())

In-place cubic spline interpolation using cached LU factorization.

Solves the tridiagonal system ONCE, then evaluates at all query points. Thread-safe: workspaces allocated from task-local pool.

Arguments

• output::AbstractVector{T}: Pre-allocated output buffer (modified in-place)
• cache::CubicSplineCache{T}: Pre-computed cache with LU factorization
• y::AbstractVector{T}: Function values at grid points
• x_query::AbstractVector{T}: Query points
• extrap::Symbol=:none: Extrapolation mode (:none, :constant, :extension, :wrap)
• deriv::Int=0: Derivative order (0=value, 1=first derivative, 2=second derivative)
• search::AbstractSearchPolicy=Binary(): Search algorithm for interval finding

cubic_interp!(output, x, y, x_query; bc=NaturalBC(), extrap=:none, autocache=true, deriv=0, search=Binary())

In-place cubic spline interpolation with optional automatic caching.

CubicInterpolant{Tg,Tv,C,P,BC}

Lightweight callable interpolant for broadcast fusion optimization. Returned by cubic_interp(x, y) (2-argument form).

Type Parameters

• Tg<:AbstractFloat: Grid type (Float32 or Float64) for x-coordinates
• Tv: Value type for y-values (can be Tg, Complex{Tg}, or other Number)
• C: CubicSplineCache type (preserves grid type info for O(1) vs O(log n) lookup)
• P: Search policy type (Binary, HintedBinary, LinearBinary, etc.)
• BC: Boundary condition type (BCPair or PeriodicBC)

Fields

• cache::C: Pre-computed CubicSplineCache (LU factorization)
• y::Vector{Tv}: y-values (function values at grid points)
• z::Vector{Tv}: Pre-computed second derivative coefficients (solves system once!)
• bc::BC: Boundary condition used for this interpolant
• extrap::ExtrapVal: Extrapolation mode (union-split for efficient dispatch)
• search_policy::P: Default search policy for interval lookup

Usage

itp = cubic_interp(x, y)
result = @. coef * itp(rho) * other_terms  # fused, zero-allocation per call
val = itp(0.5)                              # scalar (zero-allocation)

# Create with custom search policy
itp = cubic_interp(x, y; search=LinearBinary())
val = itp(0.5)                              # uses LinearBinary() by default
val = itp(0.5; search=Binary())             # override with Binary()

# Complex values
x = [0.0, 1.0, 2.0, 3.0, 4.0]
y = [1.0+2.0im, 3.0+4.0im, 5.0+6.0im, 7.0+8.0im, 9.0+10.0im]
itp = cubic_interp(x, y)
val = itp(0.5)  # returns ComplexF64

Performance Notes

• System solved ONCE at construction -> z coefficients pre-computed
• Each scalar call just evaluates cubic polynomial (zero-allocation!)
• Broadcast operations are perfectly fused (no intermediate arrays)
• Extrapolation mode uses union-splitting for near-zero overhead dispatch

CubicSplineCache{T,X,F,BC,S}

Cache structure for cubic spline interpolation with reusable Thomas factorization.

Type Parameters

• T: Float type (Float32 or Float64)
• X: Grid type (Vector{T} or AbstractRange{T})
• F: Factorization type (ThomasFactorization{T,V})
• BC: Boundary condition data type (BCPair{L,R} for derivative BC, PeriodicData{T} for periodic)
• S: Grid spacing type (ScalarSpacing{T} for Range, VectorSpacing{T} for Vector)

Fields

• x::X: Grid points (immutable after construction, can be Range or Vector)
• spacing::S: Grid spacing data (ScalarSpacing for uniform, VectorSpacing for non-uniform)
• thomas::F: Thomas factorization of tridiagonal matrix A (contains dl, du, inv_d)
• bc_config::BC: Boundary condition data (BCPair for derivative BC, PeriodicData for periodic)

Notes

The Thomas factorization depends ONLY on x geometry and can be reused for:

• Different y vectors (varying function values)
• Different x_query vectors (varying query points)

When x is an AbstractRange, O(1) index lookup is used instead of O(log n) binary search.

Thread-Safety

Workspaces (d, z) are allocated from task-local pools via @with_pool, not stored in this struct. This makes the cache thread-safe by design.

Performance

The spacing field stores grid spacing with precomputed reciprocals.

• ScalarSpacing: O(1) memory for uniform grids (Range inputs) with constant propagation
• VectorSpacing: O(N) memory for non-uniform grids (Vector inputs)

Both store inv_h for eliminating floating-point division in kernel hot paths. On ARM64 (Apple Silicon), this reduces per-evaluation latency from ~10 cycles (fdiv) to ~4 cycles (fmul) — a 2.5× speedup in the inner loop.

Boundary Conditions

• bc=NaturalBC() (default): Natural spline with z[1] = z[n+1] = 0
• bc=PeriodicBC(): Periodic spline with C2 continuity at boundaries

set_cubic_cache_size!(n::Int)

Set maximum cache size for future Julia sessions via Preferences.jl.

Requires Julia restart to take effect (cache size is a load-time constant for thread-safety and compiler optimization).

Arguments

• n::Int: Maximum cache entries (default: 16)

Example

set_cubic_cache_size!(32)  # Sets preference
# Restart Julia to apply new size
get_cubic_cache_size()     # Now returns 32

get_cubic_cache_size()

Get current maximum cache size (load-time constant).

clear_cubic_cache!()

Clear all cached x-grids.

Thread Safety (RCU)

Atomically replaces registry snapshots with empty vectors. Existing readers continue to see their old snapshots until they finish.