From 2405400d4c7a86722f66d2edc4d6fc19c720b4af Mon Sep 17 00:00:00 2001 From: ChrisW09 <50968720+ChrisW09@users.noreply.github.com> Date: Fri, 16 Oct 2026 14:13:38 +0200 Subject: [PATCH] fix(splines): evaluate the I-spline basis in closed form ISplineTransformer integrated each M-spline with a fixed 200-point trapezoid grid over the full knot range, so a basis function whose support is narrower than one grid step came out identically zero or ramped up over the wrong interval. This happens on skewed data under the default quantile placement and with a larger output_dim, and it also broke the ordering I_0 >= I_1 >= ... >= I_{K-1} of every row. Each I-spline is now the antiderivative of its B-spline (itself a B-spline of degree p + 1) divided by its analytic integral (t_{k+p+1} - t_k) / (p + 1), which is exact for any knot spacing. Values below / above the knot range stay 0 / 1. A basis function with zero-width support (a repeated knot) becomes the step function at that knot, its limit as the support shrinks. Fixes #53 --- pretab/expansion/spline/i_spline.py | 68 +++++++++---------- .../spline/test_spline_expansions.py | 55 +++++++++++++++ 2 files changed, 89 insertions(+), 34 deletions(-) diff --git a/pretab/expansion/spline/i_spline.py b/pretab/expansion/spline/i_spline.py index 85f374c..1d89f35 100644 --- a/pretab/expansion/spline/i_spline.py +++ b/pretab/expansion/spline/i_spline.py @@ -75,45 +75,45 @@ def __init__( def _feature_suffix(self) -> str: return "is" - def _ispline_basis(self, x: np.ndarray, knots: np.ndarray, basis_idx: int) -> np.ndarray: + def _ispline_columns(self, x: np.ndarray, knots: np.ndarray) -> np.ndarray: """ - Compute a single I-spline basis function. - - The M-spline is evaluated on a fine grid, integrated with the trapezoidal - rule to build a cumulative integral, then linearly interpolated at the - requested points and normalized by the full-range integral. + Evaluate every I-spline basis function at ``x`` in closed form. + + ``I_k(x)`` is the integral of the normalized M-spline ``M_k`` from the left + boundary, i.e. the antiderivative of the B-spline ``B_k`` divided by its + full-range integral ``(t_{k+p+1} - t_k) / (p + 1)``. The antiderivative of + a B-spline is itself a B-spline of degree ``p + 1``, so the basis is exact + for any knot spacing (a fixed quadrature grid cannot resolve basis + functions whose support is narrower than the grid step). Values below / + above the knot range are 0 / 1. A degenerate basis function with + zero-width support (only possible for a knot vector with a repeated + knot, e.g. an interior knot that coincides with a boundary knot on tied + data) integrates to the step function at that knot, its limit as the + support shrinks. """ - x_min_knot = knots[0] - x_max_knot = knots[-1] + n_coef = len(knots) - self.degree - 1 + lower, upper = knots[0], knots[-1] + antiderivative = BSpline(knots, np.eye(n_coef), self.degree).antiderivative() + total = (knots[self.degree + 1 : self.degree + 1 + n_coef] - knots[:n_coef]) / (self.degree + 1) + + x_clipped = np.clip(x, lower, upper) + values = antiderivative(x_clipped) - antiderivative(lower) + # Every basis function is fully integrated at the right boundary; pin it + # instead of evaluating there, where a repeated boundary knot would put + # the evaluation on a zero-width interval. + values = np.where(x_clipped[:, None] >= upper, total, values) + live = total > 0 + values[:, live] = values[:, live] / total[live] + values[:, ~live] = (x_clipped[:, None] >= knots[:n_coef][~live]).astype(float) + return np.clip(values, 0.0, 1.0) + + def _ispline_basis(self, x: np.ndarray, knots: np.ndarray, basis_idx: int) -> np.ndarray: + """Compute a single I-spline basis function (see :meth:`_ispline_columns`).""" n_coef = len(knots) - self.degree - 1 if basis_idx >= n_coef: return np.zeros(len(x)) - - grid = np.linspace(x_min_knot, x_max_knot, 200) - coef = np.zeros(n_coef) - coef[basis_idx] = 1.0 - spline = BSpline(knots, coef, self.degree, extrapolate=False) - m_values = np.nan_to_num(spline(grid), nan=0.0) - - knot_span = knots[basis_idx + self.degree + 1] - knots[basis_idx] - if knot_span > 1e-10: - m_values = m_values * (self.degree + 1) / knot_span - - cumulative = np.zeros(len(grid)) - for i in range(1, len(grid)): - dx = grid[i] - grid[i - 1] - cumulative[i] = cumulative[i - 1] + 0.5 * (m_values[i - 1] + m_values[i]) * dx - - ispline_values = np.interp(x, grid, cumulative, left=0.0, right=cumulative[-1]) - - max_integral = cumulative[-1] - if max_integral > 1e-10: - ispline_values = ispline_values / max_integral - return ispline_values + return self._ispline_columns(np.asarray(x, dtype=float), knots)[:, basis_idx] def _design_matrix(self, x: np.ndarray, knots: np.ndarray) -> np.ndarray: - n_coef = len(knots) - self.degree - 1 - design = np.zeros((len(x), n_coef)) - for i in range(n_coef): - design[:, i] = self._ispline_basis(x, knots, i) + design = self._ispline_columns(np.asarray(x, dtype=float), knots) return np.nan_to_num(design, nan=0.0) diff --git a/tests/expansion/spline/test_spline_expansions.py b/tests/expansion/spline/test_spline_expansions.py index a55d9ca..d7876f1 100644 --- a/tests/expansion/spline/test_spline_expansions.py +++ b/tests/expansion/spline/test_spline_expansions.py @@ -132,6 +132,61 @@ def test_ispline_bounded_unit_interval(): assert np.all(Xt <= 1.0 + 1e-9) +def _exact_ispline(x, knots, degree): + """Reference I-spline: the normalized antiderivative of each B-spline basis function.""" + from scipy.interpolate import BSpline + + n_coef = len(knots) - degree - 1 + columns = [] + for i in range(n_coef): + antiderivative = BSpline(knots, np.eye(n_coef)[i], degree).antiderivative() + lower, upper = antiderivative(knots[0]), antiderivative(knots[-1]) + columns.append((antiderivative(x) - lower) / (upper - lower)) + return np.column_stack(columns) + + +@pytest.mark.parametrize( + ("draw", "output_dim"), + [ + ("lognormal(0, 2)", 6), + ("lognormal(0, 2)", 10), + ("lognormal(0, 2)", 30), + ("normal", 20), + ("normal", 50), + ("lognormal(0, 1)", 20), + ("exponential", 20), + ], +) +def test_ispline_matches_exact_integral_on_tight_knots(draw, output_dim): + """Regression guard for issue #53: knot spans narrower than a fixed quadrature + grid step produced all-zero and misplaced columns.""" + rng = np.random.default_rng(0) + X = { + "lognormal(0, 2)": lambda: rng.lognormal(0, 2, (2000, 1)), + "lognormal(0, 1)": lambda: rng.lognormal(0, 1, (2000, 1)), + "normal": lambda: rng.normal(size=(2000, 1)), + "exponential": lambda: rng.exponential(size=(2000, 1)), + }[draw]() + transformer = ISplineTransformer(output_dim=output_dim).fit(X) + Xt = transformer.transform(X) + + exact = _exact_ispline(X[:, 0], transformer.knots_[0], transformer.degree) + np.testing.assert_allclose(Xt, exact, atol=1e-10) + assert (Xt.max(axis=0) > 0).all() + # Every I-spline row is ordered I_0 >= I_1 >= ... >= I_{K-1}. + assert (np.diff(Xt, axis=1) <= 1e-12).all() + + +def test_ispline_reaches_zero_and_one_at_the_range_boundaries(): + """Every column, including those with very narrow support on skewed data, rises + from exactly 0 at the left boundary to exactly 1 at the right boundary.""" + X = np.random.default_rng(1).lognormal(0, 2, (2000, 1)) + transformer = ISplineTransformer(output_dim=10).fit(X) + boundary = transformer.transform(np.array([[X.min()], [X.max()]])) + np.testing.assert_allclose(boundary[0], 0.0, atol=1e-12) + np.testing.assert_allclose(boundary[1], 1.0, atol=1e-12) + + def test_ispline_shape_multi_feature(): rng = np.random.RandomState(2) X = rng.uniform(0, 5, size=(100, 2))