Add mgcv spherical spline parity

README.md

@@ -125,12 +125,18 @@ result, and prediction interfaces.
 
 | Formula surface    | Supported terms                                                                                      |
 | ------------------ | ---------------------------------------------------------------------------------------------------- |
-| Univariate smooths | `s(..., bs='bs')`, `cr`, `cs`, `cc`, `cp`, `ds`, `gp`, `ps`, `tp`, `ts`                              |
+| Metric smooths     | `s(..., bs='bs')`, `cr`, `cs`, `cc`, `cp`, `ds`, `gp`, `ps`, `sos`, `tp`, `ts`                       |
 | Structured smooths | Markov random fields `mrf`, random effects `re`, factor smooths `fs`, sum-to-zero factor smooths `sz` |
-| Tensor products    | `te(...)` and `ti(...)` over supported numeric or MRF marginals                                      |
+| Tensor products    | `te(...)` and `ti(...)` over supported metric or MRF marginals                                       |
 | Parametric terms   | numeric and factor terms, supported interactions, intercept policies, and formula offsets            |
 | Shared smoothing   | supported `id=` groups, fixed/free smoothing parameters, `select=True`, and `pc=` on supported bases |
 
+For spherical splines, write
+`s(latitude, longitude, bs='sos')`: coordinates are degrees in latitude-first,
+longitude-second order. A spherical marginal inside `te()` or `ti()` needs
+`d=2`. Because the array API creates one smooth per input column, joint SOS
+terms use the formula interface.
+
 
 ### Shape-constrained functionality
 

docs/api/gam.rst

@@ -30,6 +30,21 @@ High-level model
    model.summary()
    model.plot()
 
+Spherical smooths
+-----------------
+
+Use ``s(latitude, longitude, bs="sos")`` for an isotropic smooth on a
+sphere.  Coordinates are supplied in degrees, latitude first and longitude
+second.  The default ``m=0`` is the second-order Wendelberger spline; integer
+orders from ``-2`` through ``4`` select the upstream Duchon or Wahba kernel
+branches.  SOS margins in ``te`` and ``ti`` must be grouped with ``d=2``.
+The array API cannot express this joint two-coordinate term, so SOS models use
+the formula interface.  The ``m=-1`` null space is four-dimensional; combining
+that order with an ``fs`` factor smooth is rejected because the corresponding
+upstream penalty split is LAPACK-orientation dependent.  Upstream SOS smooths
+also do not define derivative matrices, and NAMpy does not expose the
+hemisphere-specific ``plot.gam`` schemes 0 and 1 through its generic plotter.
+
 Shape-constrained smooths
 -------------------------
 

docs/generate_notebooks.py

@@ -1515,6 +1515,7 @@ def gam_notebook() -> dict:
 | `cc` | cyclic cubic spline for periodic covariates |
 | `ps` | P-spline with difference penalties |
 | `gp` | low-rank Gaussian-process smooth with spherical, power-exponential, or Matérn covariance |
+| `sos` | isotropic spherical spline for latitude then longitude in degrees |
 | `mrf` | region effect coupled by a neighbor graph, polygon boundary, or supplied penalty |
 | `tp`, `ts` | thin-plate regression spline; `ts` adds shrinkage |
 | `te(...)` | scale-invariant tensor product including main-effect directions |
@@ -1540,6 +1541,9 @@ def gam_notebook() -> dict:
     "gaussian_process": GAM(
         formula="demand ~ s(temperature, humidity, bs='gp', k=20, m=[3, 1.0])"
     ),
