feat(gam): add mgcv-parity Duchon splines

README.md

@@ -125,7 +125,7 @@ result, and prediction interfaces.
 
 | Formula surface    | Supported terms                                                                                      |
 | ------------------ | ---------------------------------------------------------------------------------------------------- |
-| Univariate smooths | `s(..., bs='bs')`, `cr`, `cs`, `cc`, `cp`, `ps`, `tp`, `ts`                                          |
+| Univariate smooths | `s(..., bs='bs')`, `cr`, `cs`, `cc`, `cp`, `ds`, `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            |

nampy/gam/compiler/factory.py

@@ -17,13 +17,15 @@
 from ..smooths.shape.scop import ShapeConstrainedPSplineTerm
 from ..smooths.univariate.bs import DerivativeBSplineTerm1D
 from ..smooths.univariate.cr import CubicSplineTerm
+from ..smooths.univariate.ds import DuchonSplineTerm
 from ..smooths.univariate.ps import PSplineTerm1D
 from ..specs import LinearPredictorSpec, PenaltyGroupSpec, TermSpec
 from ..specs.smooth import (
     CubicRegressionSmoothSpec,
     CubicShrinkageSmoothSpec,
     CyclicCubicRegressionSmoothSpec,
     DerivativeBSplineSmoothSpec,
+    DuchonSplineSmoothSpec,
     FactorSmoothInteractionSpec,
     PSplineSmoothSpec,
     RandomEffectSmoothSpec,
@@ -186,6 +188,25 @@ def instantiate_term(term_like: TermSpec | Any):
             metadata=metadata,
         )
 
+    if isinstance(smooth_spec, DuchonSplineSmoothSpec):
+        return DuchonSplineTerm(
+            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/smooths/__init__.py

@@ -21,6 +21,7 @@
 from .tensor.ti import InteractionTensorProductSplineTerm
 from .univariate.bs import DerivativeBSplineTerm1D
 from .univariate.cr import CubicSplineTerm
+from .univariate.ds import DuchonSplineTerm
 from .univariate.ps import PSplineTerm1D
 from .univariate.tp import ThinPlateSplineTerm
 
@@ -31,6 +32,7 @@
 
 bs = DerivativeBSplineTerm1D
 cr = cs = cc = CubicSplineTerm
+ds = DuchonSplineTerm
 cp = ps = PSplineTerm1D
 tp = ts = ThinPlateSplineTerm
 fs = FSmoothInteractionTerm
@@ -57,6 +59,7 @@
     "sync_by_state_attributes",
     "build_penalty_definition",
     "CubicSplineTerm",
+    "DuchonSplineTerm",
     "DerivativeBSplineTerm1D",
     "PSplineTerm1D",
     "ThinPlateSplineTerm",
@@ -72,6 +75,7 @@
     "cr",
     "cs",
     "cc",
+    "ds",
     "cp",
     "ps",
     "tp",

nampy/gam/smooths/categorical/fs.py

@@ -17,6 +17,7 @@
 from ..smooth_base import BaseSmoothTerm, by_values_from_new_data, column_as_object
 from ..univariate.bs import DerivativeBSplineTerm1D
 from ..univariate.cr import CubicSplineTerm
+from ..univariate.ds import DuchonSplineTerm
 from ..univariate.ps import PSplineTerm1D
 from .categorical_utils import (
     as_object_1d,
@@ -108,7 +109,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, ps, tp, ts
+    bs, cr, cs, cc, cp, ds, ps, tp, ts
     """
     base_bs = str(base_bs).lower()
     metric_features = list(metric_features)
@@ -119,16 +120,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 {"tp", "ts"}:
+    if len(metric_features) > 1 and base_bs not in {"ds", "tp", "ts"}:
         raise NotImplementedError(
             f"Current {mode} implementation supports multivariate base smooths only "
-            f"for bs in {{'tp','ts'}}, got base bs={base_bs!r}."
+            f"for bs in {{'ds','tp','ts'}}, got base bs={base_bs!r}."
         )
 
-    if xt_rest is not None and base_bs not in {"bs", "tp", "ts", "ps", "cp"}:
+    if xt_rest is not None and base_bs not in {
+        "bs",
+        "cp",
+        "ds",
+        "ps",
+        "tp",
+        "ts",
+    }:
         raise NotImplementedError(
-            "Extra xt options are currently only supported for bs/tp/ts/ps/cp base "
-            "smooths, "
+            "Extra xt options are currently only supported for bs/cp/ds/ps/tp/ts "
+            "base smooths, "
             f"got xt={xt_rest!r} with base bs={base_bs!r}."
         )
 
@@ -194,6 +202,24 @@ def _build_base_smooth_term(
             metadata=metadata,
         )
 
+    if base_bs == "ds":
+        return DuchonSplineTerm(
+            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,
@@ -216,13 +242,15 @@ def _build_base_smooth_term(
 
     raise NotImplementedError(
         f"Current {mode} implementation supports base bs in "
-        f"{{'bs','cr','cs','cc','cp','ps','tp','ts'}}, got {base_bs!r}."
+        f"{{'bs','cr','cs','cc','cp','ds','ps','tp','ts'}}, got {base_bs!r}."
     )
 
 
 def _penalty_rank_from_base_term(base_term, basis_matrix, penalty_matrix) -> int:
     if isinstance(base_term, DerivativeBSplineTerm1D):
         return int(base_term._setup.ranks[0])
+    if isinstance(base_term, DuchonSplineTerm):
+        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)
@@ -358,6 +386,28 @@ def __init__(
 
         self.skip_centering = True
 
+    @property
+    def expected_linked_penalty_count(self):
+        # The number depends on the fitted base null space (fs) or factor
+        # combinations (sz), so defer linked-group sizing until compilation.
+        return None
+
+    def _metric_knots_for_base(self):
+        knots = self.knots
+        if knots is None or not isinstance(knots, (list, tuple)):
+            return knots
+        if len(knots) != len(self._feature_names or []):
+            return knots
+        metric_names = set(self._metric_feature_names or [])
+        selected = [
+            value
+            for name, value in zip(self._feature_names, knots, strict=True)
+            if name in metric_names
+        ]
+        if len(selected) == 1:
+            return selected[0]
+        return selected
+
     @property
     def basis_train(self):
         if self._delegate_term is not None:
@@ -464,7 +514,7 @@ def _build_delegate_base_or_re(self, X, feature_names, *, default_bs, mode):
                 label=self.label,
                 fixed=self.fixed,
                 by=self.by,
-                knots=self.knots,
+                knots=self._metric_knots_for_base(),
                 xt_rest=base_spec.xt_rest,
                 outer_m=self.m,
                 mode=mode,
@@ -587,7 +637,7 @@ def fit(self, X, feature_names):
             label=self.label,
             fixed=self.fixed,
             by=None,
-            knots=self.knots,
+            knots=self._metric_knots_for_base(),
             xt_rest=base_spec.xt_rest,
             outer_m=self.m,
             mode="fs",
@@ -603,12 +653,13 @@ def fit(self, X, feature_names):
             )
 
         self._base_term = base_term
-        B0, S0, _ = self._base_constructor_fit_matrices()
-        B0 = np.asarray(B0, dtype=np.float64)
+        B_setup, S0, _ = self._base_constructor_fit_matrices()
+        B_setup = np.asarray(B_setup, dtype=np.float64)
+        B0 = np.asarray(self._base_constructor_predict_matrix(X), dtype=np.float64)
         S0 = np.asarray(S0, dtype=np.float64)
 
-        base_rank = _penalty_rank_from_base_term(base_term, B0, S0)
-        null_d = int(B0.shape[1] - base_rank)
+        base_rank = _penalty_rank_from_base_term(base_term, B_setup, S0)
+        null_d = int(B_setup.shape[1] - base_rank)
         if null_d <= 0:
             raise NotImplementedError(_fs_full_rank_base_error(base_spec.bs))
 
@@ -637,13 +688,13 @@ def fit(self, X, feature_names):
         # mgcv uses nat.param(X, S, rank, type=1): eigendecompose R^{-T} S R^{-1}
         # (R from QR of X) and normalise the range space to an identity penalty.
         rp = nat_param_type1(
-            B0,
+            B_setup,
             S0,
             rank=base_rank,
             unit_fnorm=True,
         )
-        X_reparam = rp["X"]  # (n, p0) reparameterised basis
         P_coef = rp["P"]  # (p0, p0) transform: B0 @ P_coef = X_reparam
+        X_reparam = B0 @ P_coef
         r = rp["rank"]  # penalty rank
         D = rp["D"]  # scale^2 * ones(r) after type=1 + unit_fnorm
         null_d = B0.shape[1] - r
@@ -866,7 +917,7 @@ def fit(self, X, feature_names):
             label=self.label,
             fixed=self.fixed,
             by=None,
-            knots=self.knots,
+            knots=self._metric_knots_for_base(),
             xt_rest=base_spec.xt_rest,
             outer_m=self.m,
             mode="sz",
@@ -882,8 +933,8 @@ def fit(self, X, feature_names):
             )
 
         self._base_term = base_term
-        B0, S0, _ = self._base_constructor_fit_matrices()
-        B0 = np.asarray(B0, dtype=np.float64)
+        _B_setup, S0, _ = self._base_constructor_fit_matrices()
+        B0 = np.asarray(self._base_constructor_predict_matrix(X), dtype=np.float64)
         S0 = np.asarray(S0, dtype=np.float64)
 
         level_lists = []

nampy/gam/smooths/smooth_base.py

@@ -604,24 +604,43 @@ def _linked_id_setup_matrix(self, feature_names):
             X_shared[:, feature_names.index(str(name))] = col
         return X_shared
 
-    def _linked_id_marginal_setups(self):
+    def _linked_id_marginal_setups(self, feature_groups=None):
         cols = self._linked_id_pooled_columns()
         setup = self._linked_id_setup()
         if cols is None or setup is None:
             return None
-        return [
-            {
-                "mode": "linked_id",
-                "id": str(setup.get("id")),
-                "pooled_feature_names": [str(name)],
-                "pooled_feature_values": [np.asarray(col, dtype=object).copy()],
-                "n_linked_terms": int(setup.get("n_linked_terms", 0)),
-                "linked_term_labels": list(setup.get("linked_term_labels", [])),
-            }
-            for name, col in zip(
-                setup.get("pooled_feature_names", []), cols, strict=True
+        pooled_names = [str(name) for name in setup.get("pooled_feature_names", [])]
+        if feature_groups is None:
+            groups = [(name,) for name in pooled_names]
+        else:
+            groups = [
+                tuple(group) if isinstance(group, (list, tuple)) else (group,)
+                for group in feature_groups
+            ]
+        column_by_name = {
+            name: np.asarray(column, dtype=object).copy()
+            for name, column in zip(pooled_names, cols, strict=True)
+        }
+        out = []
+        for group in groups:
+            names = [str(name) for name in group]
+            missing = [name for name in names if name not in column_by_name]
+            if missing:
+                raise KeyError(
+                    f"Linked tensor marginal features {missing!r} are absent from "
+                    "the pooled basis setup."
+                )
+            out.append(
+                {
+                    "mode": "linked_id",
+                    "id": str(setup.get("id")),
+                    "pooled_feature_names": names,
+                    "pooled_feature_values": [column_by_name[name] for name in names],
+                    "n_linked_terms": int(setup.get("n_linked_terms", 0)),
+                    "linked_term_labels": list(setup.get("linked_term_labels", [])),
+                }
             )
-        ]
+        return out
 
     @abc.abstractmethod
     def fit(self, X, feature_names):

nampy/gam/smooths/tensor/marginals.py

@@ -9,10 +9,13 @@
 from ..smooth_base import column_as_float
 from ..univariate.bs import DerivativeBSplineTerm1D
 from ..univariate.cr import CubicSplineTerm
+from ..univariate.ds import DuchonSplineTerm
 from ..univariate.ps import PSplineTerm1D
 from ..univariate.tp import ThinPlateSplineTerm
 
-TENSOR_MARGINAL_BASES = frozenset({"bs", "cr", "cs", "cc", "cp", "ps", "tp", "ts"})
+TENSOR_MARGINAL_BASES = frozenset(
+    {"bs", "cr", "cs", "cc", "cp", "ds", "ps", "tp", "ts"}
+)
 
 
 def _as_marginal_features(feature):
@@ -111,6 +114,22 @@ def make_tensor_marginal_term(
             metadata=metadata,
         )
 
+    if basis == "ds":
+        return DuchonSplineTerm(
+            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,
@@ -340,12 +359,11 @@ def build_tensor_product_components(
         )
         shared_setup = getattr(m, "shared_basis_setup", None)
         use_linked_id_predict_path = (
-            len(marginal_indices) == 1
-            and isinstance(shared_setup, dict)
+            isinstance(shared_setup, dict)
             and str(shared_setup.get("mode", "")).lower() == "linked_id"
             and shared_setup.get("pooled_feature_values")
         )
-        if use_linked_id_predict_path:
+        if use_linked_id_predict_path and len(marginal_indices) == 1:
             x_train = np.asarray(
                 shared_setup["pooled_feature_values"][0], dtype=np.float64
             ).ravel()
@@ -399,9 +417,7 @@ def tensor_predict_matrix(
         )
     centered = _normalize_bool_list(centered, len(marginals))
     blocks = []
-    for m, center_i, xp in zip(
-        marginals, centered, np_transforms, strict=True
-    ):
+    for m, center_i, xp in zip(marginals, centered, np_transforms, strict=True):
         blocks.append(
             tensor_marginal_predict_matrix(m, X_new, centered=center_i, np_transform=xp)
         )