feat(gam): add cyclic P-spline smooths

README.md

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

@@ -139,14 +139,15 @@ def instantiate_term(term_like: TermSpec | Any):
         )
 
     if isinstance(smooth_spec, PSplineSmoothSpec):
+        basis = str(smooth_spec.bs).lower()
         if len(features) != 1:
             raise NotImplementedError(
-                "Current runtime only materializes 1D s(..., bs='ps') terms."
+                f"Current runtime only materializes 1D s(..., bs={basis!r}) terms."
             )
         return PSplineTerm1D(
             feature=features[0],
             k=smooth_spec.k,
-            basis="ps",
+            basis=basis,
             m=smooth_spec.m,
             label=label,
             term_id=term_like.term_id,

nampy/gam/smooths/__init__.py

@@ -29,7 +29,7 @@
 ti = InteractionTensorProductSplineTerm
 
 cr = cs = cc = CubicSplineTerm
-ps = PSplineTerm1D
+cp = ps = PSplineTerm1D
 tp = ts = ThinPlateSplineTerm
 fs = FSmoothInteractionTerm
 sz = SZSmoothInteractionTerm
@@ -68,6 +68,7 @@
     "cr",
     "cs",
     "cc",
+    "cp",
     "ps",
     "tp",
     "ts",

nampy/gam/smooths/categorical/fs.py

@@ -106,7 +106,7 @@ def _build_base_smooth_term(
     Build the per-level base smooth used inside fs/sz.
 
     Supported base smooth classes in the current codebase:
-    cr, cs, cc, ps, tp, ts
+    cr, cs, cc, cp, ps, tp, ts
     """
     base_bs = str(base_bs).lower()
     metric_features = list(metric_features)
@@ -123,9 +123,10 @@ def _build_base_smooth_term(
             f"for bs in {{'tp','ts'}}, got base bs={base_bs!r}."
         )
 
-    if xt_rest is not None and base_bs not in {"tp", "ts", "ps"}:
+    if xt_rest is not None and base_bs not in {"tp", "ts", "ps", "cp"}:
         raise NotImplementedError(
-            f"Extra xt options are currently only supported for tp/ts/ps base smooths, "
+            "Extra xt options are currently only supported for tp/ts/ps/cp base "
+            "smooths, "
             f"got xt={xt_rest!r} with base bs={base_bs!r}."
         )
 
@@ -147,15 +148,15 @@ def _build_base_smooth_term(
             metadata=metadata,
         )
 
-    if base_bs == "ps":
+    if base_bs in {"ps", "cp"}:
         ps_m = None if xt_rest is None else xt_rest.get("m", None)
         # For fs/sz, mgcv keeps the outer basis dimension and uses xt mainly to
         # choose the base smoother family / order parameters.
         ps_k = k
         return PSplineTerm1D(
             feature=metric_features[0],
             k=ps_k,
-            basis="ps",
+            basis=base_bs,
             m=ps_m,
             label=label,
             smoothing_id=None,
@@ -191,12 +192,14 @@ def _build_base_smooth_term(
 
     raise NotImplementedError(
         f"Current {mode} implementation supports base bs in "
-        f"{{'cr','cs','cc','ps','tp','ts'}}, got {base_bs!r}."
+        f"{{'cr','cs','cc','cp','ps','tp','ts'}}, got {base_bs!r}."
     )
 
 
 def _penalty_rank_from_base_term(base_term, basis_matrix, penalty_matrix) -> int:
     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)
         # mgcv::smooth.construct.ps.smooth.spec uses rank <- bs.dim - m[2].
         penalty_order = int(base_term.m[1])
         return max(0, int(basis_matrix.shape[1]) - penalty_order)

nampy/gam/smooths/tensor/marginals.py

@@ -11,7 +11,7 @@
 from ..univariate.ps import PSplineTerm1D
 from ..univariate.tp import ThinPlateSplineTerm
 
-TENSOR_MARGINAL_BASES = frozenset({"cr", "cs", "cc", "ps", "tp", "ts"})
+TENSOR_MARGINAL_BASES = frozenset({"cr", "cs", "cc", "cp", "ps", "tp", "ts"})
 
 
 def _as_marginal_features(feature):
@@ -72,11 +72,11 @@ def make_tensor_marginal_term(
             metadata=metadata,
         )
 
-    if basis == "ps":
+    if basis in {"ps", "cp"}:
         if len(marginal_features) != 1:
             raise ValueError(
-                "Tensor marginal basis 'ps' only handles one feature; mgcv coerces "
-                "multivariate ps marginals to tp before construction."
+                f"Tensor marginal basis {basis!r} only handles one feature; mgcv "
+                "coerces multivariate ps/cp marginals to tp before construction."
             )
         return PSplineTerm1D(
             feature=marginal_features[0],

nampy/gam/smooths/univariate/__init__.py

@@ -3,8 +3,8 @@
 from .tp import ThinPlateSplineTerm
 
 cr = cs = cc = CubicSplineTerm
-ps = PSplineTerm1D
+cp = ps = PSplineTerm1D
 tp = ts = ThinPlateSplineTerm
 
 __all__ = ["CubicSplineTerm", "PSplineTerm1D", "ThinPlateSplineTerm"]
-__all__ += ["cr", "cs", "cc", "ps", "tp", "ts"]
+__all__ += ["cr", "cs", "cc", "cp", "ps", "tp", "ts"]

nampy/gam/smooths/univariate/ps.py

@@ -1,11 +1,9 @@
 """
-P-spline smooth term (``bs='ps'``).
+P-spline smooth terms (``bs='ps'`` and cyclic ``bs='cp'``).
 
 Implements the :class:`BaseSmoothTerm` interface for a P-spline: a B-spline
-basis with a discrete difference penalty on adjacent coefficients.  Unlike
-regression splines, P-splines do not require a set of knots to be chosen
-ahead of time; instead, many equally-spaced knots are used and the smoothness
-is controlled entirely by the penalty order and the smoothing parameter.
+basis with a discrete coefficient-difference penalty. Cyclic P-splines use a
+wrapped basis, a circular difference penalty, and periodic newdata mapping.
 """
 
 import numpy as np
@@ -15,7 +13,7 @@
 from ...splines.univariate.ps import (
     build_pspline_term_setup,
     predict_pspline_term,
-    pspline_predict_matrix,
+    predict_pspline_term_derivative,
 )
 from ..registry import register_smooth
 from ..smooth_base import (
@@ -29,6 +27,7 @@
 
 
 @register_smooth("ps")
+@register_smooth("cp")
 class PSplineTerm1D(BaseSmoothTerm):
     term_type = "smooth"
     basis_name = "ps"
@@ -71,22 +70,36 @@ def __init__(
         self.knots = knots
         self.null_penalty_tol = float(null_penalty_tol)
 
+        def normalize_order(value):
+            if value is None or (isinstance(value, float) and np.isnan(value)):
+                return 2
+            numeric = float(value)
+            if not np.isfinite(numeric) or numeric != np.rint(numeric):
+                raise ValueError(
+                    f"For bs={self.basis_name!r}, m entries must be integers or NA."
+                )
+            return int(numeric)
+
         if m is None:
             self.m = (2, 2)
         elif np.isscalar(m):
-            self.m = (int(m), int(m))
+            value = normalize_order(m)
+            self.m = (value, value)
         else:
-            vals = tuple(int(v) for v in m)
+            vals = tuple(normalize_order(v) for v in m)
             if len(vals) == 1:
                 self.m = (vals[0], vals[0])
-            elif len(vals) == 2:
+            elif len(vals) == 2 or (self.basis_name == "cp" and len(vals) > 2):
                 self.m = vals
             else:
-                raise ValueError("For bs='ps', m must have length 1 or 2.")
+                raise ValueError(
+                    f"For bs={self.basis_name!r}, m must have length 1 or 2."
+                )
 
-        if self.basis_name != "ps":
+        if self.basis_name not in {"ps", "cp"}:
             raise NotImplementedError(
-                f"PSplineTerm1D currently supports only basis='ps', got {basis!r}."
+                "PSplineTerm1D supports only basis in {'ps', 'cp'}, "
+                f"got {basis!r}."
             )
         if self.select and self.fixed:
             raise ValueError("select=True and fixed=True are incompatible.")
@@ -131,9 +144,11 @@ def fit(self, X, feature_names):
         else:
             x_setup_values = np.asarray(xj, dtype=np.float64).reshape(-1)
 
-        basis_order, penalty_order = self.m
+        basis_order, penalty_order = self.m[:2]
         if basis_order < 0 or penalty_order < 0:
-            raise ValueError("For bs='ps', m entries must be >= 0.")
+            raise ValueError(
+                f"For bs={self.basis_name!r}, m entries must be >= 0."
+            )
 
         shared_X = self._linked_id_setup_matrix(feature_names)
         if shared_X is not None:
@@ -149,6 +164,7 @@ def fit(self, X, feature_names):
                 bs_dim=self.k,
                 m=self.m,
                 knots=self.knots,
+                basis=self.basis_name,
             )
             setup_base = np.asarray(self._setup.basis_train, dtype=np.float64)
             base = np.asarray(predict_pspline_term(xj, self._setup), dtype=np.float64)
@@ -163,6 +179,7 @@ def fit(self, X, feature_names):
                 bs_dim=self.k,
                 m=self.m,
                 knots=self.knots,
+                basis=self.basis_name,
             )
             point_base = np.asarray(self._setup.basis_train, dtype=np.float64)
             if self._linear_functional:
@@ -309,12 +326,7 @@ def derivative_matrix(self, X_new=None, order=1):
             )
         source = self._X_train if X_new is None else X_new
         xj = column_as_numeric_array(source, self._feature_index)
-        B = pspline_predict_matrix(
-            xj,
-            self._setup.knots,
-            basis_order=self._setup.basis_order,
-            deriv=order,
-        )
+        B = predict_pspline_term_derivative(xj, self._setup, deriv=order)
         return self._apply_constraint_transform_and_by(B, source)
 
     def tensor_marginal_fit_matrices(

nampy/gam/specs/modeling.py

@@ -42,15 +42,15 @@ def make_predictor_specs(model, feature_names, *, knots=None):
                     metadata={},
                 )
             )
-        elif basis == "ps":
+        elif basis in {"ps", "cp"}:
             main_terms.append(
                 TermSpec(
                     kind="smooth",
                     features=(str(name),),
                     by_variable=None,
                     smooth_spec=build_smooth_spec(
                         special="s",
-                        bs="ps",
+                        bs=basis,
                         k=model.k,
                         m=None,
                         sp=None,
@@ -105,7 +105,7 @@ def make_predictor_specs(model, feature_names, *, knots=None):
         else:
             raise NotImplementedError(
                 "Automatic main-effect construction currently supports "
-                "{'cr','cs','cc','ps','tp','ts','re'}, "
+                "{'cr','cs','cc','cp','ps','tp','ts','re'}, "
                 f"got {model.basis!r}."
             )
 

nampy/gam/specs/smooth_build.py

@@ -86,6 +86,7 @@ def _build_s_cc(opts) -> CyclicCubicRegressionSmoothSpec:
 def _build_s_ps(opts) -> PSplineSmoothSpec:
     return PSplineSmoothSpec(
         special="s",
+        bs=str(opts["bs"]).lower(),
         k=opts["k"],
         fx=opts["fx"],
         select=opts["select"],
@@ -182,6 +183,7 @@ def _build_s_sz(opts) -> SumToZeroFactorSmoothSpec:
     "cs": _build_s_cs,
     "cc": _build_s_cc,
     "ps": _build_s_ps,
+    "cp": _build_s_ps,
     "tp": _build_s_tp,
     "ts": _build_s_ts,
     "re": _build_s_re,
@@ -274,7 +276,7 @@ def _is_vector_fx(fx) -> bool:
     return fx is not None and not np.isscalar(fx)
 
 
-_PC_SUPPORTED_S_BASES = {"cc", "cr", "cs", "ps", "tp", "ts"}
+_PC_SUPPORTED_S_BASES = {"cc", "cp", "cr", "cs", "ps", "tp", "ts"}
 
 
 def _dispatch_smooth_spec_from_options(opts) -> SmoothSpec:
@@ -290,7 +292,7 @@ def _dispatch_smooth_spec_from_options(opts) -> SmoothSpec:
             raise NotImplementedError(
                 f"pc= is not supported for s(..., bs={merged['bs']!r}); "
                 "point constraints are only supported for bs in "
-                "{'cc', 'cr', 'cs', 'ps', 'tp', 'ts'}."
+                "{'cc', 'cp', 'cr', 'cs', 'ps', 'tp', 'ts'}."
             )
         return builder(merged)
     if has_pc and special_key not in {"te", "ti"}:

nampy/gam/splines/univariate/__init__.py

@@ -11,7 +11,12 @@
     PSplineBasisSetup,
     bspline_design_matrix,
     build_pspline_term_setup,
+    cyclic_pspline_design,
+    cyclic_pspline_difference_penalty,
+    cyclic_pspline_knots,
+    cyclic_wrap,
     predict_pspline_term,
+    predict_pspline_term_derivative,
     pspline_difference_penalty,
     pspline_knots,
     pspline_predict_matrix,
@@ -23,6 +28,10 @@
     "bspline_design_matrix",
     "cyclic_cubic_bd",
     "cyclic_cubic_predict_matrix",
+    "cyclic_pspline_design",
+    "cyclic_pspline_difference_penalty",
+    "cyclic_pspline_knots",
+    "cyclic_wrap",
     "place_knots_through_values",
     "pspline_difference_penalty",
     "pspline_knots",
@@ -31,6 +40,7 @@
     "PSplineBasisSetup",
     "build_pspline_term_setup",
     "predict_pspline_term",
+    "predict_pspline_term_derivative",
     "build_tprs_term_setup",
     "predict_tprs_term",
 ]