# SoapFilm Soap film smooth for 2-D domains with complex boundaries. Usage ``` python SoapFilm( *, boundary=None, knots=None, k=30, ) ``` Implements the soap-film smoother of Wood, Bravington & Hedley (2008). Ordinary 2-D smooths such as [TPRS](TPRS.md#whittaker.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](SoapFilm.md#whittaker.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](SoapFilm.md#whittaker.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: 1. **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. 2. **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). 3. **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 ``` python 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](#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()](#basis_matrix) | Evaluate the soap film basis at `x`. | | [fit()](#fit) | Fit the soap film smooth to training data `x`. | | [identifiability_constraints()](#identifiability_constraints) | Return the mean-constraint row for the intercept. | | [null_space_dimension()](#null_space_dimension) | Return `0` by convention for the soap film penalty. | | [penalty_matrix()](#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 ``` python 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 ``` python 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](SoapFilm.md#whittaker.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 ``` python 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 ``` python null_space_dimension() ``` The stiffness matrix's true null space is spanned by the constant function, but [SoapFilm](SoapFilm.md#whittaker.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 ``` python 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.