diff --git a/README.md b/README.md index 2b8b49d8..e552e36a 100644 --- a/README.md +++ b/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 | diff --git a/nampy/gam/compiler/factory.py b/nampy/gam/compiler/factory.py index c40cb0fc..c1739629 100644 --- a/nampy/gam/compiler/factory.py +++ b/nampy/gam/compiler/factory.py @@ -17,6 +17,7 @@ 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 ( @@ -24,6 +25,7 @@ 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( diff --git a/nampy/gam/smooths/__init__.py b/nampy/gam/smooths/__init__.py index 5c86567c..44c5b6e4 100644 --- a/nampy/gam/smooths/__init__.py +++ b/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", diff --git a/nampy/gam/smooths/categorical/fs.py b/nampy/gam/smooths/categorical/fs.py index 3bcf2e4b..f674866d 100644 --- a/nampy/gam/smooths/categorical/fs.py +++ b/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 = [] diff --git a/nampy/gam/smooths/smooth_base.py b/nampy/gam/smooths/smooth_base.py index fd78ba7b..1af8e0c3 100644 --- a/nampy/gam/smooths/smooth_base.py +++ b/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): diff --git a/nampy/gam/smooths/tensor/marginals.py b/nampy/gam/smooths/tensor/marginals.py index 04afb8fb..a25906f5 100644 --- a/nampy/gam/smooths/tensor/marginals.py +++ b/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) ) diff --git a/nampy/gam/smooths/tensor/te.py b/nampy/gam/smooths/tensor/te.py index b3d5aad3..3562ee0d 100644 --- a/nampy/gam/smooths/tensor/te.py +++ b/nampy/gam/smooths/tensor/te.py @@ -115,7 +115,7 @@ def __init__( self._by_state = None def fit(self, X, feature_names): - marginal_shared_setups = self._linked_id_marginal_setups() + marginal_shared_setups = self._linked_id_marginal_setups(self.feature) marginals, _, _ = build_tensor_marginal_terms( feature=self.feature, k=self.k, diff --git a/nampy/gam/smooths/tensor/ti.py b/nampy/gam/smooths/tensor/ti.py index cf91f29c..aca048d9 100644 --- a/nampy/gam/smooths/tensor/ti.py +++ b/nampy/gam/smooths/tensor/ti.py @@ -121,7 +121,7 @@ def fit(self, X, feature_names): self._set_by_state(X, feature_names) self._mc = _normalize_mc(self.mc, len(self.feature)) - marginal_shared_setups = self._linked_id_marginal_setups() + marginal_shared_setups = self._linked_id_marginal_setups(self.feature) marginals, _, _ = build_tensor_marginal_terms( feature=self.feature, k=self.k, diff --git a/nampy/gam/smooths/univariate/__init__.py b/nampy/gam/smooths/univariate/__init__.py index 98b81bf8..87b0c1c7 100644 --- a/nampy/gam/smooths/univariate/__init__.py +++ b/nampy/gam/smooths/univariate/__init__.py @@ -1,17 +1,20 @@ from .bs import DerivativeBSplineTerm1D from .cr import CubicSplineTerm +from .ds import DuchonSplineTerm from .ps import PSplineTerm1D from .tp import ThinPlateSplineTerm bs = DerivativeBSplineTerm1D cr = cs = cc = CubicSplineTerm +ds = DuchonSplineTerm cp = ps = PSplineTerm1D tp = ts = ThinPlateSplineTerm __all__ = [ "DerivativeBSplineTerm1D", "CubicSplineTerm", + "DuchonSplineTerm", "PSplineTerm1D", "ThinPlateSplineTerm", ] -__all__ += ["bs", "cr", "cs", "cc", "cp", "ps", "tp", "ts"] +__all__ += ["bs", "cr", "cs", "cc", "cp", "ds", "ps", "tp", "ts"] diff --git a/nampy/gam/smooths/univariate/ds.py b/nampy/gam/smooths/univariate/ds.py new file mode 100644 index 00000000..6bd84fb4 --- /dev/null +++ b/nampy/gam/smooths/univariate/ds.py @@ -0,0 +1,244 @@ +"""Duchon regression-spline smooth term (``bs='ds'``).""" + +from __future__ import annotations + +import numpy as np + +from ...constraints.absorption import ( + fit_single_penalty_with_constraint_policy, + fit_single_penalty_with_setup_basis, +) +from ...penalties.algebra import penalty_rescale_factor, scale_penalty +from ...splines.univariate.ds import ( + build_duchon_spline_setup, + predict_duchon_spline, +) +from ..registry import register_smooth +from ..smooth_base import BaseSmoothTerm, _resolve_feature, columns_as_float_matrix + + +@register_smooth("ds") +class DuchonSplineTerm(BaseSmoothTerm): + term_type = "smooth" + basis_name = "ds" + supports_tensor_marginal = False + + def __init__( + self, + feature, + k=-1, + m=None, + label=None, + term_id=None, + smoothing_id=None, + by=None, + sp=None, + select=False, + fixed=False, + constraint_mode="auto", + pc=None, + knots=None, + xt=None, + null_penalty_tol=1e-10, + metadata=None, + ): + features = list(feature) if not isinstance(feature, (str, int)) else [feature] + super().__init__( + feature=features, + label=label or f"s({', '.join(map(str, features))})", + term_id=term_id, + smoothing_id=smoothing_id, + by=by, + sp=sp, + metadata=metadata, + ) + self.k = int(k) + self.m = m + self.select = bool(select) + self.fixed = bool(fixed) + self.constraint_mode = str(constraint_mode).lower() + self.pc = pc + self.knots = knots + self.xt = xt + self.null_penalty_tol = float(null_penalty_tol) + + if self.constraint_mode not in {"auto", "factor_by", "always", "never"}: + raise ValueError( + "constraint_mode must be one of " + "{'auto', 'factor_by', 'always', 'never'}." + ) + + self._feature_indices = None + self._feature_names = None + self._by_state = None + self._basis_train = None + self._penalties = None + self._setup = None + + @property + def expected_linked_penalty_count(self): + return None if self.select else 1 + + def fit(self, X, feature_names): + feature_indices = [] + resolved_names = [] + for feature in self.feature: + index, name = _resolve_feature(feature, feature_names) + feature_indices.append(index) + resolved_names.append(name) + + values = columns_as_float_matrix(X, feature_indices) + self._set_by_state(X, feature_names) + self._feature_indices = feature_indices + self._feature_names = resolved_names + self._set_resolved_features(resolved_names) + + shared_X = self._linked_id_setup_matrix(feature_names) + if shared_X is not None: + setup_values = columns_as_float_matrix(shared_X, feature_indices) + else: + setup_values = values + self._setup = build_duchon_spline_setup( + setup_values, + k=self.k, + m=self.m, + knots=self.knots, + xt=self.xt, + ) + setup_base = np.asarray(self._setup.basis_train, dtype=np.float64) + base = ( + setup_base + if shared_X is None + else np.asarray( + predict_duchon_spline(values, self._setup), dtype=np.float64 + ) + ) + raw_penalty = np.asarray(self._setup.penalty, dtype=np.float64) + penalty = scale_penalty(setup_base, raw_penalty) + self._set_penalty_rescale_factors( + [penalty_rescale_factor(setup_base, raw_penalty)] + ) + + if self.pc is not None: + constrained_basis, constrained_penalties, transform, _ = ( + self._apply_point_constraint( + base, + [penalty], + self.pc, + feature_names=self._feature_names, + point_basis_fn=lambda points: predict_duchon_spline( + points, self._setup + )[0], + fixed=self.fixed, + ) + ) + self._basis_train = np.asarray(constrained_basis, dtype=np.float64) + self._penalties = constrained_penalties + self._record_constraint_result("pc", transform, absorbed_by="runtime") + return self + + auto_constrain = bool(self._by_state.is_constant) + if shared_X is None: + result = fit_single_penalty_with_constraint_policy( + base, + penalty, + self._by_state, + constraint_mode=self.constraint_mode, + fixed=self.fixed, + auto_constrain_when=auto_constrain, + ) + else: + result = fit_single_penalty_with_setup_basis( + base, + setup_base, + penalty, + self._by_state, + constraint_mode=self.constraint_mode, + fixed=self.fixed, + auto_constrain_when=auto_constrain, + ) + self._basis_train = result.basis_train + self._penalties = result.penalties + self._record_constraint_result( + result.constraint_kind, + result.constraint_transform, + absorbed_by=( + "runtime" if result.constraint_transform is not None else None + ), + ) + return self + + def get_penalty_definitions(self): + self._require_fitted() + if not self.penalties: + return [] + metadata = { + "term_type": self.term_type, + "basis_name": self.basis_name, + "feature": list(self.feature), + "label": self.label, + "by": self.by, + "by_name": self._by_state.feature_name, + "by_is_constant": bool(self._by_state.is_constant), + "constraint_mode": self.constraint_mode, + "constraint_kind": self.constraint_kind, + "pc": self.pc, + "knots": self.knots, + "xt": self.xt, + "m": self.m, + "penalty_order": self._setup.penalty_order, + "shift_order": self._setup.shift_order, + "original_null_space_dim": self._setup.null_space_dim, + "fixed": bool(self.fixed), + } + selection_metadata = {**metadata, "is_selection_penalty": True} + metadata = self._penalty_metadata_with_scale(metadata, penalty_index=0) + return self._build_penalty_block( + self.penalties[0], + rank=int(self._setup.rank), + smooth_metadata=metadata, + selection_metadata=selection_metadata, + ) + + def transform_new(self, X_new): + self._require_fitted() + values = columns_as_float_matrix(X_new, self._feature_indices) + basis = predict_duchon_spline(values, self._setup) + return self._apply_constraint_transform_and_by(basis, X_new) + + def tensor_marginal_fit_matrices( + self, *, centered=False, apply_np=False, x_train=None + ): + del apply_np, x_train + self._require_fitted() + setup_base = np.asarray(self._setup.basis_train, dtype=np.float64) + setup_penalty = np.asarray(self._setup.penalty, dtype=np.float64) + if centered: + if ( + self._linked_id_setup() is not None + and self.constraint_transform is not None + ): + transform = np.asarray(self.constraint_transform, dtype=np.float64) + scaled = scale_penalty(setup_base, setup_penalty) + return ( + np.asarray(setup_base @ transform, dtype=np.float64), + np.asarray(transform.T @ scaled @ transform, dtype=np.float64), + None, + ) + return super().tensor_marginal_fit_matrices(centered=True) + return setup_base, setup_penalty, None + + def tensor_marginal_predict_matrix( + self, X_new, *, centered=False, np_transform=None + ): + if centered: + basis = np.asarray(self.transform_new(X_new), dtype=np.float64) + else: + values = columns_as_float_matrix(X_new, self._feature_indices) + basis = predict_duchon_spline(values, self._setup) + if np_transform is not None: + basis = basis @ np.asarray(np_transform, dtype=np.float64) + return np.asarray(basis, dtype=np.float64) + + +__all__ = ["DuchonSplineTerm"] diff --git a/nampy/gam/specs/__init__.py b/nampy/gam/specs/__init__.py index b4835990..38c68879 100644 --- a/nampy/gam/specs/__init__.py +++ b/nampy/gam/specs/__init__.py @@ -9,6 +9,7 @@ CubicShrinkageSmoothSpec, CyclicCubicRegressionSmoothSpec, DerivativeBSplineSmoothSpec, + DuchonSplineSmoothSpec, FactorSmoothInteractionSpec, PSplineSmoothSpec, RandomEffectSmoothSpec, @@ -31,6 +32,7 @@ "CubicShrinkageSmoothSpec", "CyclicCubicRegressionSmoothSpec", "DerivativeBSplineSmoothSpec", + "DuchonSplineSmoothSpec", "FactorSmoothInteractionSpec", "PSplineSmoothSpec", "RandomEffectSmoothSpec", diff --git a/nampy/gam/specs/build.py b/nampy/gam/specs/build.py index 1523f650..c73b1958 100644 --- a/nampy/gam/specs/build.py +++ b/nampy/gam/specs/build.py @@ -477,7 +477,9 @@ def _fs_by_without_factor_feature_base_spec(smooth_spec: FactorSmoothInteraction if base_bs == "ps": kwargs["m"] = None if xt_rest is None else xt_rest.get("m", None) - elif base_bs in {"tp", "ts"}: + elif base_bs in {"ds", "tp", "ts"}: + if base_bs == "ds": + kwargs["m"] = smooth_spec.m kwargs["xt"] = xt_rest elif xt_rest: raise NotImplementedError( diff --git a/nampy/gam/specs/modeling.py b/nampy/gam/specs/modeling.py index 0a77e08f..b881fd27 100644 --- a/nampy/gam/specs/modeling.py +++ b/nampy/gam/specs/modeling.py @@ -42,7 +42,7 @@ def make_predictor_specs(model, feature_names, *, knots=None): metadata={}, ) ) - elif basis in {"bs", "ps", "cp"}: + elif basis in {"bs", "ds", "ps", "cp"}: main_terms.append( TermSpec( kind="smooth", @@ -105,7 +105,7 @@ def make_predictor_specs(model, feature_names, *, knots=None): else: raise NotImplementedError( "Automatic main-effect construction currently supports " - "{'bs','cr','cs','cc','cp','ps','tp','ts','re'}, " + "{'bs','cr','cs','cc','cp','ds','ps','tp','ts','re'}, " f"got {model.basis!r}." ) diff --git a/nampy/gam/specs/smooth.py b/nampy/gam/specs/smooth.py index 299ac872..6590e60d 100644 --- a/nampy/gam/specs/smooth.py +++ b/nampy/gam/specs/smooth.py @@ -58,6 +58,15 @@ class DerivativeBSplineSmoothSpec(BaseSmoothSpec): pc: Any = None +@dataclass(frozen=True) +class DuchonSplineSmoothSpec(BaseSmoothSpec): + bs: str = "ds" + m: Any = None + xt: Any = None + constraint_mode: str = "auto" + pc: Any = None + + @dataclass(frozen=True) class ShapeConstrainedSmoothSpec(BaseSmoothSpec): """SCAM SCOP-spline specification for a named shape basis code.""" @@ -132,6 +141,7 @@ class TensorInteractionSmoothSpec(BaseSmoothSpec): CubicRegressionSmoothSpec, CyclicCubicRegressionSmoothSpec, DerivativeBSplineSmoothSpec, + DuchonSplineSmoothSpec, CubicShrinkageSmoothSpec, PSplineSmoothSpec, ShapeConstrainedSmoothSpec, diff --git a/nampy/gam/specs/smooth_build.py b/nampy/gam/specs/smooth_build.py index 94a017a2..ef1da39e 100644 --- a/nampy/gam/specs/smooth_build.py +++ b/nampy/gam/specs/smooth_build.py @@ -13,6 +13,7 @@ CubicShrinkageSmoothSpec, CyclicCubicRegressionSmoothSpec, DerivativeBSplineSmoothSpec, + DuchonSplineSmoothSpec, FactorSmoothInteractionSpec, PSplineSmoothSpec, RandomEffectSmoothSpec, @@ -113,6 +114,21 @@ def _build_s_bs(opts) -> DerivativeBSplineSmoothSpec: ) +def _build_s_ds(opts) -> DuchonSplineSmoothSpec: + return DuchonSplineSmoothSpec( + special="s", + k=opts["k"], + fx=opts["fx"], + select=opts["select"], + sp=opts["sp"], + knots=opts["knots"], + m=opts["m"], + xt=opts["xt"], + constraint_mode=opts["constraint_mode"], + pc=opts["pc"], + ) + + def _build_s_shape(opts) -> ShapeConstrainedSmoothSpec: return ShapeConstrainedSmoothSpec( special="s", @@ -202,6 +218,7 @@ def _build_s_sz(opts) -> SumToZeroFactorSmoothSpec: "ps": _build_s_ps, "cp": _build_s_ps, "bs": _build_s_bs, + "ds": _build_s_ds, "tp": _build_s_tp, "ts": _build_s_ts, "re": _build_s_re, @@ -294,7 +311,17 @@ def _is_vector_fx(fx) -> bool: return fx is not None and not np.isscalar(fx) -_PC_SUPPORTED_S_BASES = {"bs", "cc", "cp", "cr", "cs", "ps", "tp", "ts"} +_PC_SUPPORTED_S_BASES = { + "bs", + "cc", + "cp", + "cr", + "cs", + "ds", + "ps", + "tp", + "ts", +} def _dispatch_smooth_spec_from_options(opts) -> SmoothSpec: @@ -310,7 +337,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 " - "{'bs', 'cc', 'cp', 'cr', 'cs', 'ps', 'tp', 'ts'}." + "{'bs', 'cc', 'cp', 'cr', 'cs', 'ds', 'ps', 'tp', 'ts'}." ) return builder(merged) if has_pc and special_key not in {"te", "ti"}: @@ -508,11 +535,11 @@ def _default_k_for_smooth(kind, basis, features, default_k): # mgcv::te()/ti() default k to 5^d per marginal. The current # Python tensor surface supports one feature per marginal, so d = 1. return [5] * len(features) - if str(basis).lower() in {"bs", "tp", "ts"}: - # mgcv/R/smooth.r::s() leaves k = -1; smooth.construct.tp.smooth.spec - # resolves the d-dependent default M + c(8, 27, 100)[min(d, 3)] at - # construction time (mgcv/R/smooth.r:1316-1318). A flat default here - # would wrongly give k = 10 for d > 1. + if str(basis).lower() in {"bs", "ds", "tp", "ts"}: + # mgcv/R/smooth.r::s() leaves k = -1. The basis constructor then + # resolves its dimension-dependent default (TP/TS use M + 8/27/100; + # DS uses M + 10/30/100). A flat default here would be wrong in more + # than one dimension. return -1 return _default_k_for_basis(basis, default_k) diff --git a/nampy/gam/splines/univariate/__init__.py b/nampy/gam/splines/univariate/__init__.py index 75026dc4..1ab6b528 100644 --- a/nampy/gam/splines/univariate/__init__.py +++ b/nampy/gam/splines/univariate/__init__.py @@ -16,6 +16,16 @@ cyclic_cubic_predict_matrix, place_knots_through_values, ) +from .ds import ( + DuchonSplineSetup, + build_duchon_spline_setup, + default_duchon_k, + duchon_kernel, + duchon_null_space_dimension, + duchon_polynomial_basis, + normalize_duchon_orders, + predict_duchon_spline, +) from .ps import ( PSplineBasisSetup, bspline_design_matrix, @@ -40,6 +50,14 @@ "derivative_penalty_root", "normalize_bspline_orders", "predict_derivative_bspline", + "DuchonSplineSetup", + "build_duchon_spline_setup", + "default_duchon_k", + "duchon_kernel", + "duchon_null_space_dimension", + "duchon_polynomial_basis", + "normalize_duchon_orders", + "predict_duchon_spline", "add_full_rank_shrinkage", "bspline_design_matrix", "cyclic_cubic_bd", diff --git a/nampy/gam/splines/univariate/ds.py b/nampy/gam/splines/univariate/ds.py new file mode 100644 index 00000000..738faf92 --- /dev/null +++ b/nampy/gam/splines/univariate/ds.py @@ -0,0 +1,422 @@ +"""Duchon regression-spline primitives for ``mgcv``'s ``bs='ds'``.""" + +from __future__ import annotations + +import math +import warnings +from dataclasses import dataclass + +import numpy as np +from scipy.spatial import distance_matrix + +from ...linalg.qr import r_linpack_qr_no_pivot +from ..basis.tp import tp_T +from .tp import _top_eigensystem + + +def normalize_duchon_orders(m, dimension: int) -> tuple[int, float]: + """Mirror ``smooth.construct.ds.smooth.spec`` normalization of ``m``.""" + dimension = int(dimension) + if dimension < 1: + raise ValueError("Duchon smooths require at least one covariate.") + + if m is None: + raw = [np.nan, np.nan] + elif np.isscalar(m): + raw = [m, np.nan] + else: + raw = list(np.asarray(m, dtype=object).ravel()) + raw = (raw + [np.nan, np.nan])[:2] + + def _missing(value): + if value is None: + return True + if isinstance(value, str) and value.strip().upper() == "NA": + return True + try: + return bool(np.isnan(float(value))) + except (TypeError, ValueError): + return False + + try: + penalty_order = 2 if _missing(raw[0]) else int(np.rint(float(raw[0]))) + shift_order = 0.0 if _missing(raw[1]) else float(np.rint(2 * float(raw[1])) / 2) + except (TypeError, ValueError) as exc: + raise ValueError("For bs='ds', m must contain numeric values or NA.") from exc + + penalty_order = max(1, penalty_order) + if shift_order >= dimension / 2: + shift_order = (dimension - 1) / 2 + warnings.warn("s value reduced", stacklevel=2) + if shift_order <= -dimension / 2: + shift_order = -(dimension - 1) / 2 + warnings.warn("s value increased", stacklevel=2) + if penalty_order + shift_order <= dimension / 2: + shift_order = 0.5 + dimension / 2 - penalty_order + if shift_order >= dimension / 2: + raise ValueError("No suitable s (i.e. m[2]) try increasing m[1]") + warnings.warn( + "s value modified to give continuous function", + stacklevel=2, + ) + return int(penalty_order), float(shift_order) + + +def duchon_null_space_dimension(dimension: int, penalty_order: int) -> int: + """Return the polynomial null-space size ``choose(m + d - 1, d)``.""" + return math.comb(int(penalty_order) + int(dimension) - 1, int(dimension)) + + +def default_duchon_k(dimension: int, null_space_dim: int) -> int: + defaults = (10, 30, 100) + return int(null_space_dim) + defaults[min(int(dimension), len(defaults)) - 1] + + +def duchon_polynomial_basis(x, penalty_order: int): + """Port ``DuchonT`` using mgcv's thin-plate polynomial ordering.""" + values = np.asarray(x, dtype=np.float64) + if values.ndim == 1: + values = values.reshape(-1, 1) + dimension = int(values.shape[1]) + null_space_dim = duchon_null_space_dimension(dimension, penalty_order) + return np.asarray( + tp_T(values, null_space_dim, int(penalty_order), dimension), + dtype=np.float64, + ) + + +def duchon_kernel(x, knots, penalty_order: int, shift_order: float): + """Port ``DuchonE`` including its exponent-dependent sign convention.""" + values = np.asarray(x, dtype=np.float64) + knot_values = np.asarray(knots, dtype=np.float64) + if values.ndim == 1: + values = values.reshape(-1, 1) + if knot_values.ndim == 1: + knot_values = knot_values.reshape(-1, 1) + if values.shape[1] != knot_values.shape[1]: + raise ValueError("Duchon data and knots must have the same dimension.") + + distances = distance_matrix(values, knot_values) + exponent_float = 2 * int(penalty_order) + 2 * float(shift_order) - values.shape[1] + exponent = int(np.rint(exponent_float)) + if not np.isclose(exponent_float, exponent, atol=0.0, rtol=0.0): + raise ValueError( + "Duchon kernel exponent must be integral after m normalization." + ) + + if exponent % 2 == 0: + kernel = np.zeros_like(distances, dtype=np.float64) + nonzero = distances != 0.0 + kernel[nonzero] = distances[nonzero] ** exponent * np.log(distances[nonzero]) + else: + kernel = distances**exponent + sign = 1 - 2 * ((math.floor(exponent / 2) + 1) % 2) + return np.asarray(kernel * sign, dtype=np.float64) + + +def _r_linpack_qty(packed_qr, qraux, values): + """Apply ``t(Q)`` from base R's LINPACK QR representation.""" + packed = np.asarray(packed_qr, dtype=np.float64) + aux = np.asarray(qraux, dtype=np.float64) + out = np.asarray(values, dtype=np.float64).copy() + if out.ndim == 1: + out = out.reshape(-1, 1) + for j in range(min(aux.size, packed.shape[1])): + if aux[j] == 0.0: + continue + reflector = packed[j:, j].copy() + reflector[0] = aux[j] + denominator = float(reflector[0]) + for column in range(out.shape[1]): + step = -float(np.dot(reflector, out[j:, column])) / denominator + out[j:, column] += step * reflector + return np.asarray(out, dtype=np.float64) + + +def _parse_duchon_xt(xt): + max_knots = 2000 + seed = 1 + if xt is None: + return max_knots, seed + if not isinstance(xt, dict): + raise NotImplementedError( + "For bs='ds', xt must be None or a dict with optional keys " + "{'max.knots', 'seed'}." + ) + if xt.get("max.knots") is not None: + max_knots = int(xt["max.knots"]) + if xt.get("seed") is not None: + seed = int(xt["seed"]) + if max_knots < 1: + raise ValueError("For bs='ds', xt['max.knots'] must be positive.") + return max_knots, seed + + +def _normalize_duchon_knots(knots, dimension: int): + """Collect per-coordinate knot vectors like the upstream constructor.""" + if knots is None: + return None + + dimension = int(dimension) + if isinstance(knots, (list, tuple)): + if dimension == 1 and (not knots or np.isscalar(knots[0])): + return np.asarray(knots, dtype=np.float64).reshape(-1, 1) + if len(knots) != dimension or any(value is None for value in knots): + return None + columns = [np.asarray(value, dtype=np.float64).ravel() for value in knots] + if any(column.size != columns[0].size for column in columns[1:]): + raise ValueError( + "components of knots relating to a single smooth must be of same length" + ) + return np.asarray(np.column_stack(columns), dtype=np.float64) + + values = np.asarray(knots, dtype=np.float64) + if values.ndim == 1: + if dimension != 1: + return None + return values.reshape(-1, 1) + if values.ndim != 2 or values.shape[1] != dimension: + return None + return np.asarray(values, dtype=np.float64) + + +def _duchon_unique_rows(values): + """Mirror ``uniquecombs(x, ordered=TRUE)`` for numeric matrices.""" + matrix = np.asarray(values, dtype=np.float64) + if matrix.ndim == 1: + matrix = matrix.reshape(-1, 1) + if matrix.shape[1] == 1: + return np.unique(matrix[:, 0]).reshape(-1, 1) + + first_by_key = {} + for row in matrix: + key = "*".join(format(float(value), ".15g") for value in row) + first_by_key.setdefault(key, np.asarray(row, dtype=np.float64).copy()) + return np.vstack([first_by_key[key] for key in sorted(first_by_key)]) + + +class _RMersenneTwister: + """Minimal R-compatible MT19937 stream used by ``temp.seed``/``sample``.""" + + def __init__(self, seed): + state = int(seed) & 0xFFFFFFFF + for _ in range(50): + state = (69069 * state + 1) & 0xFFFFFFFF + state = (69069 * state + 1) & 0xFFFFFFFF # R stores 624 in this slot. + self._state = [] + for _ in range(624): + state = (69069 * state + 1) & 0xFFFFFFFF + self._state.append(state) + self._index = 624 + + def _twist(self): + for index in range(624): + word = (self._state[index] & 0x80000000) | ( + self._state[(index + 1) % 624] & 0x7FFFFFFF + ) + self._state[index] = ( + self._state[(index + 397) % 624] + ^ (word >> 1) + ^ (0x9908B0DF if word & 1 else 0) + ) + self._index = 0 + + def uniform(self): + if self._index >= 624: + self._twist() + word = self._state[self._index] + self._index += 1 + word ^= word >> 11 + word ^= (word << 7) & 0x9D2C5680 + word ^= (word << 15) & 0xEFC60000 + word ^= word >> 18 + return float(word & 0xFFFFFFFF) / float(2**32) + + def uniform_index(self, size): + size = int(size) + bits = int(math.ceil(math.log2(size))) + mask = (1 << bits) - 1 + while True: + value = 0 + for _ in range(0, bits + 1, 16): + value = 65536 * value + math.floor(self.uniform() * 65536) + value &= mask + if value < size: + return int(value) + + +def _r_sample_without_replacement(size, sample_size, seed): + """Mirror modern R ``sample.int(..., replace=FALSE)`` rejection sampling.""" + available = list(range(int(size))) + remaining = int(size) + rng = _RMersenneTwister(seed) + selected = [] + for _ in range(int(sample_size)): + index = rng.uniform_index(remaining) + selected.append(available[index]) + remaining -= 1 + available[index] = available[remaining] + return np.asarray(selected, dtype=np.intp) + + +def _duchon_setup_locations(values, shift, knots, *, max_knots, seed): + values = np.asarray(values, dtype=np.float64) + if knots is not None: + return np.asarray(knots, dtype=np.float64), False + unique = _duchon_unique_rows(values) - np.asarray(shift, dtype=np.float64)[None, :] + if values.shape[0] > int(max_knots) and unique.shape[0] > int(max_knots): + indices = _r_sample_without_replacement( + unique.shape[0], + int(max_knots), + int(seed), + ) + return np.asarray(unique[indices, :], dtype=np.float64), True + return np.asarray(unique, dtype=np.float64), False + + +@dataclass +class DuchonSplineSetup: + shift: np.ndarray + knots: np.ndarray + UZ: np.ndarray + penalty_order: int + shift_order: float + null_space_dim: int + rank: int + bs_dim: int + basis_train: np.ndarray + penalty: np.ndarray + used_supplied_knots: bool + used_subsampling: bool + + +def build_duchon_spline_setup(X, *, k=-1, m=None, knots=None, xt=None): + """Port ``smooth.construct.ds.smooth.spec`` and retain prediction state.""" + values = np.asarray(X, dtype=np.float64) + if values.ndim == 1: + values = values.reshape(-1, 1) + if values.ndim != 2 or values.shape[1] < 1: + raise ValueError("Duchon smooth data must be a non-empty numeric matrix.") + + n_obs, dimension = values.shape + penalty_order, shift_order = normalize_duchon_orders(m, dimension) + null_space_dim = duchon_null_space_dimension(dimension, penalty_order) + bs_dim = int(k) + if bs_dim < 0: + bs_dim = default_duchon_k(dimension, null_space_dim) + if bs_dim < null_space_dim + 1: + bs_dim = null_space_dim + 1 + warnings.warn("basis dimension reset to minimum possible", stacklevel=2) + + unique = _duchon_unique_rows(values) + if unique.shape[0] < bs_dim: + raise ValueError( + "A term has fewer unique covariate combinations than specified " + "maximum degrees of freedom" + ) + + shift = np.mean(values, axis=0) + supplied = _normalize_duchon_knots(knots, dimension) + if supplied is not None and supplied.shape[0] == 0: + supplied = None + if supplied is not None and supplied.shape[0] > n_obs: + warnings.warn( + "more knots than data in a ds term: knots ignored.", + stacklevel=2, + ) + supplied = None + if supplied is not None: + supplied = supplied - shift[None, :] + + max_knots, seed = _parse_duchon_xt(xt) + setup_knots, used_subsampling = _duchon_setup_locations( + values, + shift, + supplied, + max_knots=max_knots, + seed=seed, + ) + n_knots = int(setup_knots.shape[0]) + if n_knots < bs_dim: + raise ValueError( + "Duchon spline requires at least as many knot locations as basis " + "coefficients." + ) + + E = duchon_kernel(setup_knots, setup_knots, penalty_order, shift_order) + T = duchon_polynomial_basis(setup_knots, penalty_order) + if bs_dim < n_knots: + eigenvalues, eigenvectors = _top_eigensystem( + E, + bs_dim, + tolerance_exponent=0.5, + ) + diagonal_penalty = np.diag(eigenvalues) + constraint = (T.T @ eigenvectors).T + else: + eigenvectors = np.eye(bs_dim, dtype=np.float64) + diagonal_penalty = np.asarray(E, dtype=np.float64) + constraint = np.asarray(T, dtype=np.float64) + + packed_qr, qraux = r_linpack_qr_no_pivot(constraint) + first = _r_linpack_qty(packed_qr, qraux, diagonal_penalty) + reduced = _r_linpack_qty( + packed_qr, + qraux, + first[null_space_dim:, :].T, + )[null_space_dim:, :] + penalty = np.zeros((bs_dim, bs_dim), dtype=np.float64) + penalty[: bs_dim - null_space_dim, : bs_dim - null_space_dim] = reduced + + UZ = _r_linpack_qty( + packed_qr, + qraux, + eigenvectors.T, + )[null_space_dim:, :].T + setup = DuchonSplineSetup( + shift=np.asarray(shift, dtype=np.float64), + knots=np.asarray(setup_knots, dtype=np.float64), + UZ=np.asarray(UZ, dtype=np.float64), + penalty_order=int(penalty_order), + shift_order=float(shift_order), + null_space_dim=int(null_space_dim), + rank=int(bs_dim - null_space_dim), + bs_dim=int(bs_dim), + basis_train=np.zeros((n_obs, bs_dim), dtype=np.float64), + penalty=np.asarray(penalty, dtype=np.float64), + used_supplied_knots=bool(supplied is not None), + used_subsampling=bool(used_subsampling), + ) + setup.basis_train = predict_duchon_spline(values, setup) + return setup + + +def predict_duchon_spline(X_new, setup: DuchonSplineSetup): + """Port ``Predict.matrix.duchon.spline``.""" + values = np.asarray(X_new, dtype=np.float64) + if values.ndim == 1: + values = values.reshape(-1, 1) + shifted = values - setup.shift[None, :] + radial = duchon_kernel( + shifted, + setup.knots, + setup.penalty_order, + setup.shift_order, + ) + polynomial = duchon_polynomial_basis(shifted, setup.penalty_order) + return np.asarray( + np.column_stack([radial @ setup.UZ, polynomial]), + dtype=np.float64, + ) + + +__all__ = [ + "DuchonSplineSetup", + "build_duchon_spline_setup", + "default_duchon_k", + "duchon_kernel", + "duchon_null_space_dimension", + "duchon_polynomial_basis", + "normalize_duchon_orders", + "predict_duchon_spline", +] diff --git a/nampy/gam/splines/univariate/tp.py b/nampy/gam/splines/univariate/tp.py index 95052c56..04e752ff 100644 --- a/nampy/gam/splines/univariate/tp.py +++ b/nampy/gam/splines/univariate/tp.py @@ -285,7 +285,7 @@ def choose_tprs_setup_locations(X_shifted, knots=None, max_knots=2000, seed=1): return np.asarray(Xu, dtype=np.float64) -def _top_eigensystem(E, k): +def _top_eigensystem(E, k, *, tolerance_exponent=0.7): """ mgcv-compatible top-k eigensystem for a symmetric matrix E. @@ -304,7 +304,7 @@ def _top_eigensystem(E, k): if k > n: raise ValueError(f"k must be <= matrix dimension, got k={k}, n={n}.") - tol = float(np.finfo(np.float64).eps ** 0.7) + tol = float(np.finfo(np.float64).eps ** float(tolerance_exponent)) # mgcv/src/mat.c::Rlanczos checks convergence every # min(max((m + lm) / 2, 10), floor(n / 10)) steps. For the tp/ts setup # here lm is always zero, so mirror the same cadence exactly. diff --git a/pyproject.toml b/pyproject.toml index 628bc0f9..de0443e1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -207,6 +207,7 @@ markers = [ "smooth_cc: tests covering cyclic cubic smooths", "smooth_bs: tests covering integrated-derivative B-spline smooths", "smooth_cp: tests covering cyclic P-spline smooths", + "smooth_ds: tests covering Duchon regression spline smooths", "smooth_ps: tests covering P-spline smooths", "smooth_tp: tests covering thin plate smooths", "smooth_ts: tests covering shrinkage thin plate smooths", diff --git a/tests/SUBSYSTEM_COVERAGE.md b/tests/SUBSYSTEM_COVERAGE.md index 67225661..d5f3c782 100644 --- a/tests/SUBSYSTEM_COVERAGE.md +++ b/tests/SUBSYSTEM_COVERAGE.md @@ -15,7 +15,7 @@ local development notes rather than duplicated here. | Subsystem | Primary owner(s) | Primary tests | Notes | | --- | --- | --- | --- | | Formula/spec parsing | `nampy/gam/formula/`, `nampy/gam/specs/` | `tests/parity/test_mgcv_formula_parse_parity.py` | Direct formula parity vs `mgcv`. | -| Smooth constructors / raw basis owners | `nampy/gam/smooths/`, `nampy/gam/splines/` | `tests/smooths/test_mgcv_raw_constructor_parity.py`, `tests/smooths/test_mgcv_smoothcon_parity.py`, `tests/parity/test_gam_spec_build_owner_contracts.py`, `tests/parity/test_mgcv_cp_combinations_parity.py`, `tests/parity/test_mgcv_bs_combinations_parity.py` | Basis, penalty, constructor, and smoothCon surfaces, including cyclic P-splines (`cp`) and integrated-derivative B-splines (`bs`) across prediction, multi-penalty, tensor/factor-smooth combinations, and the upstream tensor-`m` wrong-length warning and zero fallback. | +| Smooth constructors / raw basis owners | `nampy/gam/smooths/`, `nampy/gam/splines/` | `tests/smooths/test_mgcv_raw_constructor_parity.py`, `tests/smooths/test_mgcv_smoothcon_parity.py`, `tests/parity/test_gam_spec_build_owner_contracts.py`, `tests/parity/test_mgcv_cp_combinations_parity.py`, `tests/parity/test_mgcv_bs_combinations_parity.py`, `tests/parity/test_mgcv_ds_combinations_parity.py` | Basis, penalty, constructor, and smoothCon surfaces, including cyclic P-splines (`cp`), integrated-derivative B-splines (`bs`), and multivariate Duchon splines (`ds`) across prediction, selection, linked bases, tensor/factor-smooth combinations, and the upstream tensor-`m` wrong-length warning and zero fallback. | | `pc=` / linked `id=` routing | smooth metadata + linked basis owners | `tests/smooths/test_mgcv_pc_id_parity.py` | Localizes shared-smoothing and point-constraint issues, including `te`/`ti` point constraints. | | Design / pre-fit assembly | `nampy/gam/compiler/`, `nampy/gam/fit/penalized_system.py` | `tests/optimization/test_mgcv_gam_setup_assembly_parity.py`, `tests/optimization/test_mgcv_preoptimization_blocks_parity.py`, `tests/optimization/test_mgcv_preoptimization_reparam_parity.py` | Setup, blocks, and reparameterization parity, including one global shared-component block with overlapping linear-predictor indices. | | Term wrapping / by-variable / offset routing | predictor wrapping + compiled term owners | `tests/optimization/test_gam_term_wrapping_owner_contracts.py` | Localizes wrapped predictor blocks, offset routing, and general-family block ownership before broader prediction parity. | @@ -47,7 +47,7 @@ closes the former seven-stage under-tested backlog. | Pipeline stage | Combination coverage now owned | | --- | --- | | 1. Formula parsing and canonical specs | Formula lists/shared components, transformed covariates, supported numeric/factor interactions, factor-by smooths, formula and fit offsets, intercept policies, weights, knots, `min_sp`, `drop_intercept`, and tensor-`m` warning/fallback behavior; the combined interaction recipe is rebuilt on newdata and compared with a committed `mgcv` reference fixture. | -| 2. Runtime terms and low-level bases | Train/newdata pairing for every supported basis (`cr/cs/cc/ps/tp/ts/re/fs/sz/te/ti`), including univariate and structured boundary response/SE parity. Row-permutation contracts cover linked and identified-FS terms, PS/TP/TS bases, `te`/`ti`, SZ, and factor-by smooths. Non-unique constructor representations use column-space projectors and penalized response operators instead of arbitrary coefficient orientation. | +| 2. Runtime terms and low-level bases | Train/newdata pairing for every supported basis (`bs/cr/cs/cc/cp/ds/ps/tp/ts/re/fs/sz/te/ti`), including univariate, multivariate Duchon, and structured boundary response/SE parity. Row-permutation contracts cover linked and identified-FS terms, PS/TP/TS bases, `te`/`ti`, SZ, and factor-by smooths. Non-unique constructor representations use column-space projectors and penalized response operators instead of arbitrary coefficient orientation. | | 3. Constructed/wrapped terms | Numeric-by plus linked `id=`, factor-by plus linked `id=`, tensor-by, and mixed fixed/free/select penalty ownership, numerical coefficient-map composition, and fitted response/SE parity. | | 4. Predictor/model compilation | Multi-predictor layouts with unequal intercept and offset policies, overlapping shared-component coefficient indices, and three-term linked, fixed/free, select, and rank-deficient assembly; supported `gaulss` and `gammals` layouts are also fitted and compared with `mgcv`. | | 5. Side conditions and identifiability | Repeated, nested, reverse-formula-order, tensor/main-effect, three-way, both identified near-rank regimes, zero-width, no-intercept, ordered factor-by, linked, general-family, two-predictor, SZ, and exempt random/factor smooth cases. Deletion rank and behavior replace raw pivot-column identity where QR/eigen choices are non-unique. | diff --git a/tests/TAXONOMY.md b/tests/TAXONOMY.md index 0c4e6aec..d5a6a46c 100644 --- a/tests/TAXONOMY.md +++ b/tests/TAXONOMY.md @@ -12,7 +12,7 @@ The GAM test suite is intentionally overlapping. The goal is fast subset runs an - `tests/`: shared helpers, marker inference, taxonomy registry, static reference fixtures, and parity-generation R scripts ## Taxonomy Axes -- `smooth_`: `bs`, `cr`, `cs`, `cc`, `cp`, `ps`, `tp`, `ts`, `te`, `ti`, `fs`, `sz`, `re` +- `smooth_`: `bs`, `cr`, `cs`, `cc`, `cp`, `ds`, `ps`, `tp`, `ts`, `te`, `ti`, `fs`, `sz`, `re` - `family_`: `gaussian`, `binomial`, `poisson`, `gamma`, `negbin`, `gaulss`, `gammals`, `general` - `method_`: `fixed`, `reml`, `ml`, `laml`, `gcv`, `ubre` - `link_`: `identity`, `log`, `inverse`, `logit`, `probit`, `cloglog`, `cauchit`, `sqrt` diff --git a/tests/_taxonomy_registry.py b/tests/_taxonomy_registry.py index d949bb36..ebe8d7b8 100644 --- a/tests/_taxonomy_registry.py +++ b/tests/_taxonomy_registry.py @@ -19,6 +19,7 @@ def _leaf(leaf_id: str, *nodeid_parts: str) -> LeafCoverageExpectation: "cs": "smooth_cs", "cc": "smooth_cc", "cp": "smooth_cp", + "ds": "smooth_ds", "ps": "smooth_ps", "tp": "smooth_tp", "ts": "smooth_ts", @@ -84,6 +85,7 @@ def _leaf(leaf_id: str, *nodeid_parts: str) -> LeafCoverageExpectation: "test_mgcv_output_parity.py": {"surface_output"}, "test_mgcv_smoothcon_parity.py": {"surface_smoothcon"}, "test_mgcv_raw_constructor_parity.py": {"surface_smoothcon"}, + "test_mgcv_ds_combinations_parity.py": {"surface_snapshot"}, "test_mgcv_score_hist_trace_parity.py": {"surface_trace"}, "test_mgcv_optimization_lifecycle_parity.py": {"surface_trace"}, "test_mgcv_linked_id_trace_parity.py": {"surface_trace"}, @@ -118,6 +120,7 @@ def _leaf(leaf_id: str, *nodeid_parts: str) -> LeafCoverageExpectation: "test_mgcv_pc_id_parity.py", "test_mgcv_known_gaps.py", "test_mgcv_output_parity.py", + "test_mgcv_ds_combinations_parity.py", "test_mgcv_score_hist_trace_parity.py", "test_mgcv_linked_id_trace_parity.py", "test_mgcv_score_gamma_parity.py", @@ -135,6 +138,11 @@ def _leaf(leaf_id: str, *nodeid_parts: str) -> LeafCoverageExpectation: "tests/smooths/test_mgcv_smoothcon_parity.py", "tests/smooths/test_mgcv_pc_id_parity.py", ), + "smooth_ds": ( + "tests/parity/test_mgcv_ds_combinations_parity.py", + "tests/smooths/test_mgcv_raw_constructor_parity.py", + "tests/smooths/test_mgcv_smoothcon_parity.py", + ), "smooth_cs": ( "tests/smooths/test_mgcv_pc_id_parity.py", "tests/smooths/test_mgcv_raw_constructor_parity.py", diff --git a/tests/mgcv_invariant_policy.py b/tests/mgcv_invariant_policy.py index 88eae515..581a3c89 100644 --- a/tests/mgcv_invariant_policy.py +++ b/tests/mgcv_invariant_policy.py @@ -74,7 +74,9 @@ def gam_side_uses_invariant_transform(class_name: str) -> bool: return class_name in _GAM_SIDE_INVARIANT_CLASS_NAMES -def final_fit_uses_exact_orientation_parity(model, *, skip_coef_comparison: bool) -> bool: +def final_fit_uses_exact_orientation_parity( + model, *, skip_coef_comparison: bool +) -> bool: """Return whether compiled terms have uniquely identified coefficient bases.""" if skip_coef_comparison: return False @@ -230,7 +232,15 @@ def _copy_raw_value(value): def _normalized_penalties(value): if isinstance(value, dict): - values = list(value.values()) + keys = list(value) + if "S" in value: + keys = ["S"] + sorted( + (key for key in keys if key != "S"), + key=lambda key: ( + (0, int(key)) if str(key).isdigit() else (1, str(key)) + ), + ) + values = [value[key] for key in keys] else: values = list(value) return [np.asarray(v, dtype=np.float64) for v in values] @@ -288,6 +298,17 @@ def _canonicalize_tprs_raw_state(state): return state +def _canonicalize_duchon_raw_state(state): + extra = state["extra"] + extra.pop("used_supplied_knots", False) + extra.pop("used_subsampling", False) + extra.pop("pure_knot", False) + state["S"] = [penalty_spectrum(S) for S in state["S"]] + state["X"] = matrix_self_gram(state["X"]) + extra["UZ"] = stable_column_space_projector(extra["UZ"]) + return state + + def _canonicalize_cs_raw_state(state): state["S"] = [penalty_spectrum(S) for S in state["S"]] return state @@ -342,6 +363,8 @@ def canonicalize_raw_representation_state(state: dict[str, Any]) -> dict[str, An return _canonicalize_cs_raw_state(state) if class_name in {"tprs.smooth", "ts.smooth"}: return _canonicalize_tprs_raw_state(state) + if class_name == "duchon.spline": + return _canonicalize_duchon_raw_state(state) if class_name == "fs.interaction": return _canonicalize_fs_raw_state(state) if class_name == "sz.interaction": diff --git a/tests/mgcv_parity_utils.py b/tests/mgcv_parity_utils.py index fbd3c993..d2e79681 100644 --- a/tests/mgcv_parity_utils.py +++ b/tests/mgcv_parity_utils.py @@ -1625,6 +1625,12 @@ def _run_mgcv_raw_constructor( shift = pack_vector(sm$shift, "numeric"), drop_null = isTRUE(sm$drop.null != 0) ), + "duchon.spline" = list( + knt = pack_matrix(sm$knt), + UZ = pack_matrix(sm$UZ), + shift = pack_vector(sm$shift, "numeric"), + p_order = pack_vector(sm$p.order, "numeric") + ), "random.effect" = list( C = pack_constraint(sm$C), random = isTRUE(sm$random), diff --git a/tests/parity/test_mgcv_ds_combinations_parity.py b/tests/parity/test_mgcv_ds_combinations_parity.py new file mode 100644 index 00000000..d59c117b --- /dev/null +++ b/tests/parity/test_mgcv_ds_combinations_parity.py @@ -0,0 +1,267 @@ +"""Integrated parity coverage for Duchon regression splines (``bs='ds'``).""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +import pytest + +from nampy.gam import GAM +from nampy.gam.splines.univariate.ds import build_duchon_spline_setup +from tests.mgcv_parity_utils import ( + _assert_basic_mgcv_parity, + _fit_nampy_model, + _fit_nampy_snapshot, + _run_mgcv_snapshot, +) + + +def _ds_data(seed=271, n=180): + rng = np.random.default_rng(seed) + x0 = rng.uniform(-2.0, 2.0, size=n) + x1 = rng.uniform(-1.5, 1.5, size=n) + x2 = rng.uniform(-1.8, 1.8, size=n) + x3 = rng.uniform(-1.2, 1.7, size=n) + z = 0.8 + rng.uniform(-0.4, 0.7, size=n) + f = np.asarray(["a", "b", "c"], dtype=object)[np.arange(n) % 3] + f1 = np.asarray(["u", "v"], dtype=object)[np.arange(n) % 2] + y = ( + 0.2 + + z * np.sin(1.2 * x0) + + 0.3 * x1**2 + - 0.25 * np.cos(x2) + + 0.15 * x3 + + 0.2 * (f == "b") + - 0.15 * (f1 == "v") + + rng.normal(scale=0.12, size=n) + ) + return pd.DataFrame( + { + "y": y, + "x0": x0, + "x1": x1, + "x2": x2, + "x3": x3, + "z": z, + "f": f, + "f1": f1, + } + ) + + +def _assert_snapshot_fit(actual, expected, *, atol=2e-7): + for key in ("response", "link"): + np.testing.assert_allclose( + actual["predictions"][key], + expected["predictions"][key], + atol=atol, + rtol=atol, + ) + np.testing.assert_allclose( + actual["fit"]["edf_total"], + expected["fit"]["edf_total"], + atol=atol, + rtol=atol, + ) + + +def test_ds_numeric_by_select_true_matches_mgcv(): + data = _ds_data(seed=272) + formula = 'y ~ s(x0, x1, by=z, bs="ds", k=10, m=[1,.5])' + actual = _fit_nampy_snapshot(data, formula, "gaussian", "REML", select=True) + expected = _run_mgcv_snapshot(data, formula, "gaussian", "REML", select=True) + assert len(actual["fit"]["smoothing_params"]) == 2 + _assert_basic_mgcv_parity( + actual, + expected, + pred_atol=4e-7, + pred_rtol=4e-7, + sp_log_atol=5e-6, + criterion_atol=2e-7, + ) + + +def test_ds_factor_by_fixed_sp_matches_mgcv(): + data = _ds_data(seed=273) + formula = 'y ~ s(x0, x1, by=f, bs="ds", k=10, m=[1,.5], sp=.7)' + actual = _fit_nampy_snapshot(data, formula, "gaussian", "fixed") + expected = _run_mgcv_snapshot(data, formula, "gaussian", "REML") + _assert_snapshot_fit(actual, expected) + + +def test_ds_linked_multivariate_terms_pool_basis_and_share_sp(): + data = _ds_data(seed=274) + formula = ( + 'y ~ s(x0, x1, bs="ds", k=10, m=[1,.5], id="duchon", sp=.7)' + ' + s(x2, x3, bs="ds", k=10, m=[1,.5], id="duchon")' + ) + model = _fit_nampy_model(data, formula, "gaussian", "fixed") + expected = _run_mgcv_snapshot(data, formula, "gaussian", "REML") + assert len(model.smoothing_params) == 1 + runtimes = [ + term.predict_fn.__self__ + for term in model.gam_result_.compiled_model.compiled_terms + if term.term_type == "smooth" + ] + np.testing.assert_allclose(runtimes[0]._setup.knots, runtimes[1]._setup.knots) + actual = model.parity_snapshot(X=data, include_covariances=True) + _assert_snapshot_fit(actual, expected, atol=5e-7) + + +def test_ds_point_constraint_and_fixed_basis_match_mgcv(): + data = _ds_data(seed=275) + formula = 'y ~ s(x0, x1, bs="ds", k=10, m=[1,.5], pc=[.2,-.3], sp=.8)' + actual = _fit_nampy_snapshot(data, formula, "gaussian", "fixed") + expected = _run_mgcv_snapshot(data, formula, "gaussian", "REML") + _assert_snapshot_fit(actual, expected) + + fixed_formula = 'y ~ s(x2, x3, bs="ds", k=10, m=[1,.5], fx=True)' + model = _fit_nampy_model(data, fixed_formula, "gaussian", "fixed") + assert model.gam_result_.compiled_model.compiled_penalties == () + _assert_snapshot_fit( + model.parity_snapshot(X=data, include_covariances=True), + _run_mgcv_snapshot(data, fixed_formula, "gaussian", "REML"), + ) + selected = _fit_nampy_model( + data, fixed_formula, "gaussian", "fixed", select=True + ) + assert selected.gam_result_.compiled_model.compiled_penalties == () + _assert_snapshot_fit( + selected.parity_snapshot(X=data, include_covariances=True), + _run_mgcv_snapshot( + data, fixed_formula, "gaussian", "REML", select=True + ), + ) + + +@pytest.mark.parametrize( + "special", + ["te", "ti"], +) +def test_ds_multivariate_tensor_margin_fixed_sp_matches_mgcv(special): + data = _ds_data(seed=276, n=150) + formula = ( + f'y ~ {special}(x0, x1, x2, d=[2,1], bs=["ds","cr"], ' + "k=[10,5], m=[[1,.5], None], sp=[.6,.8])" + ) + actual = _fit_nampy_snapshot(data, formula, "gaussian", "fixed") + expected = _run_mgcv_snapshot(data, formula, "gaussian", "REML") + _assert_snapshot_fit(actual, expected, atol=8e-7) + + +def test_ds_linked_multivariate_tensor_margins_share_two_sp(): + data = _ds_data(seed=280, n=140) + data["w0"] = 0.7 * data["x0"] - 0.2 * data["x3"] + data["w1"] = 0.6 * data["x1"] + 0.3 * data["x2"] + formula = ( + 'y ~ te(x0, x1, x2, d=[2,1], bs=["ds","cr"], k=[10,5], ' + 'm=[[1,.5], None], id="tensor_ds", sp=[.6,.8])' + ' + te(x3, w0, w1, d=[2,1], bs=["ds","cr"], k=[10,5], ' + 'm=[[1,.5], None], id="tensor_ds")' + ) + model = _fit_nampy_model(data, formula, "gaussian", "fixed") + assert len(model.smoothing_params) == 2 + _assert_snapshot_fit( + model.parity_snapshot(X=data, include_covariances=True), + _run_mgcv_snapshot(data, formula, "gaussian", "REML"), + atol=1e-6, + ) + + +@pytest.mark.parametrize( + "formula", + [ + 'y ~ s(f, x0, x1, bs="fs", xt="ds", k=10, m=[1,.5], sp=[.7,.9])', + 'y ~ s(f, f1, x0, x1, bs="sz", xt="ds", k=10, m=[1,.5], id="shared", sp=.7)', + ], + ids=["fs", "sz"], +) +def test_ds_multivariate_factor_smooth_base_fixed_sp_matches_mgcv(formula): + data = _ds_data(seed=277, n=150) + actual = _fit_nampy_snapshot(data, formula, "gaussian", "fixed") + expected = _run_mgcv_snapshot(data, formula, "gaussian", "REML") + _assert_snapshot_fit(actual, expected, atol=1e-6) + + +def test_ds_identified_fs_base_uses_dynamic_penalty_vector(): + data = _ds_data(seed=281, n=140) + formula = ( + 'y ~ s(f, x0, x1, bs="fs", xt="ds", k=10, m=[1,.5], ' + 'id="fs_ds", sp=[.7,.9])' + ) + model = _fit_nampy_model(data, formula, "gaussian", "fixed") + assert len(model.smoothing_params) == 2 + _assert_snapshot_fit( + model.parity_snapshot(X=data, include_covariances=True), + _run_mgcv_snapshot(data, formula, "gaussian", "REML"), + atol=1e-6, + ) + + +def test_ds_array_api_persistence_and_blocked_extrapolation(tmp_path): + data = _ds_data(seed=278, n=130) + features = data[["x0", "x1"]] + model = GAM( + family="gaussian", + basis="ds", + k=10, + optimize_smoothing=False, + smoothing_params=[0.5, 0.8], + ).fit(X=features, y=data["y"].to_numpy(dtype=np.float64)) + newdata = pd.DataFrame({"x0": [-4.0, 0.0, 4.0], "x1": [-3.0, 0.5, 3.5]}) + expected = model.predict(newdata, type="link", block_size=1) + path = tmp_path / "ds.pkl" + model.save_model(path) + restored = GAM.load_model(path) + np.testing.assert_allclose( + restored.predict(newdata, type="link", block_size=1), expected + ) + + +def test_ds_order_normalization_warnings_and_validation(): + x = np.linspace(-1.0, 1.0, 30) + X = np.column_stack([x, np.sin(x)]) + with pytest.warns(UserWarning, match="s value reduced"): + setup = build_duchon_spline_setup(X, k=8, m=[1, 3]) + assert setup.shift_order == 0.5 + with pytest.warns(UserWarning) as caught: + build_duchon_spline_setup(X, k=8, m=[1, -3]) + messages = [str(item.message) for item in caught] + assert "s value increased" in messages + assert "s value modified to give continuous function" in messages + with pytest.warns(UserWarning, match="s value modified"): + build_duchon_spline_setup(X, k=8, m=[1, 0]) + with pytest.warns(UserWarning, match="basis dimension reset"): + assert build_duchon_spline_setup(X, k=1, m=[1, 0.5]).bs_dim == 2 + with pytest.warns(UserWarning, match="more knots than data"): + ignored = build_duchon_spline_setup( + X, k=8, m=[1, 0.5], knots=[np.arange(31), np.arange(31)] + ) + assert not ignored.used_supplied_knots + + with pytest.raises(ValueError, match="same length"): + build_duchon_spline_setup( + X, + k=8, + m=[1, 0.5], + knots=[np.linspace(-1, 1, 10), np.linspace(-1, 1, 9)], + ) + with pytest.raises(ValueError, match="fewer unique covariate combinations"): + build_duchon_spline_setup(np.repeat(X[:4], 5, axis=0), k=8, m=[1, 0.5]) + with pytest.raises(ValueError, match="at least as many knot locations"): + build_duchon_spline_setup( + X, + k=8, + m=[1, 0.5], + knots=[np.arange(7), np.arange(7)], + ) + + +def test_ds_public_derivative_is_explicitly_unsupported(): + data = _ds_data(seed=279, n=80) + model = GAM( + formula='y ~ s(x0, bs="ds", k=10, sp=.7)', + optimize_smoothing=False, + ).fit(data=data) + with pytest.raises(NotImplementedError, match="derivative provider"): + model.derivative(data, smooth_number=1) diff --git a/tests/reference_fixtures/mgcv/109281e33f7454eb15886da8e37b1aa10ac1f3559e300a65be0970927e29c023.json.gz b/tests/reference_fixtures/mgcv/109281e33f7454eb15886da8e37b1aa10ac1f3559e300a65be0970927e29c023.json.gz new file mode 100644 index 00000000..eb433f95 Binary files /dev/null and b/tests/reference_fixtures/mgcv/109281e33f7454eb15886da8e37b1aa10ac1f3559e300a65be0970927e29c023.json.gz differ diff --git a/tests/reference_fixtures/mgcv/1fa370fcfa19a1de003ece7b6188e834e24c0906fed269de0ad2c2f1a2d84b74.json.gz b/tests/reference_fixtures/mgcv/1fa370fcfa19a1de003ece7b6188e834e24c0906fed269de0ad2c2f1a2d84b74.json.gz new file mode 100644 index 00000000..cbb7ad32 Binary files /dev/null and b/tests/reference_fixtures/mgcv/1fa370fcfa19a1de003ece7b6188e834e24c0906fed269de0ad2c2f1a2d84b74.json.gz differ diff --git a/tests/reference_fixtures/mgcv/2125a5695c4deb7922a1706fc3ed67ea1360efc14344263a0ef36f82f6a6f55b.json.gz b/tests/reference_fixtures/mgcv/2125a5695c4deb7922a1706fc3ed67ea1360efc14344263a0ef36f82f6a6f55b.json.gz new file mode 100644 index 00000000..02a40d0c Binary files /dev/null and b/tests/reference_fixtures/mgcv/2125a5695c4deb7922a1706fc3ed67ea1360efc14344263a0ef36f82f6a6f55b.json.gz differ diff --git a/tests/reference_fixtures/mgcv/2c27f773ad2eca61ed1842cec1f36c8a05e9bb1538be345dbdfe75cfd7587ac2.json.gz b/tests/reference_fixtures/mgcv/2c27f773ad2eca61ed1842cec1f36c8a05e9bb1538be345dbdfe75cfd7587ac2.json.gz new file mode 100644 index 00000000..853e98ab Binary files /dev/null and b/tests/reference_fixtures/mgcv/2c27f773ad2eca61ed1842cec1f36c8a05e9bb1538be345dbdfe75cfd7587ac2.json.gz differ diff --git a/tests/reference_fixtures/mgcv/2ce7b2f636eb1811403a236b997ecd2942f085f379f7ea4fa391794e6a5f3e59.json.gz b/tests/reference_fixtures/mgcv/2ce7b2f636eb1811403a236b997ecd2942f085f379f7ea4fa391794e6a5f3e59.json.gz new file mode 100644 index 00000000..8076ab6b Binary files /dev/null and b/tests/reference_fixtures/mgcv/2ce7b2f636eb1811403a236b997ecd2942f085f379f7ea4fa391794e6a5f3e59.json.gz differ diff --git a/tests/reference_fixtures/mgcv/2d2a3178a33e0dd92cbd8acf7374714fda9c31624f16aee2f47459dc4da757d6.json.gz b/tests/reference_fixtures/mgcv/2d2a3178a33e0dd92cbd8acf7374714fda9c31624f16aee2f47459dc4da757d6.json.gz new file mode 100644 index 00000000..943d4c4e Binary files /dev/null and b/tests/reference_fixtures/mgcv/2d2a3178a33e0dd92cbd8acf7374714fda9c31624f16aee2f47459dc4da757d6.json.gz differ diff --git a/tests/reference_fixtures/mgcv/2e370465f98fd5cbde02abc2da383bb001da3c5aecb21834c3c847f44ab19f09.json.gz b/tests/reference_fixtures/mgcv/2e370465f98fd5cbde02abc2da383bb001da3c5aecb21834c3c847f44ab19f09.json.gz new file mode 100644 index 00000000..eb57e055 Binary files /dev/null and b/tests/reference_fixtures/mgcv/2e370465f98fd5cbde02abc2da383bb001da3c5aecb21834c3c847f44ab19f09.json.gz differ diff --git a/tests/reference_fixtures/mgcv/32e463f3b3fea3bc98ff3b2d3ea254fb483a6256c1df2c697d4f9397ac49991e.json.gz b/tests/reference_fixtures/mgcv/32e463f3b3fea3bc98ff3b2d3ea254fb483a6256c1df2c697d4f9397ac49991e.json.gz new file mode 100644 index 00000000..edcefa24 Binary files /dev/null and b/tests/reference_fixtures/mgcv/32e463f3b3fea3bc98ff3b2d3ea254fb483a6256c1df2c697d4f9397ac49991e.json.gz differ diff --git a/tests/reference_fixtures/mgcv/32ecbb4b2249ee3552436710e71f9971ef6639f3d1faee978367c41fa974710e.json.gz b/tests/reference_fixtures/mgcv/32ecbb4b2249ee3552436710e71f9971ef6639f3d1faee978367c41fa974710e.json.gz new file mode 100644 index 00000000..254c88f9 Binary files /dev/null and b/tests/reference_fixtures/mgcv/32ecbb4b2249ee3552436710e71f9971ef6639f3d1faee978367c41fa974710e.json.gz differ diff --git a/tests/reference_fixtures/mgcv/3e98777f20fed2008c26ab798383ef914c611c8023013585c7467c4bbf5612f9.json.gz b/tests/reference_fixtures/mgcv/3e98777f20fed2008c26ab798383ef914c611c8023013585c7467c4bbf5612f9.json.gz new file mode 100644 index 00000000..4a986df1 Binary files /dev/null and b/tests/reference_fixtures/mgcv/3e98777f20fed2008c26ab798383ef914c611c8023013585c7467c4bbf5612f9.json.gz differ diff --git a/tests/reference_fixtures/mgcv/563fa6a632f670b03dba9092d2aca3677a9ef5c3b9f9b3ebebf76a210448b5dd.json.gz b/tests/reference_fixtures/mgcv/563fa6a632f670b03dba9092d2aca3677a9ef5c3b9f9b3ebebf76a210448b5dd.json.gz new file mode 100644 index 00000000..7fe85949 Binary files /dev/null and b/tests/reference_fixtures/mgcv/563fa6a632f670b03dba9092d2aca3677a9ef5c3b9f9b3ebebf76a210448b5dd.json.gz differ diff --git a/tests/reference_fixtures/mgcv/67433ec0570024ca6f6344c660655af8fec3dfdeccb40b0b85f6367c99a1098f.json.gz b/tests/reference_fixtures/mgcv/67433ec0570024ca6f6344c660655af8fec3dfdeccb40b0b85f6367c99a1098f.json.gz new file mode 100644 index 00000000..7786fa21 Binary files /dev/null and b/tests/reference_fixtures/mgcv/67433ec0570024ca6f6344c660655af8fec3dfdeccb40b0b85f6367c99a1098f.json.gz differ diff --git a/tests/reference_fixtures/mgcv/67fbdd9bd236fbd3e38f4767af5f4055a63e01ee89ca4ef2cc3c4d98c551593d.json.gz b/tests/reference_fixtures/mgcv/67fbdd9bd236fbd3e38f4767af5f4055a63e01ee89ca4ef2cc3c4d98c551593d.json.gz new file mode 100644 index 00000000..43561135 Binary files /dev/null and b/tests/reference_fixtures/mgcv/67fbdd9bd236fbd3e38f4767af5f4055a63e01ee89ca4ef2cc3c4d98c551593d.json.gz differ diff --git a/tests/reference_fixtures/mgcv/6e0df58b91a98e9148951fc5dff6771cbb33aaad32bd18e3f8bec81335ad442c.json.gz b/tests/reference_fixtures/mgcv/6e0df58b91a98e9148951fc5dff6771cbb33aaad32bd18e3f8bec81335ad442c.json.gz new file mode 100644 index 00000000..bf47d452 Binary files /dev/null and b/tests/reference_fixtures/mgcv/6e0df58b91a98e9148951fc5dff6771cbb33aaad32bd18e3f8bec81335ad442c.json.gz differ diff --git a/tests/reference_fixtures/mgcv/6e1aabe661b441300ec9faae1b17accfc4e710387a33de99795da5ef7c37d3fd.json.gz b/tests/reference_fixtures/mgcv/6e1aabe661b441300ec9faae1b17accfc4e710387a33de99795da5ef7c37d3fd.json.gz new file mode 100644 index 00000000..be19a5e9 Binary files /dev/null and b/tests/reference_fixtures/mgcv/6e1aabe661b441300ec9faae1b17accfc4e710387a33de99795da5ef7c37d3fd.json.gz differ diff --git a/tests/reference_fixtures/mgcv/7141d1b4bc886086f99d3860e765c92d7ff00dc3e579ee5edcdc1106114ff1c2.json.gz b/tests/reference_fixtures/mgcv/7141d1b4bc886086f99d3860e765c92d7ff00dc3e579ee5edcdc1106114ff1c2.json.gz new file mode 100644 index 00000000..1cc37d77 Binary files /dev/null and b/tests/reference_fixtures/mgcv/7141d1b4bc886086f99d3860e765c92d7ff00dc3e579ee5edcdc1106114ff1c2.json.gz differ diff --git a/tests/reference_fixtures/mgcv/72bb316923340a3177f52f73cffc6374fc513437f049c1766db08962de3675bb.json.gz b/tests/reference_fixtures/mgcv/72bb316923340a3177f52f73cffc6374fc513437f049c1766db08962de3675bb.json.gz new file mode 100644 index 00000000..f15336d9 Binary files /dev/null and b/tests/reference_fixtures/mgcv/72bb316923340a3177f52f73cffc6374fc513437f049c1766db08962de3675bb.json.gz differ diff --git a/tests/reference_fixtures/mgcv/7eb2f48bc03a50bc54acf889ccc1a20f3f1cb0d2cfba70cc52acab5ccf6c9ceb.json.gz b/tests/reference_fixtures/mgcv/7eb2f48bc03a50bc54acf889ccc1a20f3f1cb0d2cfba70cc52acab5ccf6c9ceb.json.gz new file mode 100644 index 00000000..7bfd75db Binary files /dev/null and b/tests/reference_fixtures/mgcv/7eb2f48bc03a50bc54acf889ccc1a20f3f1cb0d2cfba70cc52acab5ccf6c9ceb.json.gz differ diff --git a/tests/reference_fixtures/mgcv/7f015605a167c8e2aef6de40958c6fcbdd2966a5d1babc236e209bb069c6220f.json.gz b/tests/reference_fixtures/mgcv/7f015605a167c8e2aef6de40958c6fcbdd2966a5d1babc236e209bb069c6220f.json.gz new file mode 100644 index 00000000..7786fa21 Binary files /dev/null and b/tests/reference_fixtures/mgcv/7f015605a167c8e2aef6de40958c6fcbdd2966a5d1babc236e209bb069c6220f.json.gz differ diff --git a/tests/reference_fixtures/mgcv/93ce7f8090ed9f70905ce30c70bc91fe96e85cba30643a33a5b33c2f4304c2be.json.gz b/tests/reference_fixtures/mgcv/93ce7f8090ed9f70905ce30c70bc91fe96e85cba30643a33a5b33c2f4304c2be.json.gz new file mode 100644 index 00000000..a9018ee2 Binary files /dev/null and b/tests/reference_fixtures/mgcv/93ce7f8090ed9f70905ce30c70bc91fe96e85cba30643a33a5b33c2f4304c2be.json.gz differ diff --git a/tests/reference_fixtures/mgcv/987464d68c2bb0b0fef33a64cbcd511c50b30ab25e898619c0394f100794ca1c.json.gz b/tests/reference_fixtures/mgcv/987464d68c2bb0b0fef33a64cbcd511c50b30ab25e898619c0394f100794ca1c.json.gz new file mode 100644 index 00000000..c0c9a84c Binary files /dev/null and b/tests/reference_fixtures/mgcv/987464d68c2bb0b0fef33a64cbcd511c50b30ab25e898619c0394f100794ca1c.json.gz differ diff --git a/tests/reference_fixtures/mgcv/a1c30fbc57a05472809764c5a66aef8fa982daee3f214b165600dae39fd8e0e5.json.gz b/tests/reference_fixtures/mgcv/a1c30fbc57a05472809764c5a66aef8fa982daee3f214b165600dae39fd8e0e5.json.gz new file mode 100644 index 00000000..4c2bf269 Binary files /dev/null and b/tests/reference_fixtures/mgcv/a1c30fbc57a05472809764c5a66aef8fa982daee3f214b165600dae39fd8e0e5.json.gz differ diff --git a/tests/reference_fixtures/mgcv/ae134509bb4c2d6efa6d8207bc733ffdbef0b7e418cfeddd25b9aff47ac3e080.json.gz b/tests/reference_fixtures/mgcv/ae134509bb4c2d6efa6d8207bc733ffdbef0b7e418cfeddd25b9aff47ac3e080.json.gz new file mode 100644 index 00000000..0d8b3978 Binary files /dev/null and b/tests/reference_fixtures/mgcv/ae134509bb4c2d6efa6d8207bc733ffdbef0b7e418cfeddd25b9aff47ac3e080.json.gz differ diff --git a/tests/reference_fixtures/mgcv/b101d0d3317b2aee7ace5cbb5accf94e69397d42dab6e5cc40cedfbf56247faf.json.gz b/tests/reference_fixtures/mgcv/b101d0d3317b2aee7ace5cbb5accf94e69397d42dab6e5cc40cedfbf56247faf.json.gz new file mode 100644 index 00000000..e8887d05 Binary files /dev/null and b/tests/reference_fixtures/mgcv/b101d0d3317b2aee7ace5cbb5accf94e69397d42dab6e5cc40cedfbf56247faf.json.gz differ diff --git a/tests/reference_fixtures/mgcv/b150f200737f8d486eb88e7e3f9b6d97ca8b507cc1898a179f0615efa01d953b.json.gz b/tests/reference_fixtures/mgcv/b150f200737f8d486eb88e7e3f9b6d97ca8b507cc1898a179f0615efa01d953b.json.gz new file mode 100644 index 00000000..9ff9c41d Binary files /dev/null and b/tests/reference_fixtures/mgcv/b150f200737f8d486eb88e7e3f9b6d97ca8b507cc1898a179f0615efa01d953b.json.gz differ diff --git a/tests/reference_fixtures/mgcv/cc13e9a2d53fae51db18942f5077d1acc48dc2cf54774026b6f6ec2b51a8f224.json.gz b/tests/reference_fixtures/mgcv/cc13e9a2d53fae51db18942f5077d1acc48dc2cf54774026b6f6ec2b51a8f224.json.gz new file mode 100644 index 00000000..48650dce Binary files /dev/null and b/tests/reference_fixtures/mgcv/cc13e9a2d53fae51db18942f5077d1acc48dc2cf54774026b6f6ec2b51a8f224.json.gz differ diff --git a/tests/reference_fixtures/mgcv/d3445a2fb257b4478c80557265f3d50fbe867bbe9bed5de30ca05d5b676dc78a.json.gz b/tests/reference_fixtures/mgcv/d3445a2fb257b4478c80557265f3d50fbe867bbe9bed5de30ca05d5b676dc78a.json.gz new file mode 100644 index 00000000..c2e42bd7 Binary files /dev/null and b/tests/reference_fixtures/mgcv/d3445a2fb257b4478c80557265f3d50fbe867bbe9bed5de30ca05d5b676dc78a.json.gz differ diff --git a/tests/reference_fixtures/mgcv/da1c060ba90221639e8e2f2dfe8f9afae9cced4d048737f0760ea4c22496e112.json.gz b/tests/reference_fixtures/mgcv/da1c060ba90221639e8e2f2dfe8f9afae9cced4d048737f0760ea4c22496e112.json.gz new file mode 100644 index 00000000..a7f2af00 Binary files /dev/null and b/tests/reference_fixtures/mgcv/da1c060ba90221639e8e2f2dfe8f9afae9cced4d048737f0760ea4c22496e112.json.gz differ diff --git a/tests/reference_fixtures/mgcv/e860e4d5cdf5da5f097b1080a43d9019dfe3d3701d0f95696c528af4037706f2.json.gz b/tests/reference_fixtures/mgcv/e860e4d5cdf5da5f097b1080a43d9019dfe3d3701d0f95696c528af4037706f2.json.gz new file mode 100644 index 00000000..c7a3f7ef Binary files /dev/null and b/tests/reference_fixtures/mgcv/e860e4d5cdf5da5f097b1080a43d9019dfe3d3701d0f95696c528af4037706f2.json.gz differ diff --git a/tests/smooths/test_mgcv_raw_constructor_parity.py b/tests/smooths/test_mgcv_raw_constructor_parity.py index b496ba83..72b3af25 100644 --- a/tests/smooths/test_mgcv_raw_constructor_parity.py +++ b/tests/smooths/test_mgcv_raw_constructor_parity.py @@ -28,6 +28,7 @@ ) from nampy.gam.smooths.univariate.bs import DerivativeBSplineTerm1D from nampy.gam.smooths.univariate.cr import CubicSplineTerm +from nampy.gam.smooths.univariate.ds import DuchonSplineTerm from nampy.gam.smooths.univariate.ps import PSplineTerm1D from nampy.gam.smooths.univariate.tp import ThinPlateSplineTerm from nampy.gam.specs.build import build_formula_model @@ -207,6 +208,18 @@ def _build(data): return _build +def _observed_row_knots(columns, n_knots: int): + cols = [str(col) for col in columns] + + def _build(data): + return { + col: np.asarray(data[col], dtype=np.float64)[: int(n_knots)].copy() + for col in cols + } + + return _build + + def _merge_knots_factories(*builders): def _build(data): out = {} @@ -571,6 +584,61 @@ def _build_tprs_case_matrix(): return cases +def _build_duchon_case_matrix(): + return [ + _case( + "ds_1d_default_k", + _factory(_make_univariate_data, seed=150, n=90), + 'y ~ s(x, bs="ds")', + atol=2e-8, + ), + _case( + "ds_2d_default_k", + _factory(_make_gaussian_data, seed=151, n=100), + 'y ~ s(x0, x1, bs="ds")', + atol=2e-8, + ), + _case( + "ds_2d_custom_orders", + _factory(_make_gaussian_data, seed=152, n=90), + 'y ~ s(x0, x1, bs="ds", k=10, m=[1.4, .74])', + atol=2e-8, + ), + _case( + "ds_2d_supplied_truncated", + _factory(_make_gaussian_data, seed=153, n=90), + 'y ~ s(x0, x1, bs="ds", k=10, m=[1, .5])', + atol=2e-8, + knots_factory=_observed_row_knots(["x0", "x1"], 14), + ), + _case( + "ds_2d_supplied_pure_knot", + _factory(_make_gaussian_data, seed=154, n=90), + 'y ~ s(x0, x1, bs="ds", k=10, m=[1, .5])', + atol=2e-8, + knots_factory=_observed_row_knots(["x0", "x1"], 10), + ), + _case( + "ds_2d_max_knots_xt", + _factory(_make_gaussian_data, seed=155, n=70), + 'y ~ s(x0, x1, bs="ds", k=10, xt={"max.knots": 14, "seed": 7})', + atol=2e-8, + ), + _case( + "ds_3d_basic", + _factory(_make_gaussian_data_3col, seed=156, n=100), + 'y ~ s(x0, x1, x2, bs="ds", k=15)', + atol=5e-8, + ), + _case( + "ds_3d_default_k", + _factory(_make_gaussian_data_3col, seed=157, n=120), + 'y ~ s(x0, x1, x2, bs="ds")', + atol=2e-6, + ), + ] + + def _build_re_case_matrix(): penalty_multi = { "S": [ @@ -708,7 +776,7 @@ def _build_factor_smooth_case_matrix(): ) ) - for base_bs in ["tp", "ts"]: + for base_bs in ["tp", "ts", "ds"]: cases.append( _case( f"fs_2d_base_{base_bs}", @@ -726,6 +794,16 @@ def _build_factor_smooth_case_matrix(): ) ) + cases.append( + _case( + "fs_2d_base_ds_supplied_knots", + _factory(_make_factorized_gaussian_data, seed=171, n=96), + 'y ~ s(f, x0, x1, bs="fs", xt="ds", k=10, m=[1, .5])', + atol=2e-7, + knots_factory=_observed_row_knots(["x0", "x1"], 14), + ) + ) + return cases @@ -813,6 +891,18 @@ def _build_tensor_case_matrix(): 'y ~ te(x0, x1, x2, d=[2, 1], bs=["cr", "cr"], k=[10, 5])', atol=1e-7, ), + _case( + "te_d_duchon_margin", + _factory(_make_gaussian_data_3col, seed=822, n=110), + 'y ~ te(x0, x1, x2, d=[2, 1], bs=["ds", "cr"], k=[10, 5], m=[[1, .5], None])', + atol=2e-7, + ), + _case( + "ti_d_duchon_margin", + _factory(_make_gaussian_data_3col, seed=823, n=110), + 'y ~ ti(x0, x1, x2, d=[2, 1], bs=["ds", "cr"], k=[10, 5], m=[[1, .5], None])', + atol=2e-7, + ), # Supplied knots on ti() (only te had a knots case before). _case( "ti_knots_cr_cs", @@ -839,6 +929,7 @@ def _build_tensor_case_matrix(): *_build_cp_case_matrix(), *_build_bs_case_matrix(), *_build_tprs_case_matrix(), + *_build_duchon_case_matrix(), *_build_re_case_matrix(), *_build_factor_smooth_case_matrix(), *_build_tensor_case_matrix(), @@ -1035,6 +1126,27 @@ def _serialize_tprs_raw(term): ) +def _serialize_duchon_raw(term): + setup = term._setup + basis = np.asarray(setup.basis_train, dtype=np.float64) + return _common_raw_state( + "duchon.spline", + basis, + [np.asarray(setup.penalty, dtype=np.float64)], + rank=int(setup.rank), + null_space_dim=int(setup.null_space_dim), + extra={ + "knt": np.asarray(setup.knots, dtype=np.float64), + "UZ": np.asarray(setup.UZ, dtype=np.float64), + "shift": np.asarray(setup.shift, dtype=np.float64), + "p_order": [int(setup.penalty_order), float(setup.shift_order)], + "used_supplied_knots": bool(setup.used_supplied_knots), + "used_subsampling": bool(setup.used_subsampling), + "pure_knot": bool(setup.knots.shape[0] == setup.bs_dim), + }, + ) + + def _serialize_re_raw(term): B = np.asarray(term._basis_train, dtype=np.float64) q = int(B.shape[1]) @@ -1260,6 +1372,8 @@ def _serialize_term_raw(term, X): return _serialize_ps_raw(term) if isinstance(term, ThinPlateSplineTerm): return _serialize_tprs_raw(term) + if isinstance(term, DuchonSplineTerm): + return _serialize_duchon_raw(term) if isinstance(term, RandomEffectTerm): return _serialize_re_raw(term) if isinstance(term, FSmoothInteractionTerm): diff --git a/tests/smooths/test_mgcv_smoothcon_parity.py b/tests/smooths/test_mgcv_smoothcon_parity.py index 54b15e01..c34b3ef0 100644 --- a/tests/smooths/test_mgcv_smoothcon_parity.py +++ b/tests/smooths/test_mgcv_smoothcon_parity.py @@ -47,6 +47,12 @@ def _sym_rank(S: np.ndarray) -> int: return int(np.sum(ev > tol)) +def _first_penalty(penalties): + if isinstance(penalties, dict): + return next(iter(penalties.values())) + return penalties[0] + + def _assert_sz_penalty_invariants( actual_design: np.ndarray, expected_design: np.ndarray, @@ -998,6 +1004,74 @@ def test_bs_four_knot_prediction_interval_matches_mgcv(self): ) +class TestDuchonSplineSmooth: + """Duchon regression-spline smoothCon parity against mgcv.""" + + @staticmethod + def _make_data(seed=198, n=130): + rng = np.random.default_rng(seed) + x0 = rng.uniform(-2.0, 2.0, size=n) + x1 = rng.uniform(-1.5, 1.5, size=n) + y = np.sin(1.2 * x0) + 0.3 * x1**2 + rng.normal(scale=0.12, size=n) + return pd.DataFrame({"y": y, "x0": x0, "x1": x1}) + + def test_ds_1d_smoothcon_basis_and_penalty_match_mgcv(self): + data = self._make_data(seed=198) + formula = 'y ~ s(x0, bs="ds", k=11, sp=.7)' + expression = 's(x0, bs="ds", k=11, sp=.7)' + design = _compile_formula_design(data, formula) + expected_x = _run_mgcv_smoothcon_matrix(data, expression) + expected_s = _run_mgcv_smoothcon_penalties( + data, expression, absorb_cons=True, scale_penalty=True + ) + actual_x = np.asarray(design.design_matrix, dtype=np.float64) + target_x = np.asarray(expected_x["X"], dtype=np.float64) + actual_s = np.asarray(design.compiled_penalties[0].matrix, dtype=np.float64) + target_s = np.asarray(_first_penalty(expected_s["S"]), dtype=np.float64) + + _assert_allclose_up_to_column_sign(actual_x, target_x, atol=2e-8, rtol=2e-8) + np.testing.assert_allclose( + penalty_spectrum(actual_s), + penalty_spectrum(target_s), + atol=2e-8, + rtol=2e-8, + ) + np.testing.assert_allclose( + penalized_response_operator(actual_x, [actual_s]), + penalized_response_operator(target_x, [target_s]), + atol=2e-8, + rtol=2e-8, + ) + + def test_ds_2d_custom_order_smoothcon_basis_and_penalty_match_mgcv(self): + data = self._make_data(seed=199) + formula = 'y ~ s(x0, x1, bs="ds", k=10, m=[1,.5], sp=.7)' + expression = 's(x0, x1, bs="ds", k=10, m=c(1,.5), sp=.7)' + design = _compile_formula_design(data, formula) + expected_x = _run_mgcv_smoothcon_matrix(data, expression) + expected_s = _run_mgcv_smoothcon_penalties( + data, expression, absorb_cons=True, scale_penalty=True + ) + actual_x = np.asarray(design.design_matrix, dtype=np.float64) + target_x = np.asarray(expected_x["X"], dtype=np.float64) + actual_s = np.asarray(design.compiled_penalties[0].matrix, dtype=np.float64) + target_s = np.asarray(_first_penalty(expected_s["S"]), dtype=np.float64) + + _assert_allclose_up_to_column_sign(actual_x, target_x, atol=2e-8, rtol=2e-8) + np.testing.assert_allclose( + penalty_spectrum(actual_s), + penalty_spectrum(target_s), + atol=2e-8, + rtol=2e-8, + ) + np.testing.assert_allclose( + penalized_response_operator(actual_x, [actual_s]), + penalized_response_operator(target_x, [target_s]), + atol=2e-8, + rtol=2e-8, + ) + + class TestPSplineSmooth(_SharedTestPSplineSmooth): """P-spline (bs='ps') standalone parity against mgcv."""