gam: add mgcv Gaussian-process smooth parity

README.md

@@ -125,7 +125,7 @@ result, and prediction interfaces.
 
 | Formula surface    | Supported terms                                                                                      |
 | ------------------ | ---------------------------------------------------------------------------------------------------- |
-| Univariate smooths | `s(..., bs='bs')`, `cr`, `cs`, `cc`, `cp`, `ds`, `ps`, `tp`, `ts`                                    |
+| Univariate smooths | `s(..., bs='bs')`, `cr`, `cs`, `cc`, `cp`, `ds`, `gp`, `ps`, `tp`, `ts`                              |
 | Structured smooths | random effects `re`, factor smooths `fs`, sum-to-zero factor smooths `sz`                            |
 | Tensor products    | `te(...)` and `ti(...)` over supported numeric marginals                                             |
 | Parametric terms   | numeric and factor terms, supported interactions, intercept policies, and formula offsets            |

docs/generate_notebooks.py

@@ -1514,6 +1514,7 @@ def gam_notebook() -> dict:
 | `cr`, `cs` | cubic regression spline; `cs` adds null-space shrinkage |
 | `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 |
 | `tp`, `ts` | thin-plate regression spline; `ts` adds shrinkage |
 | `te(...)` | scale-invariant tensor product including main-effect directions |
 | `ti(...)` | tensor interaction with marginal main-effect directions removed |
@@ -1535,6 +1536,9 @@ def gam_notebook() -> dict:
     "shrinkage_cubic": GAM(formula="demand ~ s(temperature, bs='cs', k=10)"),
     "cyclic": GAM(formula="demand ~ s(hour, bs='cc', k=8)"),
     "p_spline": GAM(formula="demand ~ s(temperature, bs='ps', k=10)"),
+    "gaussian_process": GAM(
+        formula="demand ~ s(temperature, humidity, bs='gp', k=20, m=[3, 1.0])"
+    ),
     "thin_plate": GAM(formula="demand ~ s(temperature, bs='tp', k=10)"),
     "shrinkage_thin_plate": GAM(formula="demand ~ s(temperature, bs='ts', k=10)"),
     "tensor_surface": GAM(

docs/notebooks/01_gam.ipynb

@@ -629,6 +629,7 @@
     "| `cr`, `cs` | cubic regression spline; `cs` adds null-space shrinkage |\n",
     "| `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",
     "| `tp`, `ts` | thin-plate regression spline; `ts` adds shrinkage |\n",
     "| `te(...)` | scale-invariant tensor product including main-effect directions |\n",
     "| `ti(...)` | tensor interaction with marginal main-effect directions removed |\n",
@@ -652,6 +653,7 @@
        "['cubic',\n",
        " 'cyclic',\n",
        " 'factor_smooth',\n",
+       " 'gaussian_process',\n",
        " 'p_spline',\n",
        " 'random_effect',\n",
        " 'shrinkage_cubic',\n",
@@ -675,6 +677,9 @@
     "    \"shrinkage_cubic\": GAM(formula=\"demand ~ s(temperature, bs='cs', k=10)\"),\n",
     "    \"cyclic\": GAM(formula=\"demand ~ s(hour, bs='cc', k=8)\"),\n",
     "    \"p_spline\": GAM(formula=\"demand ~ s(temperature, bs='ps', k=10)\"),\n",
+    "    \"gaussian_process\": GAM(\n",
+    "        formula=\"demand ~ s(temperature, humidity, bs='gp', k=20, m=[3, 1.0])\"\n",
+    "    ),\n",
     "    \"thin_plate\": GAM(formula=\"demand ~ s(temperature, bs='tp', k=10)\"),\n",
     "    \"shrinkage_thin_plate\": GAM(formula=\"demand ~ s(temperature, bs='ts', k=10)\"),\n",
     "    \"tensor_surface\": GAM(\n",

nampy/gam/compiler/factory.py

@@ -18,6 +18,7 @@
 from ..smooths.univariate.bs import DerivativeBSplineTerm1D
 from ..smooths.univariate.cr import CubicSplineTerm
 from ..smooths.univariate.ds import DuchonSplineTerm
