fix(splines): close the final B-spline knot span

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)

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):

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]

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])

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)