API

interpolate(method, points, samples; metric = Euclidean(), returnRBFmatrix = false, smooth = false, linsolve = nothing, sparse = :auto, neighbors = 50)

Create an interpolation of the data in samples sampled at the locations defined in points based on the interpolation method method. metric is any of the metrics defined by the Distances package. The RBF matrix used for solving the weights can be returned with the boolean returnRBFmatrix. Note that this option is only valid for RadialBasisFunction interpolations. linsolve selects the linear-solver algorithm from LinearSolve.jl used to solve for the weights; nothing uses LinearSolve's default algorithm.

points should be an $n×k$ matrix, where $n$ is dimension of the sampled space and $k$ is the number of points. This means that each column in the matrix defines one point.

samples is an $k×m$ array, where $k$ is the number of sampled points (same as for points) and $m$ is the dimension of the sampled data.

The RadialBasisFunction interpolation supports the use of unique RBF functions and widths for each sampled point by supplying method with a vector of interpolation methods of length $k$.

The RadialBasisFunction interpolation also supports smoothing of the data points using ridge regression. All points can be smoothed equally supplying a scalar value, alternatively each point can be smoothed independently by supplying a vector of smoothing values. Note that it is no longer interpolating when using smoothing.

The returned ScatteredInterpolant object can be passed to evaluate to interpolate the data to new points.

For compactly supported basis functions (Wendland), sparse controls a sparse interpolation path that assembles only the nonzero entries of the RBF matrix and solves with sparse Cholesky (CHOLMODFactorization): :auto (default) uses it when the problem is large and the matrix sufficiently sparse, true forces it (erroring if the kernel or metric makes sparsity impossible), false always uses the dense path. Globally supported kernels (Gaussian, Multiquadratic, …) are nonzero at every distance and can never be sparse — use PartitionOfUnity for locality with those.

A Wendland kernel constructed without ε gets its support radius from the data: 1/ε is set to the median distance to the neighbors-th nearest neighbor. An explicit ε always wins, in which case neighbors is ignored.

evaluate(itp, points)

Evaluate an interpolation object itp at the locations defined in points.

points should be an $n×k$ matrix, where $n$ is dimension of the sampled space and $k$ is the number of points. This means that each column in the matrix defines one point.

addpoints!(itp, points, samples)

Add new data points to an existing PartitionOfUnity interpolant without a full rebuild: only the local systems of the patches covering the new points are re-solved (using the same smooth and linsolve settings as the original build, and re-tuning the shape parameter of affected patches when the interpolant was built with tune = :loocv).

The patch grid is fixed at construction, so every new point must lie inside the bounding box of the original data and inside at least one existing patch; otherwise an ArgumentError suggests rebuilding with interpolate. Interpolants built with a per-point smoothing vector are not supported (the smoothing values of the new points would be ambiguous) — rebuild instead.

points and samples follow the same layout as in interpolate; samples must have the same second dimension as the original sample data. Returns itp.

Multiquadratic(ε = 1)

Define a Multiquadratic Radial Basis Function

\[ϕ(r) = \sqrt{1 + (εr)^2}\]

InverseMultiquadratic(ε = 1)

Define an Inverse Multiquadratic Radial Basis Function

\[ϕ(r) = \frac{1}{\sqrt{1 + (εr)^2}}\]

Gaussian(ε = 1)

Define a Gaussian Radial Basis Function

\[ϕ(r) = e^{-(εr)^2}\]

InverseQuadratic(ε = 1)

Define an Inverse Quadratic Radial Basis Function

\[ϕ(r) = \frac{1}{1 + (εr)^2}\]

Polyharmonic(k = 1)

Define a Polyharmonic Spline Radial Basis Function

\[ϕ(r) = r^k, k = 1, 3, 5, ... \\ ϕ(r) = r^k ln(r), k = 2, 4, 6, ...\]

Polyharmonic splines are only conditionally positive definite. The plain interpolation matrix has a zero diagonal ($ϕ(0) = 0$) and is indefinite, so the system is not guaranteed to be solvable and the interpolant does not reproduce polynomial trends. For a well-posed system that also reproduces low-order polynomials, use GeneralizedPolyharmonic, which augments the system with a polynomial term.

ThinPlate()

Define a Thin Plate Spline Radial Basis Function

\[ϕ(r) = r^2 ln(r)\]

