gam_terms/basis/sphere_spec.rs
1//! Intrinsic S² (sphere) smooth specification types.
2//!
3//! These configuration types own the user-facing description of an intrinsic
4//! sphere smooth (construction method, reproducing kernel, penalty order, and
5//! realized-design identifiability policy). The sibling `sphere_kernels` and
6//! `sphere_spectral` modules consume them to assemble the basis and penalty.
7
8use ndarray::Array2;
9use serde::{Deserialize, Serialize};
10
11use super::CenterStrategy;
12
13/// Which intrinsic S² smooth construction to use.
14///
15/// - `Wahba`: reproducing-kernel basis on a center set selected by
16/// `center_strategy`. The kernel function itself is chosen via
17/// `wahba_kernel` (Sobolev = true Σ [l(l+1)]^{-m} P_l kernel,
18/// Pseudo = Wahba's 1981 closed-form pseudo-spline).
19/// - `Harmonic`: real spherical-harmonic truncation up to `max_degree` with
20/// the Laplace-Beltrami eigenvalue penalty `[l(l+1)]^m`. Basis
21/// dimension is `max_degree * (max_degree + 2)`; centers are ignored.
22#[derive(Debug, Default, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
23pub enum SphereMethod {
24 #[default]
25 Wahba,
26 Harmonic,
27}
28
29/// Which reproducing kernel to use for `SphereMethod::Wahba`.
30///
31/// Both options yield positive-definite reproducing kernels on S² with
32/// the same family of `penalty_order m ∈ {1, 2, 3, 4}`. They define
33/// *different* RKHS, however:
34///
35/// - **Sobolev** (default, more correct): true Wahba/Sobolev kernel
36/// `K_m(γ) = (1/4π) Σ_{l ≥ 1} (2l+1) · [l(l+1)]^{-m} · P_l(cos γ)`,
37/// the reproducing kernel of `H^m(S²)` under the Laplace–Beltrami
38/// inner product. Penalty quadratic form recovers
39/// `‖f‖²_{H^m} = Σ_l [l(l+1)]^m · |f̂_l|²`.
40/// For m=1, 2, 3 this is evaluated via the closed forms from
41/// Beatson & zu Castell (2018) "Thinplate Splines on the Sphere"
42/// (SIGMA 14 (2018), 083) using elementary functions plus the
43/// di/trilogarithm. For m=4 we fall back to the spectral Legendre
44/// series (96 terms ⇒ truncation error ≲ 1e-12).
45///
46/// - **Pseudo**: Wahba 1981's "pseudo-spline" kernel with Legendre
47/// weights `2 / [(l+1)(l+2)···(l+m+1)]` (decaying as `l^{-(m+1)}`,
48/// different from Sobolev's `l^{-2m}`). Faster to evaluate (one
49/// elementary polynomial in `sin(γ/2)`, `log`), and matches mgcv's
50/// `bs="sos"` exactly. At `m=4` the pseudo-spline produces
51/// numerically tiny kernel values (`K(p,p) ≈ 3e-4`) and the basis
52/// pipeline historically collapsed the smooth contribution to zero
53/// — the cure is the REML scale-invariance fix in the solver, not
54/// the kernel itself.
55///
56/// Default is **Sobolev** because it matches the canonical "Wahba
57/// spline of order m on S²" interpretation. Use `Pseudo` for exact
58/// mgcv compatibility.
59#[derive(Debug, Default, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
60pub enum SphereWahbaKernel {
61 /// True Sobolev `H^m(S²)` reproducing kernel via closed-form
62 /// (elementary + Li_n) for m=1,2,3 and a deep (L=4096..96) spectral
63 /// series for m=4. This is the user-facing default.
64 #[default]
65 Sobolev,
66 /// Wahba 1981 closed-form pseudo-spline (mgcv `bs="sos"` compatible).
67 Pseudo,
68 /// Finite truncated-spectral Sobolev kernel evaluated by the same
69 /// Legendre 3-term recurrence the GPU `s2_wahba_legendre_colmajor`
70 /// kernel runs in registers:
71 ///
72 /// `K_L^{Sobolev}(γ) = Σ_{ℓ=1..L} c_ℓ · P_ℓ(cos γ)`,
73 /// `c_0 = 0`, `c_ℓ = (2ℓ+1) / (4π · [ℓ(ℓ+1)]^m)`.
74 ///
75 /// This is the explicit parity target for the GPU kernel: at any
76 /// finite `lmax` the closed-form (or m=4-spectral-at-fixed-deep-L)
77 /// variants differ by the documented truncation error, while the GPU
78 /// matches this variant exactly to roundoff. Single-source: both CPU
79 /// and GPU use the same Legendre recurrence and the same c_ℓ array.
80 SobolevTruncated {
81 /// Largest Legendre degree retained (≥ 1). Practical range
82 /// 5..200; GPU kernel uses LMAX as a compile-time `#define`.
83 lmax: u16,
84 },
85 /// Finite truncated-spectral Wahba-1981 pseudo-spline kernel:
86 ///
87 /// `K_L^{Pseudo}(γ) = Σ_{ℓ=1..L} c_ℓ · P_ℓ(cos γ)`,
88 /// `c_0 = 0`, `c_ℓ = 2 / (4π · Π_{k=1..m+1}(ℓ + k))`.
89 ///
90 /// Same role as [`Self::SobolevTruncated`] for the pseudo branch.
91 PseudoTruncated {
92 /// Largest Legendre degree retained (≥ 1).
93 lmax: u16,
94 },
95}
96
97/// Intrinsic S² (sphere) smooth configuration.
98#[derive(Debug, Clone, Serialize, Deserialize)]
99pub struct SphericalSplineBasisSpec {
100 /// Center/knot selection strategy in latitude/longitude coordinates
101 /// (used only when `method == Wahba`).
102 pub center_strategy: CenterStrategy,
103 /// Sphere roughness penalty order m ∈ {1, 2, 3, 4}.
104 ///
105 /// - **Wahba method, Sobolev kernel (default)**: the reproducing-kernel
106 /// norm is `Σ_l [l(l+1)]^m · |f̂_l|²` — the canonical `H^m(S²)` Sobolev
107 /// norm. `m=2` is the usual curvature (thin-plate analogue) penalty.
108 /// - **Wahba method, Pseudo kernel**: the order maps to the Wahba 1981
109 /// pseudo-spline kernel with Legendre weights
110 /// `2 / [(l+1)(l+2)···(l+m+1)]`. Same `m=2` is the TPS pseudo-spline.
111 /// - **Harmonic method**: raises the Laplace-Beltrami eigenvalue
112 /// penalty to `[l(l+1)]^m` (same as the Sobolev kernel norm above).
113 pub penalty_order: usize,
114 /// Add a ridge-like shrinkage penalty.
115 pub double_penalty: bool,
116 /// Interpret latitude/longitude in radians instead of degrees. Default is
117 /// degrees (Earth/data-frame convention); set true for radians.
118 #[serde(default)]
119 pub radians: bool,
120 /// Construction method. Default Wahba (reproducing kernel); Harmonic uses
121 /// real spherical harmonics with a Laplace-Beltrami eigenvalue penalty
122 /// (Wood §5.6.2, mgcv `bs="sos"`).
123 #[serde(default)]
124 pub method: SphereMethod,
125 /// Maximum spherical-harmonic degree L when `method == Harmonic`. Basis
126 /// dimension is `L * (L + 2)`. Ignored for Wahba.
127 #[serde(default)]
128 pub max_degree: Option<usize>,
129 /// When `method == Wahba`, which reproducing kernel to use:
130 /// Sobolev (true `H^m(S²)`, the default) or Pseudo
131 /// (Wahba 1981 / mgcv `bs="sos"`). See `SphereWahbaKernel` docs.
132 ///
133 /// Deserialised specs that omit this field fall back to the standard
134 /// `Default` (Sobolev). Saved models that pre-date this field are not
135 /// supported — they must be re-fit with the current encoding.
136 #[serde(default)]
137 pub wahba_kernel: SphereWahbaKernel,
138 /// Realized-design identifiability policy for the Wahba sphere (#532).
139 ///
140 /// The finite-center Wahba kernel design `K(data, centers) · z` can span a
141 /// near-constant direction on the data rows even though the continuous
142 /// kernel omits the l=0 mode.
143 /// In any model with a parametric intercept this collides with the global
144 /// intercept — the #531 constant-vs-intercept rank-1 collision class. The
145 /// fix (mirroring `MaternIdentifiability::FrozenTransform`) lets the global
146 /// identifiability pipeline compose a parametric-orthogonalization onto `z`
147 /// and freeze the composed transform here, so the predict-time rebuild
148 /// reuses the exact fit-time realized transform instead of silently dropping
149 /// the orthogonalization.
150 #[serde(default)]
151 pub identifiability: SphericalSplineIdentifiability,
152}
153
154/// Realized-design identifiability policy for the Wahba spherical spline (#532).
155#[derive(Debug, Clone, Serialize, Deserialize, Default)]
156pub enum SphericalSplineIdentifiability {
157 /// Fit-time default: keep the raw center coefficient chart, then let the
158 /// global identifiability pipeline residualize the realized design against
159 /// the parametric block.
160 #[default]
161 CenterSumToZero,
162 /// Predict-time replay: use this frozen realized-design transform directly
163 /// (the composed raw-chart/parametric transform captured at fit time).
164 /// `transform.nrows()` equals the raw Wahba basis width after its
165 /// finite-center kernel / low-degree harmonic decomposition;
166 /// `transform.ncols()` is the constrained smooth dimension.
167 FrozenTransform { transform: Array2<f64> },
168}
169
170impl Default for SphericalSplineBasisSpec {
171 fn default() -> Self {
172 Self {
173 center_strategy: CenterStrategy::FarthestPoint { num_centers: 50 },
174 penalty_order: 2,
175 double_penalty: true,
176 radians: false,
177 method: SphereMethod::Wahba,
178 max_degree: None,
179 wahba_kernel: SphereWahbaKernel::Sobolev,
180 identifiability: SphericalSplineIdentifiability::CenterSumToZero,
181 }
182 }
183}