From 903b755333f66103efa4337bbf17126f5d13cff0 Mon Sep 17 00:00:00 2001 From: ChrisW09 <50968720+ChrisW09@users.noreply.github.com> Date: Mon, 27 Jul 2026 22:00:19 +0200 Subject: [PATCH] fix(splines): close the final B-spline knot span The degree-0 base of the Cox-de Boor recursion used the half-open span ``[k_j, k_{j+1})``, so the largest value in the data belonged to no span and every basis function evaluated to zero there. The training argmax encoded as an all-zero row, and because the default minmax scaler pins the range, so did every value at or beyond the training maximum at inference time. Treat the last span of non-zero width as closed. Picking the last *positive width* span matters: a clamped knot vector repeats its boundary knot ``degree + 1`` times, so the trailing spans are degenerate and closing one of those is cancelled by the zero denominators further up the recursion. Also clip inputs into the fitted range in both ``transform`` methods, matching ``BaseSplineTransformer.transform``, so out-of-range values evaluate on the boundary instead of falling off the basis entirely. The two verbatim copies of ``bspline_basis`` (in ``pspline`` and ``tensor_product``) are consolidated into ``pretab.core.knots`` so the fix cannot diverge again; both modules re-export the name for compatibility. Partition of unity now holds across the full range for degrees 1-3 in both families, where previously the final row summed to 0. Closes #12 Co-Authored-By: Claude Opus 5 (1M context) --- pretab/core/knots.py | 60 +++++++++++++++++++ pretab/transformers/splines/pspline.py | 25 ++++---- pretab/transformers/splines/tensor_product.py | 23 ++++--- tests/test_pspline_transformer.py | 33 ++++++++++ tests/test_tensorproduct_transformer.py | 35 +++++++++++ 5 files changed, 148 insertions(+), 28 deletions(-) diff --git a/pretab/core/knots.py b/pretab/core/knots.py index 233f7bf..8c7bf25 100644 --- a/pretab/core/knots.py +++ b/pretab/core/knots.py @@ -18,6 +18,7 @@ __all__ = [ "basis_to_knots", + "bspline_basis", "generate_internal_knots", "quantile_knots", "select_knots", @@ -26,6 +27,65 @@ ] +def _last_positive_span(knots: np.ndarray) -> int: + """Index of the final knot span with non-zero width, or ``-1`` if none. + + A clamped knot vector repeats its boundary knots ``degree + 1`` times, so + the trailing spans are degenerate (``knots[j] == knots[j + 1]``). The span + that actually reaches the right boundary is the last one with positive + width. + """ + for j in range(len(knots) - 2, -1, -1): + if knots[j] < knots[j + 1]: + return j + return -1 + + +def bspline_basis(x: np.ndarray, knots: np.ndarray, degree: int, i: int) -> np.ndarray: + """Evaluate the ``i``-th B-spline basis function of ``degree`` over ``knots``. + + Standard Cox-de Boor recursion, with one deliberate departure: the final + knot span of non-zero width is treated as **closed** rather than half-open. + Under the usual ``[k_j, k_{j+1})`` convention the largest value in the data + belongs to no span, so every basis function evaluates to zero there and the + row loses all of its signal. Closing that span keeps the basis a partition + of unity across the whole fitted range, right endpoint included. + + Parameters + ---------- + x : ndarray + Points at which to evaluate the basis function. + knots : ndarray + Full (padded) knot vector. + degree : int + Degree of the basis function. + i : int + Index of the basis function. + """ + return _bspline_basis(x, knots, degree, i, _last_positive_span(knots)) + + +def _bspline_basis(x: np.ndarray, knots: np.ndarray, degree: int, i: int, last: int) -> np.ndarray: + """Cox-de Boor recursion with ``last`` naming the span to treat as closed.""" + if degree == 0: + if i == last: + return ((knots[i] <= x) & (x <= knots[i + 1])).astype(float) + return ((knots[i] <= x) & (x < knots[i + 1])).astype(float) + + denom1 = knots[i + degree] - knots[i] + denom2 = knots[i + degree + 1] - knots[i + 1] + # A zero denominator marks a degenerate span from a repeated boundary knot; + # its contribution is zero by convention. + zero = np.zeros_like(x, dtype=float) + term1 = zero if denom1 == 0 else (x - knots[i]) / denom1 * _bspline_basis(x, knots, degree - 1, i, last) + term2 = ( + zero + if denom2 == 0 + else (knots[i + degree + 1] - x) / denom2 * _bspline_basis(x, knots, degree - 1, i + 1, last) + ) + return term1 + term2 + + def basis_to_knots(n_basis: int, degree: int) -> int: """Number of internal knots implied by ``n_basis`` basis functions of ``degree``.""" return max(0, n_basis - degree - 1) diff --git a/pretab/transformers/splines/pspline.py b/pretab/transformers/splines/pspline.py index 2c6fbf4..60e030a 100644 --- a/pretab/transformers/splines/pspline.py +++ b/pretab/transformers/splines/pspline.py @@ -3,23 +3,14 @@ from sklearn.utils.validation import check_is_fitted from ...core.exceptions import InvalidParamError +from ...core.knots import bspline_basis from ...core.params import UNSET from .mixins import SplineBasisMixin - -def bspline_basis(x, knots, degree, i): - if degree == 0: - return np.where((x >= knots[i]) & (x < knots[i + 1]), 1.0, 0.0) - else: - denom1 = knots[i + degree] - knots[i] - denom2 = knots[i + degree + 1] - knots[i + 1] - - term1 = 0.0 if denom1 == 0 else (x - knots[i]) / denom1 * bspline_basis(x, knots, degree - 1, i) - term2 = ( - 0.0 if denom2 == 0 else (knots[i + degree + 1] - x) / denom2 * bspline_basis(x, knots, degree - 1, i + 1) - ) - - return term1 + term2 +# ``bspline_basis`` used to be defined here and in ``tensor_product``; the two +# copies are now a single implementation in ``pretab.core.knots``, re-exported +# under the original name so existing imports keep working. +__all__ = ["PSplineTransformer", "bspline_basis"] class PSplineTransformer(SplineBasisMixin, TransformerMixin, BaseEstimator): @@ -194,7 +185,11 @@ def transform(self, X): all_basis = [] for i in range(X.shape[1]): - x = X[:, i] + # Clip into the fitted range so values outside it evaluate on the + # boundary rather than falling off every knot span and producing an + # all-zero row. Matches BaseSplineTransformer.transform. + knots = self.knots_[i] + x = np.clip(X[:, i], knots[0], knots[-1]) nb = len(self.knots_[i]) - self.degree - 1 basis = np.zeros((len(x), nb)) for j in range(nb): diff --git a/pretab/transformers/splines/tensor_product.py b/pretab/transformers/splines/tensor_product.py index ec5dd42..5eaf6fd 100644 --- a/pretab/transformers/splines/tensor_product.py +++ b/pretab/transformers/splines/tensor_product.py @@ -3,21 +3,14 @@ from sklearn.utils.validation import check_is_fitted from ...core.exceptions import InvalidParamError +from ...core.knots import bspline_basis from ...core.params import UNSET from .mixins import SplineBasisMixin - -def bspline_basis(x, knots, degree, i): - if degree == 0: - return ((knots[i] <= x) & (x < knots[i + 1])).astype(float) - else: - denom1 = knots[i + degree] - knots[i] - denom2 = knots[i + degree + 1] - knots[i + 1] - term1 = 0.0 if denom1 == 0 else (x - knots[i]) / denom1 * bspline_basis(x, knots, degree - 1, i) - term2 = ( - 0.0 if denom2 == 0 else (knots[i + degree + 1] - x) / denom2 * bspline_basis(x, knots, degree - 1, i + 1) - ) - return term1 + term2 +# ``bspline_basis`` used to be defined here and in ``pspline``; the two copies +# are now a single implementation in ``pretab.core.knots``, re-exported under +# the original name so existing imports keep working. +__all__ = ["TensorProductSplineTransformer", "bspline_basis"] class TensorProductSplineTransformer(SplineBasisMixin, TransformerMixin, BaseEstimator): @@ -226,7 +219,11 @@ def transform(self, X): bases = [] for d in range(self.dim_): - basis = self._basis_matrix(X[:, d], self.knots_[d]) + # Clip into the fitted range so values outside it evaluate on the + # boundary rather than falling off every knot span and producing an + # all-zero marginal (which would zero the whole tensor product). + knots = self.knots_[d] + basis = self._basis_matrix(np.clip(X[:, d], knots[0], knots[-1]), knots) bases.append(basis) n_samples = X.shape[0] diff --git a/tests/test_pspline_transformer.py b/tests/test_pspline_transformer.py index 8036e75..0a5a54d 100644 --- a/tests/test_pspline_transformer.py +++ b/tests/test_pspline_transformer.py @@ -81,3 +81,36 @@ def test_pspline_transform_requires_fit(): transformer.transform(np.random.rand(5, 1)) with pytest.raises(NotFittedError): transformer.get_penalty_matrix() + + +# --------------------------------------------------------------------------- # +# The B-spline basis must cover the whole fitted range, right endpoint included. +# +# The degree-0 recursion base used a half-open span, so the largest observed +# value belonged to no span and its entire row evaluated to zero. +# --------------------------------------------------------------------------- # +@pytest.mark.parametrize("degree", [1, 2, 3]) +def test_pspline_is_partition_of_unity_including_max(degree): + X = np.linspace(0, 1, 30).reshape(-1, 1) + Xt = PSplineTransformer(output_dim=8, degree=degree).fit_transform(X) + + np.testing.assert_allclose(Xt.sum(axis=1), 1.0) + + +def test_pspline_max_row_is_not_all_zero(): + X = np.linspace(0, 1, 30).reshape(-1, 1) + Xt = PSplineTransformer(output_dim=8).fit_transform(X) + + assert np.abs(Xt[-1]).sum() > 0 + + +def test_pspline_clips_out_of_range_input(): + X = np.linspace(0, 1, 30).reshape(-1, 1) + transformer = PSplineTransformer(output_dim=8).fit(X) + + out = transformer.transform(np.array([[-0.5], [0.5], [1.5]])) + + # Out-of-range values evaluate on the boundary instead of vanishing. + np.testing.assert_allclose(out.sum(axis=1), 1.0) + np.testing.assert_allclose(out[0], transformer.transform(np.array([[0.0]]))[0]) + np.testing.assert_allclose(out[2], transformer.transform(np.array([[1.0]]))[0]) diff --git a/tests/test_tensorproduct_transformer.py b/tests/test_tensorproduct_transformer.py index 9a5e31c..90b5d64 100644 --- a/tests/test_tensorproduct_transformer.py +++ b/tests/test_tensorproduct_transformer.py @@ -71,3 +71,38 @@ def test_tensorproduct_transform_requires_fit(): transformer = TensorProductSplineTransformer() with pytest.raises(NotFittedError): transformer.transform(np.random.rand(5, 2)) + + +# --------------------------------------------------------------------------- # +# Same endpoint guarantee as the p-spline: both families share one basis +# implementation in ``pretab.core.knots``. +# --------------------------------------------------------------------------- # +@pytest.mark.parametrize("degree", [1, 2, 3]) +def test_tensorproduct_is_partition_of_unity_including_max(degree): + X = np.linspace(0, 1, 30).reshape(-1, 1) + Xt = TensorProductSplineTransformer(output_dim=8, degree=degree).fit_transform(X) + + np.testing.assert_allclose(Xt.sum(axis=1), 1.0) + + +def test_tensorproduct_max_row_is_not_all_zero(): + X = np.linspace(0, 1, 30).reshape(-1, 1) + Xt = TensorProductSplineTransformer(output_dim=8).fit_transform(X) + + assert np.abs(Xt[-1]).sum() > 0 + + +def test_tensorproduct_clips_out_of_range_input(): + X = np.linspace(0, 1, 30).reshape(-1, 1) + transformer = TensorProductSplineTransformer(output_dim=8).fit(X) + + out = transformer.transform(np.array([[-0.5], [0.5], [1.5]])) + + np.testing.assert_allclose(out.sum(axis=1), 1.0) + + +def test_tensorproduct_multivariate_still_products_the_marginals(): + rng = np.random.default_rng(0) + transformer = TensorProductSplineTransformer(output_dim=5).fit(rng.random((20, 2))) + + assert transformer.transform(rng.random((7, 2))).shape == (7, 25)