Polyharmonic Spline API

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

Create an N-dimensional polyharmonic spline interpolant.

Arguments

• grids: NTuple{N, AbstractVector} — one grid vector per dimension. Each axis must be uniformly spaced: the stencil geometry uses a single spacing per axis, so non-uniform axes are accepted but give inaccurate results.
• data: AbstractArray{Tv, N} — 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).
• 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).
• extrap=NoExtrap(): Extrapolation mode (scalar or per-axis tuple).
• search=AutoSearch(): Search policy (scalar or per-axis tuple; used for OOB checking).
• reference_interp=nothing: If provided, enables the log-density smoothing transform: data is stored as log(ρ/ρ₀) where ρ₀ values come from reference_data (if given) or from evaluating reference_interp at each grid node. Evaluation returns ρ̃ = ρ₀ * exp(f) with correct derivative transforms.
• reference_data=nothing: Pre-computed ρ₀ values at all grid nodes (same shape as data). When provided alongside reference_interp, avoids evaluating reference_interp at every grid node during construction — useful when reference_interp is a nested PHS (expensive per-node) but the raw ρ₀ array is already available. Example: ```juliaitp_rho0 is a log-space PHS — accurate derivatives but slow to query 592K nodesitprho0 = phsinterp(grids, rho0; reference_interp = ConstantRef(1.0))Pass rho0 directly so construction stays O(N·log N); eval uses itp_rho0 for ∂ρ₀itpphs = phsinterp(grids, rho; referenceinterp = itprho0, reference_data = rho0) ```

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 type (Float, Complex, duck-typed scalar)
• 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 PHSLogTransform for log-density mode
• 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.

PHSLogTransform{N, Tr}

Optional log-density smoothing transformation container. When active, data stored in PHSInterpolantND contains log(ρ/ρ₀) and evaluation applies the inverse transform (Eqs. 21–23 from the paper) to recover the interpolated density and its derivatives.

Tr can be any callable supporting ref(query) (value) and ref(query; deriv=ops) (derivative), including:

• Any AbstractInterpolantND (cubic spline, PHS, etc.)
• A log-space PHS built with ConstantRef(1.0) for accurate near-nucleus derivatives from 3D grid data (see ConstantRef docstring)
• A custom SumOfRadials type when atomic contributions are available

Fields

• reference: callable providing ρ₀(x), ∂ρ₀/∂xξ, and ∂²ρ₀/∂xξ∂xζ

ConstantRef(val)

A callable that returns val for value queries and zero(val) for any derivative query. Use as reference_interp when building a log-density PHS whose reference density is a constant — typically ConstantRef(1.0) so that the stored data becomes log(data) and evaluations return exp(f) with accurate derivatives via the PHS chain rule.

This enables accurate ρ₀ derivatives from 3D grid data by building a log-space PHS of ρ₀ (where log(ρ₀) is smooth near nuclei) instead of a plain cubic spline (which oscillates near nuclei):

# log-space PHS of ρ₀ — stores log(ρ₀), evals return ρ₀ and ∂ρ₀/∂x accurately
itp_rho0 = phs_interp(grids, rho0; stencil_size=8, degree=3,
                      reference_interp = ConstantRef(1.0))
# main PHS of ρ: reference derivatives now come from the accurate log-space PHS
itp_phs  = phs_interp(grids, rho;  stencil_size=8, degree=3,
                      reference_interp = itp_rho0)