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