+    "spherical_spline": GAM(
+        formula="demand ~ s(latitude, longitude, bs='sos', k=20, m=0)"
+    ),
     "markov_random_field": GAM(
         formula=(
             "demand ~ s(region, bs='mrf', "

docs/notebooks/01_gam.ipynb

@@ -630,6 +630,7 @@
     "| `cc` | cyclic cubic spline for periodic covariates |\n",
     "| `ps` | P-spline with difference penalties |\n",
     "| `gp` | low-rank Gaussian-process smooth with spherical, power-exponential, or Mat\u00e9rn covariance |\n",
+    "| `sos` | isotropic spherical spline for latitude then longitude in degrees |\n",
     "| `mrf` | region effect coupled by a neighbor graph, polygon boundary, or supplied penalty |\n",
     "| `tp`, `ts` | thin-plate regression spline; `ts` adds shrinkage |\n",
     "| `te(...)` | scale-invariant tensor product including main-effect directions |\n",
@@ -682,6 +683,9 @@
     "    \"gaussian_process\": GAM(\n",
     "        formula=\"demand ~ s(temperature, humidity, bs='gp', k=20, m=[3, 1.0])\"\n",
     "    ),\n",
+    "    \"spherical_spline\": GAM(\n",
+    "        formula=\"demand ~ s(latitude, longitude, bs='sos', k=20, m=0)\"\n",
+    "    ),\n",
     "    \"markov_random_field\": GAM(\n",
     "        formula=(\n",
     "            \"demand ~ s(region, bs='mrf', \"\n",

nampy/gam/compiler/factory.py

@@ -21,6 +21,7 @@
 from ..smooths.univariate.ds import DuchonSplineTerm
 from ..smooths.univariate.gp import GaussianProcessTerm
 from ..smooths.univariate.ps import PSplineTerm1D
+from ..smooths.univariate.sos import SphericalSplineTerm
 from ..specs import LinearPredictorSpec, PenaltyGroupSpec, TermSpec
 from ..specs.smooth import (
     CubicRegressionSmoothSpec,
@@ -34,6 +35,7 @@
     PSplineSmoothSpec,
     RandomEffectSmoothSpec,
     ShapeConstrainedSmoothSpec,
+    SphericalSplineSmoothSpec,
     SumToZeroFactorSmoothSpec,
     TensorInteractionSmoothSpec,
     TensorProductSmoothSpec,
@@ -230,6 +232,25 @@ def instantiate_term(term_like: TermSpec | Any):
             metadata=metadata,
         )
 
+    if isinstance(smooth_spec, SphericalSplineSmoothSpec):
+        return SphericalSplineTerm(
+            feature=features,
+            k=smooth_spec.k,
+            m=smooth_spec.m,
+            label=label,
+            term_id=term_like.term_id,
+            smoothing_id=smoothing_id,
+            by=by,
+            sp=smooth_spec.sp,
+            select=smooth_spec.select,
+            fixed=smooth_spec.fx,
+            constraint_mode=smooth_spec.constraint_mode,
+            pc=smooth_spec.pc,
+            knots=smooth_spec.knots,
+            xt=smooth_spec.xt,
+            metadata=metadata,
+        )
+
     if isinstance(smooth_spec, MarkovRandomFieldSmoothSpec):
         return MarkovRandomFieldTerm(
             feature=features,

nampy/gam/diagnostics/plots.py

@@ -190,6 +190,12 @@ def _prepare_smooth(model, tb, *, se, n, n2, n3, xlab, ylab, main, label,
     basis_name = str(getattr(tb, "basis_name", "")).lower()
     X_train = np.asarray(model.X_)
 
+    if basis_name == "sos":
+        raise NotImplementedError(
+            "plot() for bs='sos' requires mgcv's rotated hemisphere projection; "
+            "the generic rectangular 2D plot is not equivalent."
+        )
+
     if basis_name == "re":
         # plot.random.effect (plots.r:357-367): X is the identity; the plot
         # is a normal QQ plot of the estimated effects.

nampy/gam/smooths/__init__.py

@@ -25,6 +25,7 @@
 from .univariate.ds import DuchonSplineTerm
 from .univariate.gp import GaussianProcessTerm
 from .univariate.ps import PSplineTerm1D
+from .univariate.sos import SphericalSplineTerm
 from .univariate.tp import ThinPlateSplineTerm
 
 # mgcv-facing smooth aliases keep formulas/tests readable without reintroducing
@@ -36,6 +37,7 @@
 cr = cs = cc = CubicSplineTerm
 ds = DuchonSplineTerm
 gp = GaussianProcessTerm
+sos = SphericalSplineTerm
 mrf = MarkovRandomFieldTerm
 cp = ps = PSplineTerm1D
 tp = ts = ThinPlateSplineTerm
@@ -65,6 +67,7 @@
     "CubicSplineTerm",
     "DuchonSplineTerm",
     "GaussianProcessTerm",
+    "SphericalSplineTerm",
     "MarkovRandomFieldTerm",
     "DerivativeBSplineTerm1D",
     "PSplineTerm1D",
@@ -83,6 +86,7 @@
     "cc",
     "ds",
     "gp",
+    "sos",
     "mrf",
     "cp",
     "ps",

nampy/gam/smooths/categorical/fs.py

@@ -20,6 +20,7 @@
 from ..univariate.ds import DuchonSplineTerm
 from ..univariate.gp import GaussianProcessTerm
 from ..univariate.ps import PSplineTerm1D
+from ..univariate.sos import SphericalSplineTerm
 from .categorical_utils import (
     as_object_1d,
     factor_indicator_matrix,
@@ -111,7 +112,7 @@ def _build_base_smooth_term(
     Build the per-level base smooth used inside fs/sz.
 
     Supported base smooth classes in the current codebase:
-    bs, cr, cs, cc, cp, ds, gp, mrf, ps, tp, ts
+    bs, cr, cs, cc, cp, ds, gp, mrf, ps, sos, tp, ts
     """
     base_bs = str(base_bs).lower()
     metric_features = list(metric_features)
@@ -122,10 +123,10 @@ def _build_base_smooth_term(
     if mode == "fs" and base_bs in {"cs", "ts"}:
         raise NotImplementedError(_fs_full_rank_base_error(base_bs))
 
-    if len(metric_features) > 1 and base_bs not in {"ds", "gp", "tp", "ts"}:
+    if len(metric_features) > 1 and base_bs not in {"ds", "gp", "sos", "tp", "ts"}:
         raise NotImplementedError(
             f"Current {mode} implementation supports multivariate base smooths only "
-            f"for bs in {{'ds','gp','tp','ts'}}, got base bs={base_bs!r}."
+            f"for bs in {{'ds','gp','sos','tp','ts'}}, got base bs={base_bs!r}."
         )
 
     if xt_rest is not None and base_bs not in {
@@ -135,11 +136,12 @@ def _build_base_smooth_term(
         "gp",
         "mrf",
         "ps",
+        "sos",
         "tp",
         "ts",
     }:
         raise NotImplementedError(
-            "Extra xt options are currently only supported for bs/cp/ds/gp/mrf/ps/tp/ts "
+            "Extra xt options are currently only supported for bs/cp/ds/gp/mrf/ps/sos/tp/ts "
             "base smooths, "
             f"got xt={xt_rest!r} with base bs={base_bs!r}."
         )
@@ -242,6 +244,24 @@ def _build_base_smooth_term(
             metadata=metadata,
         )
 
+    if base_bs == "sos":
+        return SphericalSplineTerm(
+            feature=metric_features,
+            k=k,
+            m=outer_m,
+            label=label,
+            smoothing_id=None,
+            by=by,
+            sp=None,
+            select=bool(select),
+            fixed=bool(fixed),
+            constraint_mode=str(constraint_mode),
+            pc=None,
+            knots=knots,
+            xt=xt_rest,
+            metadata=metadata,
+        )
+
     if base_bs == "mrf":
         return MarkovRandomFieldTerm(
             feature=metric_features[0],
@@ -280,7 +300,7 @@ def _build_base_smooth_term(
 
     raise NotImplementedError(
         f"Current {mode} implementation supports base bs in "
-        f"{{'bs','cr','cs','cc','cp','ds','gp','mrf','ps','tp','ts'}}, got {base_bs!r}."
+        f"{{'bs','cr','cs','cc','cp','ds','gp','mrf','ps','sos','tp','ts'}}, got {base_bs!r}."
     )
 
 
@@ -291,6 +311,8 @@ def _penalty_rank_from_base_term(base_term, basis_matrix, penalty_matrix) -> int
         return int(base_term._setup.rank)
     if isinstance(base_term, GaussianProcessTerm):
         return int(base_term._setup.rank)
+    if isinstance(base_term, SphericalSplineTerm):
+        return int(base_term._setup.rank)
     if isinstance(base_term, MarkovRandomFieldTerm):
         return int(base_term._setup.rank)
     if isinstance(base_term, PSplineTerm1D) and len(base_term.penalties) > 0:
@@ -697,6 +719,15 @@ def fit(self, X, feature_names):
             )
 
         self._base_term = base_term
+        if (
+            isinstance(base_term, SphericalSplineTerm)
+            and int(base_term._setup.null_space_dim) > 1
+        ):
+            raise NotImplementedError(
+                "bs='fs' with an SOS m=-1 base is not enabled: mgcv's four-way "
+                "repeated null eigenspace receives separate penalties whose "
+                "orientation is LAPACK-dependent. Use another SOS order."
+            )
         if (
             isinstance(base_term, MarkovRandomFieldTerm)
             and base_term._setup.used_low_rank

nampy/gam/smooths/tensor/marginals.py

@@ -13,10 +13,11 @@
 from ..univariate.ds import DuchonSplineTerm
 from ..univariate.gp import GaussianProcessTerm
 from ..univariate.ps import PSplineTerm1D
+from ..univariate.sos import SphericalSplineTerm
 from ..univariate.tp import ThinPlateSplineTerm
 
 TENSOR_MARGINAL_BASES = frozenset(
-    {"bs", "cr", "cs", "cc", "cp", "ds", "gp", "mrf", "ps", "tp", "ts"}
+    {"bs", "cr", "cs", "cc", "cp", "ds", "gp", "mrf", "ps", "sos", "tp", "ts"}
 )
 
 
@@ -147,6 +148,27 @@ def make_tensor_marginal_term(
             metadata=metadata,
         )
 
+    if basis == "sos":
+        if len(marginal_features) != 2:
+            raise ValueError(
+                "Tensor marginal basis 'sos' requires a two-feature group "
+                "(latitude, longitude); supply d=2 for that marginal."
+            )
+        return SphericalSplineTerm(
+            feature=marginal_features,
+            k=k,
+            m=m,
+            xt=xt,
+            label=str(feature),
+            smoothing_id=None,
+            by=None,
+            select=False,
+            fixed=False,
+            constraint_mode=constraint_mode,
+            knots=knots,
+            metadata=metadata,
+        )
+
     if basis == "mrf":
         if len(marginal_features) != 1:
             raise ValueError("Tensor marginal basis 'mrf' only handles one feature.")

nampy/gam/smooths/univariate/__init__.py

@@ -3,12 +3,14 @@
 from .ds import DuchonSplineTerm
 from .gp import GaussianProcessTerm
 from .ps import PSplineTerm1D
+from .sos import SphericalSplineTerm
 from .tp import ThinPlateSplineTerm
 
 bs = DerivativeBSplineTerm1D
 cr = cs = cc = CubicSplineTerm
 ds = DuchonSplineTerm
 gp = GaussianProcessTerm
+sos = SphericalSplineTerm
 cp = ps = PSplineTerm1D
 tp = ts = ThinPlateSplineTerm
 
@@ -17,7 +19,8 @@
     "CubicSplineTerm",
     "DuchonSplineTerm",
     "GaussianProcessTerm",
+    "SphericalSplineTerm",
     "PSplineTerm1D",
     "ThinPlateSplineTerm",
 ]
-__all__ += ["bs", "cr", "cs", "cc", "cp", "ds", "gp", "ps", "tp", "ts"]
+__all__ += ["bs", "cr", "cs", "cc", "cp", "ds", "gp", "sos", "ps", "tp", "ts"]