SoapFilm

Soap film smooth for 2-D domains with complex boundaries.

Usage

SoapFilm(
    *,
    boundary=None,
    knots=None,
    k=30,
)

Implements the soap-film smoother of Wood, Bravington & Hedley (2008). Ordinary 2-D smooths such as TPRS or tensor-product splines treat the covariate domain as if it were convex and unobstructed; when the true domain has holes, concave coastlines, peninsulas, or other complicated shapes, such smooths will “leak” information across boundaries that are close in Euclidean distance but far apart when you have to travel around the domain (e.g. two points on opposite banks of a narrow bay). SoapFilm avoids this by finite-element-discretizing the domain itself: the smooth is represented as a piecewise-linear function on a triangulation of the actual (possibly non-convex, possibly multiply-connected) region, so its value at any point only depends on interior knots reachable through the domain. Use SoapFilm instead of a standard 2-D smooth whenever the covariates are genuinely spatial coordinates and the region they live in is not simply a rectangle — for example coastal, riverine, or other geographically constrained data.

Parameters

boundary: list[NDArray] | None = None: List of boundary loops, each an (m, 2) array of ordered vertices tracing a closed polygon. The first loop is the outer boundary of the domain; any additional loops are holes cut out of it (e.g. islands or excluded regions). If not supplied, a padded rectangular bounding box around the training data is used, which reduces the smooth to an ordinary (simply-connected, convex) domain — supply an explicit boundary whenever the domain has a non-trivial shape.
knots: NDArray | None = None: Interior knot locations as an (nk, 2) array; these become the nodes of the finite-element triangulation and directly determine the basis dimension. If not supplied, a roughly square grid of candidate points is generated over the bounding box, filtered to those lying inside the domain (outer boundary minus holes), and then subsampled down to at most k points. Supplying knots explicitly gives more control over their placement, which matters near sharp domain features.
k: int = 30: Target number of basis functions. If knots is not supplied, this determines the density of the automatically generated knot grid, and the actual number of basis functions equals the number of interior knots retained (which may be less than k). If knots is supplied directly, k is not used to size the basis; the number of basis functions equals len(knots). The default is 30.

Notes

Fitting proceeds in three stages:

Triangulation. Interior knots are combined with points sampled along the boundary loops (each boundary segment contributes its endpoint and midpoint) and a Delaunay triangulation is built over all of these points.
Finite-element assembly. On each triangle, standard piecewise-linear (barycentric) basis functions phi_i are used to assemble a stiffness matrix K (with entries K_{ij} = integral of grad(phi_i) . grad(phi_j) over the domain) and a mass matrix M (with entries M_{ij} = integral of phi_i * phi_j), then both are restricted to the rows and columns corresponding to interior knots (boundary points are not free parameters).
Basis evaluation. At an arbitrary evaluation point, the containing triangle is located (via Delaunay.find_simplex) and the point’s barycentric coordinates within that triangle give its basis-function weights; points outside the triangulation fall back to a nearest-knot indicator.

The penalty matrix is exactly the interior-restricted stiffness matrix,

\mathbf{S} = \mathbf{K}_{\text{interior}}, \qquad \boldsymbol{\beta}^\top \mathbf{S} \boldsymbol{\beta} = \int_\Omega \lVert \nabla f \rVert^2 \, dA,

the discretized Dirichlet energy (membrane/thin-film bending energy) of the fitted surface over the domain Omega, consistent with the “soap film” interpretation: the fitted surface behaves like a soap film stretched across the (possibly perforated) domain boundary. S is positive semi-definite; its null space corresponds to the constant function, so null_space_dimension() is 0 here by convention (the constant is handled via identifiability_constraints() rather than being excluded from the penalty). Because the basis is piecewise-linear on a triangulation rather than smooth in the classical sense, the fitted surface is continuous but only once-differentiable, and the quality of the fit is sensitive to the density and placement of interior knots relative to the sharpness of the domain’s boundary features — very thin or highly concave regions may need denser knots near the constriction to avoid leakage.

Examples

import numpy as np
from whittaker.smooths import SoapFilm

rng = np.random.default_rng(0)
x = rng.uniform(0, 1, (100, 2))

basis = SoapFilm(k=20).fit(x)
B = basis.basis_matrix(x)
S = basis.penalty_matrix()
B.shape, S.shape

((100, 20), (20, 20))

Attributes

Name	Description
n_basis	Number of basis functions in the fitted soap film basis.

n_basis

Number of basis functions in the fitted soap film basis.

n_basis: int

Equal to the number of interior knots used during fit(), which may be fewer than the k requested at construction time if automatic knot placement retained fewer points, or exactly len(knots) if knots were supplied explicitly.

Methods

Name	Description
basis_matrix()	Evaluate the soap film basis at `x`.
fit()	Fit the soap film smooth to training data `x`.
identifiability_constraints()	Return the mean-constraint row for the intercept.
null_space_dimension()	Return `0` by convention for the soap film penalty.
penalty_matrix()	Return the `k x k` penalty matrix `S`, the interior-restricted stiffness matrix.

basis_matrix()

Evaluate the soap film basis at x.

Usage

Source

basis_matrix(x)

For each point, locates the containing triangle of the fitted triangulation and returns its barycentric coordinates as basis-function weights against the interior knots; points that fall outside the triangulation (e.g. slightly outside the training domain) fall back to a one-hot indicator of the nearest interior knot.

Parameters

x: NDArray: Evaluation points. Shape (n, 2).

Returns

NDArray: Design matrix of shape (n, k) where k is the number of interior knots.

fit()

Fit the soap film smooth to training data x.

Usage

Source

fit(x)

Determines the domain boundary and interior knots (if not already supplied), triangulates the domain, and assembles the finite-element stiffness and mass matrices.

Parameters

x: NDArray: Training covariates. Shape (n, 2); SoapFilm only supports exactly two covariates (spatial x/y coordinates).

Returns

SoapFilm: Returns self for method chaining.

Raises

ValueError: If x does not have exactly 2 columns.

identifiability_constraints()

Return the mean-constraint row for the intercept.

Usage

Source

identifiability_constraints()

Returns

NDArray: A (1, k) matrix (uniform weights 1 / k) whose product with the coefficient vector is zero when the smooth has mean zero over the interior knots.

null_space_dimension()

Return 0 by convention for the soap film penalty.

Usage

Source

null_space_dimension()

The stiffness matrix’s true null space is spanned by the constant function, but SoapFilm handles the constant via identifiability_constraints() rather than excluding it from the penalty, so this method reports 0.

Returns

int: Always 0.

penalty_matrix()

Return the k x k penalty matrix S, the interior-restricted stiffness matrix.

Usage

Source

penalty_matrix()

S is the finite-element discretization of the Dirichlet energy integral of ||grad(f)||^2 over the domain, restricted to the interior knot degrees of freedom.

Returns

NDArray: Shape (k, k), positive semi-definite.