+from ..smooths.univariate.gp import GaussianProcessTerm
 from ..smooths.univariate.ps import PSplineTerm1D
 from ..specs import LinearPredictorSpec, PenaltyGroupSpec, TermSpec
 from ..specs.smooth import (
@@ -27,6 +28,7 @@
     DerivativeBSplineSmoothSpec,
     DuchonSplineSmoothSpec,
     FactorSmoothInteractionSpec,
+    GaussianProcessSmoothSpec,
     PSplineSmoothSpec,
     RandomEffectSmoothSpec,
     ShapeConstrainedSmoothSpec,
@@ -207,6 +209,25 @@ def instantiate_term(term_like: TermSpec | Any):
             metadata=metadata,
         )
 
+    if isinstance(smooth_spec, GaussianProcessSmoothSpec):
+        return GaussianProcessTerm(
+            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, ShapeConstrainedSmoothSpec):
         if len(features) == 2:
             return BivariateShapePSplineTerm(

nampy/gam/compiler/linked_basis.py

@@ -92,6 +92,17 @@ def attach_shared_basis_metadata(predictor_specs, X, feature_names):
                 "different feature sets; mgcv 1.9-4 fails while constructing "
                 "the shared tensor basis."
             )
+        has_scalar_pc = any(
+            term.smooth_spec is not None
+            and str(term.smooth_spec.special).lower() == "s"
+            and getattr(term.smooth_spec, "pc", None) is not None
+            for term in group_terms
+        )
+        if has_scalar_pc and len(feature_tuples) > 1:
+            raise NotImplementedError(
+                "pc= is not supported across id-linked s() terms with different "
+                "feature sets; mgcv 1.9-4 fails while constructing the shared basis."
+            )
         for term in group_terms[1:]:
             _clone_linked_smooth_spec(base_term, term)
 

nampy/gam/smooths/__init__.py

@@ -22,6 +22,7 @@
 from .univariate.bs import DerivativeBSplineTerm1D
 from .univariate.cr import CubicSplineTerm
 from .univariate.ds import DuchonSplineTerm
+from .univariate.gp import GaussianProcessTerm
 from .univariate.ps import PSplineTerm1D
 from .univariate.tp import ThinPlateSplineTerm
 
@@ -33,6 +34,7 @@
 bs = DerivativeBSplineTerm1D
 cr = cs = cc = CubicSplineTerm
 ds = DuchonSplineTerm
+gp = GaussianProcessTerm
 cp = ps = PSplineTerm1D
 tp = ts = ThinPlateSplineTerm
 fs = FSmoothInteractionTerm
@@ -60,6 +62,7 @@
     "build_penalty_definition",
     "CubicSplineTerm",
     "DuchonSplineTerm",
+    "GaussianProcessTerm",
     "DerivativeBSplineTerm1D",
     "PSplineTerm1D",
     "ThinPlateSplineTerm",
@@ -76,6 +79,7 @@
     "cs",
     "cc",
     "ds",
+    "gp",
     "cp",
     "ps",
     "tp",

nampy/gam/smooths/categorical/fs.py

@@ -18,6 +18,7 @@
 from ..univariate.bs import DerivativeBSplineTerm1D
 from ..univariate.cr import CubicSplineTerm
 from ..univariate.ds import DuchonSplineTerm