This is a shorthand for Polyharmonic(2). As with Polyharmonic, the plain thin plate system is only conditionally positive definite; use GeneralizedPolyharmonic for a well-posed system with polynomial reproduction.

GeneralizedMultiquadratic(ε, β, degree)

Define a generalized Multiquadratic Radial Basis Function

\[ϕ(r) = (1 + (ε*r)^2)^β\]

Results in a positive definite system for a 'degree' of ⌈β⌉ or higher.

GeneralizedPolyharmonic(k, degree)

Define a generalized Polyharmonic Radial Basis Function

\[ϕ(r) = r^k, k = 1, 3, 5, ... \\ ϕ(r) = r^k ln(r), k = 2, 4, 6, ...\]

Results in a positive definite system for a 'degree' of ⌈k/2⌉ or higher for k = 1, 3, 5, ... and of exactly k + 1 for k = 2, 4, 6, ...

Wendland(dim, degree; ε = nothing)

Define a compactly supported Wendland Radial Basis Function. The resulting kernel is positive definite for data of dimension up to dim, is 2*degree times continuously differentiable, and vanishes identically for r ≥ 1/ε.

\[ϕ(r) = \max(1 - εr, 0)^{α(v, m)} \, f_{v,m}(εr)\]

where $v$ = degree ∈ {0, 1, 2, 3}, $m$ = dim, and $f_{v,m}$ is a polynomial of degree $v$. The kernel evaluation is provided by KernelFunctions.PiecewisePolynomialKernel; see its documentation for the exact coefficients. For dim ∈ {2, 3} and degree = 1 this is the classic C² function $ϕ(r) = (1 - εr)_+^4 (4εr + 1)$.

ε may be omitted (or passed explicitly as nothing), leaving the shape parameter unset. An unset ε cannot be evaluated or queried for support_radius directly; it must first be resolved from data via resolveshape (used internally by interpolate with the neighbors keyword), which sets ε to the reciprocal of the median distance to each point's neighbors-th nearest neighbor.

CompactSupportRBFInterpolant

Interpolant produced by interpolate for a compactly supported radial basis function (e.g. Wendland). Holds the solved weights w, the training points, the rbf kernel, the metric used to compute distances, and the tree (a NearestNeighbors.KDTree) used for the sparse assembly and reused by evaluate.

PartitionOfUnity(method; pointsperpatch = 80, overlap = 1.5, weight = nothing,
                 tune = :none)

Radial basis function partition of unity method (RBF-PUM) for large datasets. The data bounding box is covered with a regular grid of overlapping spherical patches; a small dense RBF interpolant using method (any radial basis function) is solved per patch, and patches are blended with smooth partition-of-unity weights, preserving exact interpolation at the data points.

pointsperpatch sets the targeted average number of points per patch (the trade-off between many small solves and few large ones). overlap inflates the patch radius relative to the grid cell half-diagonal and must be greater than 1 so the patches cover the whole bounding box; larger values increase smoothness of the blend regions at a higher evaluation cost. weight is the partition-of-unity weight function ψ(r) with support on [0, 1]; the default nothing selects the C² Wendland function of the data dimension at interpolate time.

tune selects automatic per-patch shape-parameter tuning: :none uses method as given, while :loocv chooses each patch's shape parameter ε by exact leave-one-out cross-validation (Rippa's method), searching candidates scaled to the patch's mean nearest-neighbor spacing plus the shape parameter stored in method itself (used as a fallback reference, not ignored); the kernel must have a shape parameter to tune (Polyharmonic, ThinPlate and GeneralizedPolyharmonic do not).

Performance: tune = :loocv is substantially slower to build than :none — each candidate evaluated costs about as much as solving the whole patch (Rippa's formula needs the full diagonal of the inverse, which is the same asymptotic cost as the factorization), and several candidates are evaluated per patch. Measured on a target-scale benchmark (n = 100,000, 3D): build time increased by roughly two orders of magnitude (~150×) relative to :none. Only use it where a single fixed ε is genuinely mis-scaled somewhere in the domain (see the manual for a worked example); on data where a well-chosen fixed ε is already close to optimal everywhere, tuning buys little to no accuracy for a large build-time cost.

Only the Euclidean metric is supported. The smooth and linsolve keywords of interpolate are forwarded to every per-patch solve. Additional points can be added to the returned interpolant with addpoints!.

Shepard(P = 2)

Standard Shepard interpolation with power parameter P.

NearestNeighbor()

Nearest neighbor interpolation.