Cubic Spline API

cubic_interp(cache, y, x_query; extrap=:none, deriv=0) -> 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)

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

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

Cubic spline interpolation with optional automatic caching.

Arguments

• deriv::Int=0: Derivative order (0=value, 1=first derivative, 2=second derivative)

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

cubic_interp(cache, y; extrap=:none) -> 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)).

Thread-Safety

Uses task-local pool for workspace allocation. Pool memory is copied by the CubicInterpolant constructor.

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!(output, cache, y, x_query; extrap=:none, deriv=0)

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)

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

In-place cubic spline interpolation with optional automatic caching.

CubicInterpolant{T,C}

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

Type Parameters

• T: Float type (Float32 or Float64)
• C: CubicSplineCache type (preserves grid type info for O(1) vs O(log n) lookup)

Fields

• cache::C: Pre-computed CubicSplineCache (LU factorization)
• y::Vector{T}: y-values (function values at grid points)
• z::Vector{T}: Pre-computed second derivative coefficients (solves system once!)
• extrap::ExtrapVal: Extrapolation mode (union-split for efficient dispatch)

Usage

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

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 LU factorization.

Type Parameters

• T: Float type (Float32 or Float64)
• X: Grid type (Vector{T} or AbstractRange{T})
• F: LU factorization type
• BC: Boundary condition data type (BCPair{T,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)
• lu_factor::F: LU factorization of tridiagonal matrix A
• bc_config::BC: Boundary condition data (BCPair for derivative BC, PeriodicData for periodic)

Notes

The LU 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.

deriv1(itp::AbstractInterpolant)
deriv2(itp::AbstractInterpolant)
deriv3(itp::AbstractInterpolant)

Create a callable, zero-allocation derivative view of the interpolant for the 1st, 2nd, or 3rd derivative.

DerivativeView is a lightweight wrapper that delegates all evaluation calls to the underlying interpolant using the deriv keyword argument (e.g., itp(xq; deriv=1)). This enables a more functional syntax without the overhead of full object copying.

Features

• Functional API: Pass derivative views into high-order functions like map or find_zeros.
• Broadcast Fusion: Enables evaluation without temporary allocations when combined with other operations using dot syntax (e.g., @. d1(xs) + d2(xs)).
• Unified Interface: Supports both single-series and multi-series (SeriesInterpolant) objects.
  • For AbstractInterpolant, returns a scalar derivative.
  • For AbstractSeriesInterpolant, returns a vector of derivatives.

Notes

• Order 3: For cubic splines, the 3rd derivative is piecewise constant and may jump at knot points. For lower-order interpolants, it returns zero.
• In-place: Supports in-place evaluation for SeriesInterpolant via d(out, xq).

Example

itp = cubic_interp(x, y)
d1 = deriv1(itp)

# Evaluate at a point
val = d1(0.5)

# Use in higher-order functions
using Roots
root = find_zero(d1, 0.5)

# Fused broadcast (no temporary allocations)
query_pts = range(0.0, 1.0, 100)
results = @. itp(query_pts) + 0.1 * d1(query_pts)

deriv1(itp::AbstractInterpolant)
deriv2(itp::AbstractInterpolant)
deriv3(itp::AbstractInterpolant)

Create a callable, zero-allocation derivative view of the interpolant for the 1st, 2nd, or 3rd derivative.

DerivativeView is a lightweight wrapper that delegates all evaluation calls to the underlying interpolant using the deriv keyword argument (e.g., itp(xq; deriv=1)). This enables a more functional syntax without the overhead of full object copying.

Features

• Functional API: Pass derivative views into high-order functions like map or find_zeros.
• Broadcast Fusion: Enables evaluation without temporary allocations when combined with other operations using dot syntax (e.g., @. d1(xs) + d2(xs)).
• Unified Interface: Supports both single-series and multi-series (SeriesInterpolant) objects.
  • For AbstractInterpolant, returns a scalar derivative.
  • For AbstractSeriesInterpolant, returns a vector of derivatives.

Notes

• Order 3: For cubic splines, the 3rd derivative is piecewise constant and may jump at knot points. For lower-order interpolants, it returns zero.
• In-place: Supports in-place evaluation for SeriesInterpolant via d(out, xq).

Example

itp = cubic_interp(x, y)
d1 = deriv1(itp)

# Evaluate at a point
val = d1(0.5)

# Use in higher-order functions
using Roots
root = find_zero(d1, 0.5)

# Fused broadcast (no temporary allocations)
query_pts = range(0.0, 1.0, 100)
results = @. itp(query_pts) + 0.1 * d1(query_pts)

deriv1(itp::AbstractInterpolant)
deriv2(itp::AbstractInterpolant)
deriv3(itp::AbstractInterpolant)

Create a callable, zero-allocation derivative view of the interpolant for the 1st, 2nd, or 3rd derivative.

DerivativeView is a lightweight wrapper that delegates all evaluation calls to the underlying interpolant using the deriv keyword argument (e.g., itp(xq; deriv=1)). This enables a more functional syntax without the overhead of full object copying.

Features

• Functional API: Pass derivative views into high-order functions like map or find_zeros.
• Broadcast Fusion: Enables evaluation without temporary allocations when combined with other operations using dot syntax (e.g., @. d1(xs) + d2(xs)).
• Unified Interface: Supports both single-series and multi-series (SeriesInterpolant) objects.
  • For AbstractInterpolant, returns a scalar derivative.
  • For AbstractSeriesInterpolant, returns a vector of derivatives.

Notes

• Order 3: For cubic splines, the 3rd derivative is piecewise constant and may jump at knot points. For lower-order interpolants, it returns zero.
• In-place: Supports in-place evaluation for SeriesInterpolant via d(out, xq).

Example

itp = cubic_interp(x, y)
d1 = deriv1(itp)

# Evaluate at a point
val = d1(0.5)

# Use in higher-order functions
using Roots
root = find_zero(d1, 0.5)

# Fused broadcast (no temporary allocations)
query_pts = range(0.0, 1.0, 100)
results = @. itp(query_pts) + 0.1 * d1(query_pts)