+from ..univariate.gp import GaussianProcessTerm
 from ..univariate.ps import PSplineTerm1D
 from .categorical_utils import (
     as_object_1d,
@@ -109,7 +110,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, ps, tp, ts
+    bs, cr, cs, cc, cp, ds, gp, ps, tp, ts
     """
     base_bs = str(base_bs).lower()
     metric_features = list(metric_features)
@@ -120,22 +121,23 @@ 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", "tp", "ts"}:
+    if len(metric_features) > 1 and base_bs not in {"ds", "gp", "tp", "ts"}:
         raise NotImplementedError(
             f"Current {mode} implementation supports multivariate base smooths only "
-            f"for bs in {{'ds','tp','ts'}}, got base bs={base_bs!r}."
+            f"for bs in {{'ds','gp','tp','ts'}}, got base bs={base_bs!r}."
         )
 
     if xt_rest is not None and base_bs not in {
         "bs",
         "cp",
         "ds",
+        "gp",
         "ps",
         "tp",
         "ts",
     }:
         raise NotImplementedError(
-            "Extra xt options are currently only supported for bs/cp/ds/ps/tp/ts "
+            "Extra xt options are currently only supported for bs/cp/ds/gp/ps/tp/ts "
             "base smooths, "
             f"got xt={xt_rest!r} with base bs={base_bs!r}."
         )
@@ -220,6 +222,24 @@ def _build_base_smooth_term(
             metadata=metadata,
         )
 
+    if base_bs == "gp":
+        return GaussianProcessTerm(
+            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 in {"tp", "ts"}:
         return make_smooth_term(
             base_bs,
@@ -242,7 +262,7 @@ def _build_base_smooth_term(
 
     raise NotImplementedError(
         f"Current {mode} implementation supports base bs in "
-        f"{{'bs','cr','cs','cc','cp','ds','ps','tp','ts'}}, got {base_bs!r}."
+        f"{{'bs','cr','cs','cc','cp','ds','gp','ps','tp','ts'}}, got {base_bs!r}."
     )
 
 
@@ -251,6 +271,8 @@ def _penalty_rank_from_base_term(base_term, basis_matrix, penalty_matrix) -> int
         return int(base_term._setup.ranks[0])
     if isinstance(base_term, DuchonSplineTerm):
         return int(base_term._setup.rank)
+    if isinstance(base_term, GaussianProcessTerm):
+        return int(base_term._setup.rank)
     if isinstance(base_term, PSplineTerm1D) and len(base_term.penalties) > 0:
         if str(base_term.basis_name).lower() == "cp":
             return int(base_term._setup.rank)

nampy/gam/smooths/tensor/marginals.py

@@ -10,11 +10,12 @@
 from ..univariate.bs import DerivativeBSplineTerm1D
 from ..univariate.cr import CubicSplineTerm
 from ..univariate.ds import DuchonSplineTerm
+from ..univariate.gp import GaussianProcessTerm
 from ..univariate.ps import PSplineTerm1D
 from ..univariate.tp import ThinPlateSplineTerm
 
 TENSOR_MARGINAL_BASES = frozenset(
-    {"bs", "cr", "cs", "cc", "cp", "ds", "ps", "tp", "ts"}
+    {"bs", "cr", "cs", "cc", "cp", "ds", "gp", "ps", "tp", "ts"}
 )
 
 
@@ -130,6 +131,22 @@ def make_tensor_marginal_term(
             metadata=metadata,
         )
 
+    if basis == "gp":
+        return GaussianProcessTerm(
+            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 in {"tp", "ts"}:
         return ThinPlateSplineTerm(
             feature=marginal_features,

nampy/gam/smooths/tensor/te.py

@@ -114,6 +114,12 @@ def __init__(
 
         self._by_state = None
 
+    @property
+    def expected_linked_penalty_count(self):
+        if self.select:
+            return None
+        return int(np.sum(~np.asarray(self.fixed_flags, dtype=bool)))
+
     def fit(self, X, feature_names):
         marginal_shared_setups = self._linked_id_marginal_setups(self.feature)
         marginals, _, _ = build_tensor_marginal_terms(

nampy/gam/smooths/tensor/ti.py

@@ -117,6 +117,12 @@ def __init__(
 
         self._by_state = None
 
+    @property
+    def expected_linked_penalty_count(self):
+        if self.select:
+            return None
+        return int(np.sum(~np.asarray(self.fixed_flags, dtype=bool)))
+
     def fit(self, X, feature_names):
         self._set_by_state(X, feature_names)
 

nampy/gam/smooths/univariate/__init__.py

@@ -1,20 +1,23 @@
 from .bs import DerivativeBSplineTerm1D
 from .cr import CubicSplineTerm
 from .ds import DuchonSplineTerm
+from .gp import GaussianProcessTerm
 from .ps import PSplineTerm1D
 from .tp import ThinPlateSplineTerm
 
 bs = DerivativeBSplineTerm1D
 cr = cs = cc = CubicSplineTerm
 ds = DuchonSplineTerm
+gp = GaussianProcessTerm
 cp = ps = PSplineTerm1D
 tp = ts = ThinPlateSplineTerm
 
 __all__ = [
     "DerivativeBSplineTerm1D",
     "CubicSplineTerm",
     "DuchonSplineTerm",
+    "GaussianProcessTerm",
     "PSplineTerm1D",
     "ThinPlateSplineTerm",
 ]
-__all__ += ["bs", "cr", "cs", "cc", "cp", "ds", "ps", "tp", "ts"]
+__all__ += ["bs", "cr", "cs", "cc", "cp", "ds", "gp", "ps", "tp", "ts"]