Polyharmonic Spline API

FastInterpolations.Experimental

Features whose API and behaviour may change, or be removed, in any release. Names here are not covered by semantic versioning; bring them in explicitly with using FastInterpolations.Experimental.

Names exported here:

• phs_interp, phs_interp!, PHSInterpolantND — polyharmonic spline interpolation on uniform rectilinear grids.

Experimental features reached through existing names (nothing to import from here):

• GPU arrays: the existing API on device arrays, active once a GPU package is loaded (Adapt / GPUArrays extensions). See the GPU Arrays guide for what may change.

phs_interp(grids, data; kwargs...) -> PHSInterpolantND

Create an N-dimensional polyharmonic spline interpolant.

Experimental: load with using FastInterpolations.Experimental; the API and behaviour may change in any release. Inputs it cannot evaluate correctly raise an ArgumentError.

Arguments

• grids: NTuple{N, AbstractVector} — one grid vector per dimension. Each axis must be uniformly spaced: ranges are accepted as is, other vectors must match their mean spacing up to rounding.
• data: AbstractArray{Tv, N} — real data values at grid nodes

Keyword Arguments

• stencil_size::Int = 8: Number of stencil nodes per axis (total = stencil_size^N). Reduce for high dimensions (e.g. 4 for N≥4). The stencil must not span more nodes than an axis has.
• degree::Int = 3: PHS radial function degree (odd positive integer: 1, 3, 5, …). Higher degree → smoother interpolant, larger condition number.
• blend_factor::Real = 1.0: Blend range = blendfactor × maxgrid_spacing. Larger values → wider blending neighbourhood → smoother but more expensive. Default 1.0 provides good balance (3× faster than 2.0, ~2× error increase). Values that leave cell centres outside every blend radius are rejected (in 4-D and above this includes the default; pass e.g. 2.0).
• extrap=NoExtrap(): Extrapolation mode (scalar or per-axis tuple): NoExtrap, FillExtrap or InBounds.
• search=AutoSearch(): Search policy (scalar or per-axis tuple; used for OOB checking).
• log_reference=nothing: Enables the log transform: data is stored as log(ρ/ρ₀) and evaluation returns ρ₀ · exp(f) with matching derivatives. ρ₀ is a nonzero constant, used as given and with the sign of the data, or a callable ρ₀(q) on an NTuple{N} point. Value queries only call ρ₀(q); derivative queries also call ρ₀(q; deriv = ops), as for a FastInterpolations interpolant, so a callable reference answers them through the deriv keyword. Every data/ρ₀ must be positive and finite. Arrays are not accepted: a reference sampled on the data grid adds no information.

Returns

PHSInterpolantND{Tg, Tv, N, degree} — callable interpolant.

Examples

x = range(0.0, 1.0, 20)
y = range(0.0, 1.0, 20)
data = [sin(xi) * cos(yj) for xi in x, yj in y]

itp = phs_interp((x, y), data)
itp((0.5, 0.3))                              # scalar query
itp((0.5, 0.3); deriv=DerivOp(1, 0))        # ∂f/∂x
itp(([0.1, 0.5, 0.9], [0.2, 0.4, 0.6]))     # batch SoA

phs_interp(grids, data, query::NTuple{N,Real}; kwargs...) -> scalar

One-shot N-dimensional PHS interpolation at a single query point.

See phs_interp(grids, data) for keyword argument documentation.

phs_interp(grids, data, queries; kwargs...) -> Vector

One-shot N-dimensional PHS interpolation at a batch of query points. queries is any query-protocol-compatible container (SoA tuple, AoS vector, etc.).

Builds a temporary interpolant (same construction cost as phs_interp(grids, data)), then allocates and fills the output vector.

phs_interp!(out, grids, data, queries; kwargs...)

In-place one-shot N-dimensional PHS interpolation. Writes results into pre-allocated out.

PHSInterpolantND{Tg, Tv, N, K}

N-dimensional local polyharmonic spline interpolant.

Implements the method from the paper, combining:

1. Local stencil-based PHS interpolation (φ(r) = r^K, K odd)
2. Weighted blending across neighbouring base-node interpolants for C² continuity
3. Optional log-density smoothing transform

A single canonical stencil geometry (and its Φ⁻¹) is precomputed once from the mean grid spacings. At boundary nodes the same Φ⁻¹ is reused with clamped data indices — identical to the reference Fortran implementation.

Type Parameters

• Tg: Grid float type
• Tv: Value float type (real floating point only)
• N: Number of dimensions
• K: PHS degree (1, 3, 5, …)

Fields

• grids: Per-axis grid vectors
• data: N-D data array (or log(ρ/ρ₀) when transform is active)
• stencil_offsets: Single canonical stencil: stencil_size^N integer offsets from origin
• phi_inv: Single Φ⁻¹ matrix for the canonical stencil
• hs: Per-axis mean grid spacing used to build stencil_offsets/phi_inv
• blend_a: Blending range parameter (≥ max grid spacing × blend_factor)
• blend_r_idx: Per-axis half-width of blend neighbourhood in index space
• transform: nothing, or the log-transform state built from log_reference
• extraps: Per-axis extrapolation modes
• searches: Per-axis search policies (used for OOB checking only)

Performance

• Construction: O(M³) for one Φ⁻¹ (M = stencil_size^N + N + 1)
• Query: O(nblend × Nstencil × M) where nblend = number of neighbours within blenda
• Memory: O(M²) for Φ⁻¹ plus O(prod(grid_sizes)) for data

Thread-Safety

Safe to call concurrently from externally-threaded loops: the mutable per-thread coefficient cache is indexed by Threads.threadid() so each thread uses its own slot.