From 303c7cdec4524ae7baa16f274c3dfd6a6cc3cb4b Mon Sep 17 00:00:00 2001 From: ChrisW09 <50968720+ChrisW09@users.noreply.github.com> Date: Fri, 9 Oct 2026 13:08:05 +0200 Subject: [PATCH 01/18] fix(splines): evaluate the I-spline basis in closed form ISplineTransformer integrated each M-spline with a fixed 200-point trapezoid grid, so basis functions narrower than one grid step came out identically zero or ramped over the wrong interval. Each I-spline is now the B-spline antiderivative normalized by its analytic integral, which is exact for any knot spacing. Fixes #53 --- pretab/expansion/spline/i_spline.py | 68 +++++++++---------- .../spline/test_spline_expansions.py | 36 ++++++++++ 2 files changed, 70 insertions(+), 34 deletions(-) diff --git a/pretab/expansion/spline/i_spline.py b/pretab/expansion/spline/i_spline.py index 85f374c..b95c609 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. one fitted before boundary knots were de-duplicated) + 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..f686460 100644 --- a/tests/expansion/spline/test_spline_expansions.py +++ b/tests/expansion/spline/test_spline_expansions.py @@ -132,6 +132,42 @@ 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("output_dim", [6, 10, 30]) +def test_ispline_matches_exact_integral_on_tight_knots(output_dim): + """Regression guard for issue #53: knot spans narrower than a fixed quadrature + grid step produced all-zero and misplaced columns.""" + X = np.random.default_rng(0).lognormal(0, 2, (2000, 1)) + 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(): + X = np.random.default_rng(1).exponential(size=(500, 1)) + transformer = ISplineTransformer(output_dim=8).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)) From 4b2becd77c578c99d20b2ad43c3b3338118aef82 Mon Sep 17 00:00:00 2001 From: ChrisW09 <50968720+ChrisW09@users.noreply.github.com> Date: Fri, 9 Oct 2026 13:08:49 +0200 Subject: [PATCH 02/18] fix(splines): drop absolute clipping thresholds from the M-spline basis MSplineTransformer zeroed spans below 1e-6 and clipped the design to +/-1e3. M-spline values scale like 1/span, so these absolute limits made the basis depend on the feature's units: small-scale features were flattened and tiny-range ones became all zero. The clips are removed; only a basis function with degenerate support relative to the knot range is zeroed. Fixes #63 --- pretab/expansion/spline/m_spline.py | 21 +++++++++++-------- .../spline/test_spline_expansions.py | 21 +++++++++++++++++++ 2 files changed, 33 insertions(+), 9 deletions(-) diff --git a/pretab/expansion/spline/m_spline.py b/pretab/expansion/spline/m_spline.py index f98682e..beac100 100644 --- a/pretab/expansion/spline/m_spline.py +++ b/pretab/expansion/spline/m_spline.py @@ -74,7 +74,14 @@ def _feature_suffix(self) -> str: return "ms" def _mspline_basis(self, x: np.ndarray, knots: np.ndarray, basis_idx: int) -> np.ndarray: - """Compute a single M-spline basis function.""" + """Compute a single M-spline basis function. + + ``M_k = (p + 1) / (t_{k+p+1} - t_k) * B_k``, which integrates to one. The + values scale like ``1 / span`` and are left unclipped, so the basis does + not depend on the units of the feature. Only a basis function whose + support is degenerate relative to the knot range (zero width, possible + only for a knot vector with a repeated knot) is returned as zeros. + """ n_coef = len(knots) - self.degree - 1 coef = np.zeros(n_coef) coef[basis_idx] = 1.0 @@ -82,17 +89,13 @@ def _mspline_basis(self, x: np.ndarray, knots: np.ndarray, basis_idx: int) -> np values = np.nan_to_num(spline(x), nan=0.0) knot_span = knots[basis_idx + self.degree + 1] - knots[basis_idx] - if knot_span > 1e-6: - values = values * (self.degree + 1) / knot_span - values = np.clip(values, 0.0, 1e6) - else: - values = np.zeros(len(x)) - return np.maximum(values, 0.0) + if knot_span <= np.finfo(float).eps * (knots[-1] - knots[0]): + return np.zeros(len(x)) + return np.maximum(values * (self.degree + 1) / knot_span, 0.0) 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._mspline_basis(x, knots, i) - design = np.nan_to_num(design, nan=0.0, posinf=0.0, neginf=0.0) - return np.clip(design, -1e3, 1e3) + 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 f686460..38ea023 100644 --- a/tests/expansion/spline/test_spline_expansions.py +++ b/tests/expansion/spline/test_spline_expansions.py @@ -116,6 +116,27 @@ def test_mspline_handles_nan(): assert np.isfinite(Xt).all() +@pytest.mark.parametrize("scale", [1.0, 1e-2, 1e-3, 1e-7]) +def test_mspline_basis_integrates_to_one_at_any_feature_scale(scale): + """Regression guard for issue #63: absolute clipping thresholds flattened the + basis on small-scale features and zeroed it on tiny-range ones.""" + from scipy.integrate import trapezoid + + X = np.random.default_rng(0).uniform(0, 1, (2000, 1)) * scale + transformer = MSplineTransformer(output_dim=8, placement_strategy="uniform").fit(X) + grid = np.linspace(X.min(), X.max(), 200_001) + integrals = trapezoid(transformer.transform(grid.reshape(-1, 1)), grid, axis=0) + np.testing.assert_allclose(integrals, 1.0, atol=1e-3) + + +def test_mspline_basis_is_scale_equivariant(): + X = np.random.default_rng(0).lognormal(0, 2, (2000, 1)) + unit = MSplineTransformer(output_dim=8).fit(X).transform(X) + scaled = MSplineTransformer(output_dim=8).fit(X * 1e-4).transform(X * 1e-4) + # M-spline values scale like 1 / span, so rescaling x by c rescales M by 1 / c. + np.testing.assert_allclose(scaled * 1e-4, unit, rtol=1e-8) + + def test_ispline_monotonic_increasing(): X = np.linspace(0, 10, 200).reshape(-1, 1) transformer = ISplineTransformer(output_dim=8, include_bias=False) From 9a9d173f5752a9f5665a80af412d5a461a1d1c9d Mon Sep 17 00:00:00 2001 From: ChrisW09 <50968720+ChrisW09@users.noreply.github.com> Date: Fri, 9 Oct 2026 13:12:11 +0200 Subject: [PATCH 03/18] fix(splines): keep knots unique and strictly interior on tied data Quantile knots on top-coded, zero-inflated or low-cardinality features repeat and land on the range boundary. The B/M/I base clipped knots into the closed range, and the natural-cubic, cubic-regression and tensor-product paths used them verbatim, so basis functions collapsed to a point: dead or duplicate columns, and all-zero B/M-spline rows at x_max. A shared supplement_interior_knots helper now drops boundary and repeated knots and tops the set up from strictly interior quantile and uniform candidates, keeping every existing knot. Both spline construction paths use it. Fixes #56 --- pretab/core/knots.py | 41 ++++++ pretab/expansion/spline/base.py | 36 ++---- pretab/expansion/spline/mixins.py | 17 ++- .../spline/test_tied_knot_placement.py | 118 ++++++++++++++++++ 4 files changed, 183 insertions(+), 29 deletions(-) create mode 100644 tests/expansion/spline/test_tied_knot_placement.py diff --git a/pretab/core/knots.py b/pretab/core/knots.py index 3a5b7b4..3bf81c2 100644 --- a/pretab/core/knots.py +++ b/pretab/core/knots.py @@ -23,6 +23,7 @@ "quantile_knots", "select_knots", "spanning_knots", + "supplement_interior_knots", "uniform_knots", ] @@ -142,3 +143,43 @@ def select_knots(knots: np.ndarray, count: int) -> np.ndarray: return knots idx = np.linspace(0, len(knots) - 1, count).round().astype(int) return knots[idx] + + +def supplement_interior_knots(x: np.ndarray, knots: np.ndarray, count: int) -> np.ndarray: + """Return unique knots strictly inside the range of ``x``, topped up to ``count``. + + Knots equal to (or outside) ``min(x)`` / ``max(x)`` and repeated knots give a + spline basis function with zero-width support -- a dead column -- so they are + dropped first. Every remaining knot is kept; when fewer than ``count`` remain, + the shortfall is filled with evenly spread quantile and uniform candidates + that are themselves strictly interior and not already present. On heavily tied + data the quantile candidates coincide with the tied values (often the range + boundary), so the uniform candidates guarantee the top-up for any feature with + a positive range. Never down-samples: callers trim an overfull set themselves. + + Parameters + ---------- + x : ndarray + Values of a single feature (finite). + knots : ndarray + Candidate interior knots, e.g. quantile knots or selector locations. + count : int + Minimum number of interior knots to return. + + Returns + ------- + ndarray + Sorted, unique knots strictly inside ``(min(x), max(x))``. + """ + x = np.asarray(x, dtype=float) + x_min, x_max = x.min(), x.max() + knots = np.asarray(knots, dtype=float) + knots = np.unique(knots[(knots > x_min) & (knots < x_max)]) + missing = count - len(knots) + if missing <= 0: + return knots + + candidates = np.unique(np.concatenate([quantile_knots(x, count), uniform_knots(x, count)])) + candidates = candidates[(candidates > x_min) & (candidates < x_max)] + candidates = np.setdiff1d(candidates, knots) + return np.unique(np.concatenate([knots, select_knots(candidates, missing)])) diff --git a/pretab/expansion/spline/base.py b/pretab/expansion/spline/base.py index ab365bf..394429f 100644 --- a/pretab/expansion/spline/base.py +++ b/pretab/expansion/spline/base.py @@ -23,9 +23,8 @@ from ...core.knots import ( basis_to_knots, generate_internal_knots, - quantile_knots, select_knots, - uniform_knots, + supplement_interior_knots, ) from ...core.parameters import UNSET, validate_placement from ...core.policy import RepresentationPolicy, resolve_out_of_range @@ -195,30 +194,22 @@ def _generate_knots(self, x: np.ndarray, n_knots: int, strategy: str) -> np.ndar def _adjust_internal_knots( self, x: np.ndarray, internal_knots: np.ndarray, min_knots: int, max_knots: int ) -> np.ndarray: - """Clip, deduplicate and rebalance internal knots to the allowed count.""" - internal_knots = np.clip(internal_knots, x.min(), x.max()) - internal_knots = np.unique(np.sort(internal_knots)) - if len(internal_knots) < min_knots: - internal_knots = self._supplement_knots(x, internal_knots, min_knots) + """Restrict, deduplicate and rebalance internal knots to the allowed count. + + Interior knots must lie strictly inside the data range and be unique: a + knot on the boundary or a repeated knot collapses a basis function's + support to a point, giving a dead column (and, at the boundary, all-zero + B/M-spline rows at ``x_max``). Such knots are dropped and the set is + topped up to ``min_knots`` before an overfull set is trimmed. + """ + internal_knots = self._supplement_knots(x, internal_knots, min_knots) if len(internal_knots) > max_knots: internal_knots = select_knots(internal_knots, max_knots) return internal_knots def _supplement_knots(self, x: np.ndarray, internal_knots: np.ndarray, target_count: int) -> np.ndarray: - """Add quantile and uniform candidates until ``target_count`` knots exist.""" - if target_count <= len(internal_knots): - return internal_knots - - candidates = [internal_knots] - if target_count > 0: - candidates.append(quantile_knots(x, target_count)) - candidates.append(uniform_knots(x, target_count)) - - combined = np.unique(np.concatenate(candidates)) - combined = np.sort(combined) - if len(combined) < target_count: - combined = uniform_knots(x, target_count) - return select_knots(np.asarray(combined), target_count) + """Keep the strictly interior, unique knots and top them up to ``target_count``.""" + return supplement_interior_knots(x, internal_knots, target_count) def _column_knots( self, @@ -255,9 +246,6 @@ def _column_knots( internal_knots = self._generate_knots(x_valid, n_internal, strategy) internal_knots = self._adjust_internal_knots(x_valid, internal_knots, min_knots, max_knots) - internal_knots = np.clip(internal_knots, x_min, x_max) - internal_knots = np.unique(np.sort(internal_knots)) - boundary_left = np.repeat(x_min, self.degree + 1) boundary_right = np.repeat(x_max, self.degree + 1) return np.concatenate([boundary_left, internal_knots, boundary_right]) diff --git a/pretab/expansion/spline/mixins.py b/pretab/expansion/spline/mixins.py index b325181..a73048c 100644 --- a/pretab/expansion/spline/mixins.py +++ b/pretab/expansion/spline/mixins.py @@ -13,7 +13,7 @@ import numpy as np from ...core.base import BasePreTabTransformer -from ...core.knots import generate_internal_knots, select_knots, spanning_knots +from ...core.knots import generate_internal_knots, select_knots, spanning_knots, supplement_interior_knots from ...exceptions import IncompatibleParamsError, PretabDataError @@ -84,13 +84,19 @@ def _place_spanning_knots(self, x, y, n_basis, strategy, selector, task, min_int window before bracketing. """ x, y = self._finite_column(x, y) + x_min, x_max = x.min(), x.max() if selector is not None: interior = self._place_interior_knots( x, y, n_basis - 2, strategy, selector, task, min_interior, max_interior ) - x_min, x_max = x.min(), x.max() return np.concatenate([[x_min], interior, [x_max]]) - return spanning_knots(x, n_basis, strategy) + knots = spanning_knots(x, n_basis, strategy) + if len(knots) <= 2: + return knots + # On tied data quantile knots repeat or land on the endpoints; keep the + # endpoints once and make the interior unique and strictly inside. + interior = supplement_interior_knots(x, knots[1:-1], n_basis - 2) + return np.concatenate([[x_min], interior, [x_max]]) def _place_interior_knots(self, x, y, n_interior, strategy, selector, task, min_interior=None, max_interior=None): """Return the interior knots (endpoints excluded) for one feature. @@ -106,7 +112,8 @@ def _place_interior_knots(self, x, y, n_interior, strategy, selector, task, min_ On the adaptive selector path ``min_interior`` / ``max_interior`` clamp the data-driven count into that window instead. Without a ``selector``, ``n_interior`` knots are placed with - :func:`pretab.core.knots.generate_internal_knots`. + :func:`pretab.core.knots.generate_internal_knots` and made unique and + strictly interior (see :func:`pretab.core.knots.supplement_interior_knots`). """ x, y = self._finite_column(x, y) if selector is not None: @@ -121,7 +128,7 @@ def _place_interior_knots(self, x, y, n_interior, strategy, selector, task, min_ min_interior = max_interior = n_interior selected = self._clamp_interior_knots(x, selected, min_interior, max_interior, strategy) return selected - return generate_internal_knots(x, n_interior, strategy) + return supplement_interior_knots(x, generate_internal_knots(x, n_interior, strategy), n_interior) def _clamp_interior_knots(self, x, knots, min_count, max_count, strategy): """Clamp a data-driven set of interior knots into ``[min_count, max_count]``. diff --git a/tests/expansion/spline/test_tied_knot_placement.py b/tests/expansion/spline/test_tied_knot_placement.py new file mode 100644 index 0000000..1345bc2 --- /dev/null +++ b/tests/expansion/spline/test_tied_knot_placement.py @@ -0,0 +1,118 @@ +"""Knot placement on tied data (top-coded, zero-inflated, low-cardinality features). + +Regression guards for issue #56: quantile knots on tied data coincide with each +other and with the range boundary, which collapsed basis functions to a point +(dead columns) and, for B/M-splines, left every row at ``x_max`` all zero. +""" + +import numpy as np +import pytest + +from pretab.core.knots import supplement_interior_knots +from pretab.transformers import ( + BSplineTransformer, + CubicRegressionSplineTransformer, + ISplineTransformer, + MSplineTransformer, + NaturalCubicSplineTransformer, + TensorProductSplineTransformer, +) + + +@pytest.fixture +def top_coded(): + # ~40% of the rows sit exactly on the cap of 100. + return np.minimum(np.random.default_rng(0).uniform(0, 160, (1000, 1)), 100.0) + + +@pytest.fixture +def zero_inflated(): + rng = np.random.default_rng(1) + return np.where(rng.random((1000, 1)) < 0.6, 0.0, rng.exponential(5.0, (1000, 1))) + + +def _interior(knots, degree): + return knots[degree + 1 : len(knots) - degree - 1] + + +def _dense_grid(X, n=2001): + return np.linspace(X.min(), X.max(), n).reshape(-1, 1) + + +@pytest.mark.parametrize("cls", [BSplineTransformer, MSplineTransformer, ISplineTransformer]) +@pytest.mark.parametrize("data", ["top_coded", "zero_inflated"]) +def test_bmi_interior_knots_are_unique_and_strictly_inside(cls, data, request): + X = request.getfixturevalue(data) + transformer = cls(output_dim=6).fit(X) + knots = transformer.knots_[0] + interior = _interior(knots, transformer.degree) + + assert len(interior) == 6 - transformer.degree - 1 + assert len(np.unique(interior)) == len(interior) + assert (interior > X.min()).all() and (interior < X.max()).all() + # Boundary multiplicity stays exactly degree + 1. + assert (knots == X.min()).sum() == transformer.degree + 1 + assert (knots == X.max()).sum() == transformer.degree + 1 + assert (np.abs(transformer.transform(_dense_grid(X))).max(axis=0) > 0).all() + + +def test_bspline_is_a_partition_of_unity_at_the_cap(top_coded): + transformer = BSplineTransformer(output_dim=6).fit(top_coded) + basis = transformer.transform(top_coded) + np.testing.assert_allclose(basis.sum(axis=1), 1.0, atol=1e-9) + # Values on and beyond the cap (clipped onto x_max) keep a full basis row. + np.testing.assert_allclose(transformer.transform(np.array([[100.0], [150.0]])).sum(axis=1), 1.0, atol=1e-9) + + +def test_ispline_equals_one_at_the_cap(top_coded): + transformer = ISplineTransformer(output_dim=6).fit(top_coded) + np.testing.assert_allclose(transformer.transform(np.array([[100.0]])), 1.0, atol=1e-12) + + +@pytest.mark.parametrize("cls", [NaturalCubicSplineTransformer, CubicRegressionSplineTransformer]) +def test_cubic_families_place_unique_knots_on_zero_inflated_data(cls, zero_inflated): + transformer = cls(output_dim=6, placement_strategy="quantile").fit(zero_inflated) + knots = transformer.knots_[0] + assert len(np.unique(knots)) == len(knots) + basis = transformer.transform(_dense_grid(zero_inflated)) + assert np.linalg.matrix_rank(basis) == basis.shape[1] + assert (np.abs(basis).max(axis=0) > 0).all() + + +def test_natural_spline_keeps_both_endpoints_once(zero_inflated): + knots = NaturalCubicSplineTransformer(output_dim=6, placement_strategy="quantile").fit(zero_inflated).knots_[0] + assert len(knots) == 7 + assert knots[0] == zero_inflated.min() and knots[-1] == zero_inflated.max() + assert (knots[1:-1] > knots[0]).all() and (knots[1:-1] < knots[-1]).all() + + +def test_tensor_product_has_no_degenerate_columns_on_zero_inflated_marginal(zero_inflated): + X = np.hstack([zero_inflated, np.random.default_rng(2).normal(size=(len(zero_inflated), 1))]) + transformer = TensorProductSplineTransformer(output_dim=8, placement_strategy="quantile").fit(X) + g0 = np.linspace(X[:, 0].min(), X[:, 0].max(), 200) + g1 = np.linspace(X[:, 1].min(), X[:, 1].max(), 200) + grid = np.array(np.meshgrid(g0, g1)).reshape(2, -1).T + assert (np.abs(transformer.transform(grid)).max(axis=0) > 0).all() + + +@pytest.mark.parametrize("levels", [[0.0, 1.0], [0.0, 1.0, 2.0], [1.0, 2.0, 3.0, 4.0]]) +def test_low_cardinality_feature_keeps_a_full_rank_bspline_basis(levels): + X = np.random.default_rng(3).choice(levels, size=(300, 1)) + transformer = BSplineTransformer(output_dim=6).fit(X) + basis = transformer.transform(_dense_grid(X)) + assert np.linalg.matrix_rank(basis) == 6 + + +def test_supplement_interior_knots_drops_boundary_and_duplicate_knots(): + x = np.array([0.0, 0.0, 0.0, 1.0, 2.0, 10.0, 10.0]) + knots = supplement_interior_knots(x, np.array([0.0, 2.0, 2.0, 10.0, 12.0]), 3) + assert len(knots) == 3 + assert 2.0 in knots # an existing interior knot is always kept + assert (knots > 0.0).all() and (knots < 10.0).all() + assert len(np.unique(knots)) == 3 + + +def test_supplement_interior_knots_never_down_samples(): + x = np.linspace(0, 1, 50) + knots = np.array([0.1, 0.2, 0.3, 0.4]) + np.testing.assert_array_equal(supplement_interior_knots(x, knots, 2), knots) From 1fd2a367abc1199dabc4bf38a70f6403a6b8e77c Mon Sep 17 00:00:00 2001 From: ChrisW09 <50968720+ChrisW09@users.noreply.github.com> Date: Fri, 9 Oct 2026 13:13:35 +0200 Subject: [PATCH 04/18] fix(placement): rank CART split candidates by impurity decrease before spacing CARTLocationSelector returned its split candidates in ascending location order and the minimum-spacing filter keeps the first of two close candidates, so a weak split just below the dominant one evicted it. Candidates are now ordered by weighted impurity decrease (as the LightGBM selector already does with gains), and the same ranking drives the over-max trim. Fixes #70 --- pretab/core/selectors.py | 87 ++++++++++++--------------- tests/core/test_location_selectors.py | 37 ++++++++++++ 2 files changed, 76 insertions(+), 48 deletions(-) diff --git a/pretab/core/selectors.py b/pretab/core/selectors.py index f4cef30..1cff033 100644 --- a/pretab/core/selectors.py +++ b/pretab/core/selectors.py @@ -20,7 +20,7 @@ """ from abc import ABC, abstractmethod -from typing import Literal +from typing import Literal, cast import numpy as np from sklearn.tree import DecisionTreeClassifier, DecisionTreeRegressor @@ -116,10 +116,11 @@ def select( def _ordered_candidates(self, x_valid: np.ndarray, y_valid: np.ndarray, task: Task) -> tuple[list[float], object]: """Fit a model and return candidate locations plus trimming context. - The candidates must be returned in the selector's preferred order (the - order :meth:`_enforce_spacing` should honour): location order for a - single tree, gain-descending order for a boosted ensemble. ``context`` is - an opaque object passed straight through to :meth:`_trim_over_max`. + The candidates must be returned in the selector's preferred order -- the + order :meth:`_enforce_spacing` honours, which keeps the earlier of two + candidates that are too close: impurity-decrease order for a single tree, + gain-descending order for a boosted ensemble. ``context`` is an opaque + object passed straight through to :meth:`_trim_over_max`. """ raise NotImplementedError @@ -131,9 +132,10 @@ def _trim_over_max(self, points: list[float], context: object, max_count: int) - def _enforce_spacing(self, split_points: list[float], x: np.ndarray) -> list[float]: """Drop locations closer than ``min_location_spacing`` of the range. - Compares each candidate against every already-kept location so the filter - is order-independent and honours both ascending (CART) and gain-descending - (LightGBM) ordering. + Compares each candidate against every already-kept location, so of two + candidates that are too close the one earlier in ``split_points`` -- the + more important one, given the importance-descending order both selectors + produce -- is kept. """ if len(split_points) <= 1: return split_points @@ -172,8 +174,9 @@ class CARTLocationSelector(BaseLocationSelector): A ``DecisionTreeRegressor`` or ``DecisionTreeClassifier`` is fitted to the feature against the target, and its split thresholds become the candidate - locations. Candidates are spaced out, and if there are too many they are - ranked by weighted impurity decrease so the most informative splits are kept. + locations. Candidates are ranked by weighted impurity decrease, so the most + informative splits are kept both when two candidates are too close to each + other and when there are too many. Parameters ---------- @@ -222,62 +225,50 @@ def _ordered_candidates(self, x_valid: np.ndarray, y_valid: np.ndarray, task: Ta ) tree.fit(x_valid, y_valid) - split_points = self._extract_split_points(tree, x_valid) - return split_points, tree + importance = self._split_importance(tree, x_valid) + # Most informative split first, so that when two candidates are closer + # than ``min_location_spacing`` the spacing filter keeps the stronger one. + return self._rank(importance, importance), importance def _trim_over_max(self, points: list[float], context: object, max_count: int) -> list[float]: - return self._select_top_locations(points, context, max_count) + importance = cast(dict[float, float], context) + return sorted(self._rank(points, importance)[:max_count]) - def _extract_split_points(self, tree, x: np.ndarray) -> list[float]: - """Collect in-range split thresholds from a fitted decision tree.""" - tree_structure = tree.tree_ - split_points = [] - - x_min, x_max = float(x.min()), float(x.max()) - - for node_id in range(tree_structure.node_count): - is_split = tree_structure.children_left[node_id] != tree_structure.children_right[node_id] - if not is_split: - continue - if tree_structure.feature[node_id] != 0: - continue - threshold = tree_structure.threshold[node_id] - if x_min < threshold < x_max: - split_points.append(threshold) - - return sorted(set(split_points)) + @staticmethod + def _rank(points, importance: dict[float, float]) -> list[float]: + """Order ``points`` by decreasing impurity decrease (ties by location).""" + return sorted(points, key=lambda point: (-importance[point], point)) - def _select_top_locations(self, candidates: list[float], tree, max_count: int) -> list[float]: - """Keep the locations whose splits reduce impurity the most.""" - if len(candidates) <= max_count: - return candidates + def _split_importance(self, tree, x: np.ndarray) -> dict[float, float]: + """Map each in-range split threshold to its weighted impurity decrease. + In a single-feature tree every threshold occurs at most once: after a split + at ``t`` no descendant holds samples on both sides of ``t``. + """ tree_structure = tree.tree_ - split_importance = {} + x_min, x_max = float(x.min()), float(x.max()) + split_importance: dict[float, float] = {} for node_id in range(tree_structure.node_count): - is_split = tree_structure.children_left[node_id] != tree_structure.children_right[node_id] - if not is_split: + left_child = tree_structure.children_left[node_id] + right_child = tree_structure.children_right[node_id] + if left_child == right_child or tree_structure.feature[node_id] != 0: continue - - threshold = tree_structure.threshold[node_id] - if threshold not in candidates: + threshold = float(tree_structure.threshold[node_id]) + if not x_min < threshold < x_max: continue n_samples = tree_structure.n_node_samples[node_id] impurity = tree_structure.impurity[node_id] - - left_child = tree_structure.children_left[node_id] - right_child = tree_structure.children_right[node_id] n_left = tree_structure.n_node_samples[left_child] n_right = tree_structure.n_node_samples[right_child] impurity_left = tree_structure.impurity[left_child] impurity_right = tree_structure.impurity[right_child] + split_importance[threshold] = float( + n_samples * impurity - (n_left * impurity_left + n_right * impurity_right) + ) - split_importance[threshold] = n_samples * impurity - (n_left * impurity_left + n_right * impurity_right) - - top = sorted(split_importance, key=lambda k: split_importance[k], reverse=True)[:max_count] - return sorted(top) + return split_importance class LightGBMLocationSelector(BaseLocationSelector): diff --git a/tests/core/test_location_selectors.py b/tests/core/test_location_selectors.py index e1e6e7f..ad00909 100644 --- a/tests/core/test_location_selectors.py +++ b/tests/core/test_location_selectors.py @@ -160,3 +160,40 @@ def test_supplement_preserves_high_end_split_end_to_end(): locations = CARTLocationSelector().select(xs.reshape(-1, 1), ys, task="regression", min_count=6, max_count=6) assert locations.max() > 5.0, f"expected a high-end location to survive supplementing, got {locations}" + + +def _single_step(seed): + rng = np.random.default_rng(seed) + x = rng.uniform(0, 10, 1000) + step = rng.uniform(2, 8) + y = np.where(x > step, 1.0, 0.0) + 0.5 * rng.normal(size=1000) + return x.reshape(-1, 1), y + + +def test_cart_spacing_keeps_the_dominant_split_over_a_weak_neighbour(): + """Regression guard for issue #70: candidates were spaced in ascending location + order, so a weak split just below the root split evicted it.""" + from sklearn.tree import DecisionTreeRegressor + + X, y = _single_step(34) + root = DecisionTreeRegressor(max_depth=1).fit(X, y).tree_.threshold[0] + locations = CARTLocationSelector().select(X, y, task="regression", min_count=6, max_count=6) + assert np.isclose(locations, root).any(), f"root split {root:.3f} missing from {locations}" + + +@pytest.mark.parametrize("seed", range(10)) +def test_cart_keeps_the_root_split_of_a_single_step_target(seed): + from sklearn.tree import DecisionTreeRegressor + + X, y = _single_step(seed) + root = DecisionTreeRegressor(max_depth=1).fit(X, y).tree_.threshold[0] + locations = CARTLocationSelector().select(X, y, task="regression", min_count=6, max_count=6) + assert np.isclose(locations, root).any() + + +def test_cart_candidates_are_ordered_by_impurity_decrease(data): + X, y = data + selector = CARTLocationSelector() + candidates, importance = selector._ordered_candidates(X, y, "regression") + gains = [importance[c] for c in candidates] + assert gains == sorted(gains, reverse=True) From 81a024be54688756f3a0da1552feb74e68b42b6f Mon Sep 17 00:00:00 2001 From: ChrisW09 <50968720+ChrisW09@users.noreply.github.com> Date: Fri, 9 Oct 2026 13:18:04 +0200 Subject: [PATCH 05/18] fix(placement): top up target-aware locations with distinct interior points on tied data On discrete or heavily tied features the selector's quantile top-up and its quantile fallbacks collapsed onto the tied values and the range boundary: too few locations (narrower than output_dim, breaking CrossFittedTransformer), duplicates, and locations on x_min/x_max that gave PLE a constant column. The top-up and both fallbacks now keep every selector-found location and fill the shortfall with strictly interior, spaced quantile candidates first, then uniform ones. LightGBM's +/-1e-35 split-at-zero sentinel is mapped to the midpoint of the values it separates, and PLE treats any positive bin width as non-degenerate instead of using an absolute 1e-10 tolerance. Fixes #57 --- pretab/core/selectors.py | 85 ++++++++++++++++++++----- pretab/encoding/numerical/ple.py | 4 +- tests/core/test_cross_fitted.py | 14 ++++ tests/core/test_feature_map_selector.py | 14 ++++ tests/core/test_location_selectors.py | 67 +++++++++++++++++++ tests/core/test_ple_selector.py | 12 ++++ 6 files changed, 177 insertions(+), 19 deletions(-) diff --git a/pretab/core/selectors.py b/pretab/core/selectors.py index 1cff033..dba14d9 100644 --- a/pretab/core/selectors.py +++ b/pretab/core/selectors.py @@ -26,10 +26,14 @@ from sklearn.tree import DecisionTreeClassifier, DecisionTreeRegressor from ..exceptions import IncompatibleParamsError, OptionalDependencyError -from .knots import quantile_knots, select_knots +from .knots import quantile_knots, select_knots, uniform_knots Task = Literal["regression", "classification"] +# Magnitude bound of LightGBM's "split at zero" sentinel threshold (its +# ``kZeroThreshold`` is 1e-35, reported as the float32 value ~1.0000000180e-35). +_LIGHTGBM_ZERO_THRESHOLD = 1e-34 + class BaseLocationSelector(ABC): """Abstract base class for count-based, target-aware location selectors. @@ -97,11 +101,11 @@ def select( y_valid = y[valid_mask] if len(x_valid) < self.min_samples_floor: - return quantile_knots(x_valid, min_count) + return np.array(self._supplement([], x_valid, min_count)) points, context = self._ordered_candidates(x_valid, y_valid, task) if len(points) == 0: - return quantile_knots(x_valid, min_count) + return np.array(self._supplement([], x_valid, min_count)) locations = self._enforce_spacing(points, x_valid) @@ -151,22 +155,44 @@ def _enforce_spacing(self, split_points: list[float], x: np.ndarray) -> list[flo return spaced def _supplement(self, existing: list[float], x: np.ndarray, target_count: int) -> list[float]: - """Top up an under-filled location set with quantile locations. + """Top up an under-filled location set to ``target_count`` distinct locations. Keeps every existing (selector-found) location and fills only the - shortfall with quantile candidates, rather than truncating the union - (which would preferentially drop the largest existing values). + shortfall, rather than truncating the union (which would preferentially + drop the largest existing values). Candidates lie strictly inside the + range of ``x`` and keep ``min_location_spacing`` from every kept location: + evenly spread quantile locations first, then uniform ones. On a tied or + discrete feature the quantiles collapse onto a few values (often the range + boundary), and the uniform candidates still provide distinct interior + locations for any feature with a positive range. Only when the requested + count is too dense for the spacing are the remaining uniform locations + added without it. """ missing = target_count - len(existing) if missing <= 0: return existing - existing_set = set(existing) - candidates = [c for c in quantile_knots(x, target_count).tolist() if c not in existing_set] - combined = sorted(existing_set | set(candidates[:missing])) - if len(combined) > target_count: - combined = select_knots(np.array(combined), target_count).tolist() - return combined + x = np.asarray(x, dtype=float).ravel() + x_min, x_max = float(x.min()), float(x.max()) + if x_max <= x_min: + # A zero-range feature has no interior; keep the historical repeated + # location so the requested count (and output width) still holds. + return sorted([*existing, *quantile_knots(x, missing).tolist()]) + + min_distance = self.min_location_spacing * (x_max - x_min) + kept = sorted(float(location) for location in existing) + for candidates in (quantile_knots(x, target_count), uniform_knots(x, target_count)): + eligible = [] + for candidate in np.unique(candidates[(candidates > x_min) & (candidates < x_max)]): + if all(abs(candidate - other) >= min_distance for other in (*kept, *eligible)): + eligible.append(float(candidate)) + kept = sorted(kept + select_knots(np.array(eligible), missing).tolist()) + missing = target_count - len(kept) + if missing <= 0: + return kept + + leftovers = np.setdiff1d(uniform_knots(x, target_count), kept) + return sorted(kept + select_knots(leftovers, missing).tolist()) class CARTLocationSelector(BaseLocationSelector): @@ -361,16 +387,39 @@ def _trim_over_max(self, points: list[float], context: object, max_count: int) - def _extract_split_points_with_gains(self, model, x: np.ndarray) -> dict[float, float]: """Collect split thresholds and their cumulative gains from a model.""" + x = np.asarray(x, dtype=float).ravel() x_min, x_max = float(x.min()), float(x.max()) split_importance: dict[float, float] = {} + # LightGBM reports a split at zero as the sentinel threshold +/-1e-35 + # rather than as a bin midpoint. Map it to the midpoint between the + # values the split separates, the location CART would report. + zero_splits = { + 1.0: self._midpoint(x[x <= 0], x[x > 0]), + -1.0: self._midpoint(x[x < 0], x[x >= 0]), + } + model_dict = model.dump_model() for tree_info in model_dict["tree_info"]: - self._traverse_tree(tree_info["tree_structure"], split_importance, x_min, x_max) + self._traverse_tree(tree_info["tree_structure"], split_importance, x_min, x_max, zero_splits) return split_importance - def _traverse_tree(self, node: dict, split_importance: dict, x_min: float, x_max: float): + @staticmethod + def _midpoint(left: np.ndarray, right: np.ndarray) -> float | None: + """Midpoint between the largest ``left`` and the smallest ``right`` value.""" + if left.size == 0 or right.size == 0: + return None + return float((left.max() + right.min()) / 2.0) + + def _traverse_tree( + self, + node: dict, + split_importance: dict, + x_min: float, + x_max: float, + zero_splits: dict[float, float | None] | None = None, + ): """Recursively accumulate split gains for the single feature.""" if "split_feature" not in node: return @@ -378,10 +427,12 @@ def _traverse_tree(self, node: dict, split_importance: dict, x_min: float, x_max if node["split_feature"] == 0: threshold = node["threshold"] gain = node.get("split_gain", 0.0) - if x_min < threshold < x_max: + if zero_splits is not None and abs(threshold) <= _LIGHTGBM_ZERO_THRESHOLD: + threshold = zero_splits[1.0 if threshold >= 0 else -1.0] + if threshold is not None and x_min < threshold < x_max: split_importance[threshold] = split_importance.get(threshold, 0.0) + gain if "left_child" in node: - self._traverse_tree(node["left_child"], split_importance, x_min, x_max) + self._traverse_tree(node["left_child"], split_importance, x_min, x_max, zero_splits) if "right_child" in node: - self._traverse_tree(node["right_child"], split_importance, x_min, x_max) + self._traverse_tree(node["right_child"], split_importance, x_min, x_max, zero_splits) diff --git a/pretab/encoding/numerical/ple.py b/pretab/encoding/numerical/ple.py index 81a639c..cca1e86 100644 --- a/pretab/encoding/numerical/ple.py +++ b/pretab/encoding/numerical/ple.py @@ -276,7 +276,7 @@ def _apply_piecewise_linear_vectorized( if len(thresholds) == 0: lower, upper = edges[0], edges[-1] width = upper - lower - if width > 1e-10: + if width > 0: values = np.clip((feature - lower) / width, 0.0, 1.0) else: values = np.full(n_samples, 0.5) @@ -297,7 +297,7 @@ def _apply_piecewise_linear_vectorized( upper_edge = edges[bin_idx + 1] bin_width = upper_edge - lower_edge - if bin_width > 1e-10: + if bin_width > 0: ple_encoded[mask, bin_idx] = np.clip((values - lower_edge) / bin_width, 0.0, 1.0) else: ple_encoded[mask, bin_idx] = 0.5 diff --git a/tests/core/test_cross_fitted.py b/tests/core/test_cross_fitted.py index 723d971..21d9480 100644 --- a/tests/core/test_cross_fitted.py +++ b/tests/core/test_cross_fitted.py @@ -104,3 +104,17 @@ def test_invalid_task_rejected(data): X, y = data with pytest.raises(InvalidParamError, match="task"): CrossFittedTransformer(PLETransformer(), task="classificaton").fit_transform(X, y) + + +@pytest.mark.parametrize("seed", range(5)) +def test_cross_fitting_a_discrete_feature_keeps_a_fixed_width(seed): + """Regression guard for issue #57: tied features gave fold-dependent widths.""" + from pretab.transformers import RBFExpansionTransformer + + rng = np.random.default_rng(seed) + x = rng.integers(0, 6, size=200).astype(float).reshape(-1, 1) + y = np.sin(x[:, 0]) + rng.normal(0, 0.5, size=200) + out = CrossFittedTransformer( + RBFExpansionTransformer(output_dim=10, target_aware=True), n_folds=5, random_state=0 + ).fit_transform(x, y) + assert out.shape == (200, 10) diff --git a/tests/core/test_feature_map_selector.py b/tests/core/test_feature_map_selector.py index 2b886c6..20d4ac5 100644 --- a/tests/core/test_feature_map_selector.py +++ b/tests/core/test_feature_map_selector.py @@ -100,3 +100,17 @@ def test_lightgbm_selector_places_centers(Cls, data): t = Cls(output_dim=5, target_aware=True, task="regression", placement_strategy="lightgbm").fit(X, y) assert all(len(c) == 5 for c in t.centers_) assert all(np.all(np.diff(c) > 0) for c in t.centers_) + + +@pytest.mark.parametrize("strategy", ["cart", "lightgbm"]) +def test_target_aware_centers_keep_output_dim_on_a_discrete_feature(strategy): + """Regression guard for issue #57: tied features returned fewer, duplicate centers.""" + if strategy == "lightgbm": + pytest.importorskip("lightgbm") + rng = np.random.default_rng(0) + x = rng.integers(0, 5, size=300).astype(float).reshape(-1, 1) + y = (x[:, 0] >= 1) + rng.normal(0, 0.1, size=300) + transformer = RBFExpansionTransformer(output_dim=10, target_aware=True, placement_strategy=strategy).fit(x, y) + centers = transformer.centers_[0] + assert transformer.transform(x).shape[1] == 10 + assert len(np.unique(centers)) == 10 diff --git a/tests/core/test_location_selectors.py b/tests/core/test_location_selectors.py index ad00909..49dcf01 100644 --- a/tests/core/test_location_selectors.py +++ b/tests/core/test_location_selectors.py @@ -197,3 +197,70 @@ def test_cart_candidates_are_ordered_by_impurity_decrease(data): candidates, importance = selector._ordered_candidates(X, y, "regression") gains = [importance[c] for c in candidates] assert gains == sorted(gains, reverse=True) + + +# --- tied / discrete features (issue #57) ------------------------------------ + + +def _discrete(levels, n=300, seed=0): + rng = np.random.default_rng(seed) + x = rng.integers(0, levels, size=n).astype(float) + return x.reshape(-1, 1), (x >= 2) + rng.normal(0, 0.1, size=n) + + +def _assert_distinct_interior(locations, X, count): + assert len(locations) == count + assert len(np.unique(locations)) == count + assert locations.min() > X.min() and locations.max() < X.max() + + +@pytest.mark.parametrize("selector_cls", [CARTLocationSelector, LightGBMLocationSelector]) +@pytest.mark.parametrize("levels", [2, 3, 5]) +def test_selector_tops_up_tied_features_with_distinct_interior_locations(selector_cls, levels): + """Regression guard for issue #57: the quantile top-up collapsed onto the tied + values and the range boundary, returning too few / duplicate locations.""" + if selector_cls is LightGBMLocationSelector: + pytest.importorskip("lightgbm") + X, y = _discrete(levels) + locations = selector_cls().select(X, y, task="regression", min_count=8, max_count=8) + _assert_distinct_interior(locations, X, 8) + + +def test_supplement_on_tied_data_respects_spacing_and_count(): + x = np.repeat([0.0, 1.0, 2.0], 50) + selector = CARTLocationSelector() + result = selector._supplement([0.5], x, 6) + assert len(result) == 6 + assert 0.5 in result + assert min(np.diff(result)) >= selector.min_location_spacing * 2.0 + assert min(result) > 0.0 and max(result) < 2.0 + + +def test_supplement_fills_dense_requests_beyond_the_spacing(): + x = np.linspace(0.0, 1.0, 1000) + result = CARTLocationSelector()._supplement([], x, 150) + assert len(result) == 150 + assert len(np.unique(result)) == 150 + + +def test_small_sample_fallback_on_tied_data_is_distinct_and_interior(): + X = np.array([[0.0], [0.0], [0.0], [1.0], [1.0]]) + locations = CARTLocationSelector().select(X, np.arange(5.0), min_count=3, max_count=3) + _assert_distinct_interior(locations, X, 3) + + +def test_constant_feature_fallback_keeps_the_requested_count(): + X = np.full((40, 1), 3.0) + locations = CARTLocationSelector().select(X, np.arange(40.0), min_count=4, max_count=4) + np.testing.assert_array_equal(locations, np.full(4, 3.0)) + + +def test_lightgbm_zero_split_sentinel_maps_to_the_bin_midpoint(): + pytest.importorskip("lightgbm") + rng = np.random.default_rng(0) + x = rng.integers(0, 5, size=300).astype(float) + y = (x >= 1) + rng.normal(0, 0.1, size=300) + locations = LightGBMLocationSelector().select(x.reshape(-1, 1), y, min_count=4, max_count=4) + # The 0-vs-positive split is reported by LightGBM as threshold 1e-35. + assert 0.5 in locations + assert (np.abs(locations) > 1e-30).all() diff --git a/tests/core/test_ple_selector.py b/tests/core/test_ple_selector.py index 498b9ac..497e09c 100644 --- a/tests/core/test_ple_selector.py +++ b/tests/core/test_ple_selector.py @@ -95,3 +95,15 @@ def test_lightgbm_selector_places_thresholds(data): t = PLETransformer(output_dim=5, task="regression", placement_strategy="lightgbm").fit(X, y) assert t.n_bins_per_feature_ == [5] assert np.all(np.diff(t.thresholds_[0]) > 0) + + +@pytest.mark.parametrize("levels", [2, 3, 4, 5, 6]) +def test_ple_keeps_output_dim_and_no_constant_columns_on_discrete_features(levels): + """Regression guard for issue #57: thresholds on the range boundary produced a + constant column and a narrower-than-requested output.""" + rng = np.random.default_rng(levels) + x = rng.integers(0, levels, size=400).astype(float).reshape(-1, 1) + y = np.sin(x[:, 0]) + rng.normal(0, 0.3, size=400) + out = PLETransformer(output_dim=7).fit(x, y).transform(x) + assert out.shape[1] == 7 + assert (out.std(axis=0) > 0).all() From c732a7ec115b495b01b4a12f15a81817353c635d Mon Sep 17 00:00:00 2001 From: ChrisW09 <50968720+ChrisW09@users.noreply.github.com> Date: Fri, 9 Oct 2026 13:18:47 +0200 Subject: [PATCH 06/18] fix(placement): label-encode classification targets for LightGBM placement LightGBMLocationSelector always trained the binary objective on the raw labels. Labels other than {0, 1} looked single-class (silently falling back to quantiles), multiclass targets only found the class-0 boundary, and string labels raised. Targets are now label-encoded and trained with the binary or multiclass objective as appropriate, matching the CART selector. Fixes #61 --- pretab/core/selectors.py | 16 ++++++++++++++-- tests/core/test_location_selectors.py | 26 ++++++++++++++++++++++++++ 2 files changed, 40 insertions(+), 2 deletions(-) diff --git a/pretab/core/selectors.py b/pretab/core/selectors.py index dba14d9..ffe86b7 100644 --- a/pretab/core/selectors.py +++ b/pretab/core/selectors.py @@ -355,9 +355,21 @@ def _import_lightgbm(): def _ordered_candidates(self, x_valid: np.ndarray, y_valid: np.ndarray, task: Task) -> tuple[list[float], object]: lgb = self._import_lightgbm() + if task == "regression": + objective = {"objective": "regression", "metric": "rmse"} + else: + # LightGBM needs integer class codes: its binary objective treats every + # label > 0 as positive and it cannot read string labels. + classes, y_valid = np.unique(y_valid, return_inverse=True) + if len(classes) < 2: + return [], None + if len(classes) == 2: + objective = {"objective": "binary", "metric": "binary_logloss"} + else: + objective = {"objective": "multiclass", "metric": "multi_logloss", "num_class": len(classes)} + params = { - "objective": "regression" if task == "regression" else "binary", - "metric": "rmse" if task == "regression" else "binary_logloss", + **objective, "num_leaves": 2**self.max_depth, "max_depth": self.max_depth, "learning_rate": self.learning_rate, diff --git a/tests/core/test_location_selectors.py b/tests/core/test_location_selectors.py index 49dcf01..e8d93ee 100644 --- a/tests/core/test_location_selectors.py +++ b/tests/core/test_location_selectors.py @@ -264,3 +264,29 @@ def test_lightgbm_zero_split_sentinel_maps_to_the_bin_midpoint(): # The 0-vs-positive split is reported by LightGBM as threshold 1e-35. assert 0.5 in locations assert (np.abs(locations) > 1e-30).all() + + +@pytest.mark.parametrize( + "make_target", + [ + pytest.param(lambda x: np.where(x < 0.9, 1, 2), id="labels-1-2"), + pytest.param(lambda x: np.digitize(x, [0.1, 0.9]), id="three-classes"), + pytest.param(lambda x: np.where(x < 0.9, "no", "yes"), id="string-labels"), + ], +) +def test_lightgbm_classification_places_a_location_at_the_class_boundary(make_target): + """Regression guard for issue #61: the binary objective on raw labels lost the + supervision for non-0/1, multiclass and string targets.""" + pytest.importorskip("lightgbm") + x = np.random.default_rng(0).uniform(0, 1, size=1000) + locations = LightGBMLocationSelector().select( + x.reshape(-1, 1), make_target(x), task="classification", min_count=2, max_count=2 + ) + assert np.abs(locations - 0.9).min() < 0.01, locations + + +def test_lightgbm_classification_with_a_single_class_falls_back(): + pytest.importorskip("lightgbm") + x = np.linspace(0, 1, 200).reshape(-1, 1) + locations = LightGBMLocationSelector().select(x, np.ones(200), task="classification", min_count=2, max_count=2) + assert len(locations) == 2 From da41ecad6d84d596bdfdf3183c61c9cd93e873d1 Mon Sep 17 00:00:00 2001 From: ChrisW09 <50968720+ChrisW09@users.noreply.github.com> Date: Fri, 9 Oct 2026 13:21:46 +0200 Subject: [PATCH 07/18] fix(splines): search the resolved knot window in target-aware spline placement The B/M/I, cubic-regression and natural-cubic splines asked their target-aware selector for a fixed 3-15 basis-function window and then reduced the result to the needed count with select_knots, which keeps evenly spaced indices of the sorted list. That positional trim routinely discarded the split where the target changes, and capped adaptive widths at 15 (14 / 12 for the cubic families). The spline adapter now accepts a per-call knot window; the transformers pass the window resolved from output_dim or [min_output_dim, max_output_dim], so the selector's importance ranking keeps the most informative splits, and a short set is only topped up. The adaptive tutorial output is updated accordingly. Fixes #69 --- docs/tutorials/adaptive_resolution.md | 13 +++--- pretab/expansion/spline/base.py | 7 +++- pretab/expansion/spline/mixins.py | 29 ++++++------- pretab/placement/adapters.py | 22 +++++++--- .../spline/test_spline_expansions.py | 42 +++++++++++++++++++ .../test_spline_placement_adapter.py | 15 +++++++ 6 files changed, 98 insertions(+), 30 deletions(-) diff --git a/docs/tutorials/adaptive_resolution.md b/docs/tutorials/adaptive_resolution.md index 9159e01..6f7870a 100644 --- a/docs/tutorials/adaptive_resolution.md +++ b/docs/tutorials/adaptive_resolution.md @@ -57,15 +57,16 @@ for name, y in [("simple", simple), ("wiggly", wiggly)]: ``` ```text -simple -> selected width 15 -wiggly -> selected width 15 +simple -> selected width 20 +wiggly -> selected width 20 ``` Both widths land inside the `[5, 20]` window without you having to guess a number up front. -The two happen to match here because the underlying CART selector's split count is governed -more by its own tree depth and minimum-samples settings than by how wiggly the signal looks; -with noisier or smaller data, or a narrower window, the two searches can land on different -widths. The bound is what you control directly, the exact count inside it is data-driven. +Here both reach the upper bound: with 3000 samples the CART selector finds more informative +splits than the window admits, so it keeps the most informative ones. Its split count is governed +more by its own tree depth and minimum-samples settings than by how wiggly the signal looks; with +smaller data or a wider window the two searches can land on different widths. The bound is what +you control directly, the exact count inside it is data-driven. ```{note} Fitting a target-aware transformer directly like this, outside a `Pipeline`, normally emits a diff --git a/pretab/expansion/spline/base.py b/pretab/expansion/spline/base.py index 394429f..a6cfa61 100644 --- a/pretab/expansion/spline/base.py +++ b/pretab/expansion/spline/base.py @@ -239,7 +239,12 @@ def _column_knots( ) internal_knots = self._adjust_internal_knots(x_valid, np.asarray(self.knot_locations), min_knots, max_knots) elif selector is not None: - selected = selector.get_knot_locations(x_valid.reshape(-1, 1), y_valid, task=self.task) + # Search exactly the window this feature needs, so the selector's own + # importance ranking picks the knots instead of a positional trim. + min_knots = min(min_knots, max_knots) + selected = selector.get_knot_locations( + x_valid.reshape(-1, 1), y_valid, task=self.task, min_knots=min_knots, max_knots=max_knots + ) internal_knots = self._adjust_internal_knots(x_valid, np.asarray(selected), min_knots, max_knots) else: n_internal = self._basis_to_knots(n_basis) diff --git a/pretab/expansion/spline/mixins.py b/pretab/expansion/spline/mixins.py index a73048c..f288697 100644 --- a/pretab/expansion/spline/mixins.py +++ b/pretab/expansion/spline/mixins.py @@ -119,13 +119,18 @@ def _place_interior_knots(self, x, y, n_interior, strategy, selector, task, min_ if selector is not None: if y is None: raise IncompatibleParamsError("A knot selector requires y during fit for target-aware knot placement.") - selected = np.asarray(selector.get_knot_locations(x.reshape(-1, 1), y, task=task), dtype=float) - x_min, x_max = x.min(), x.max() - selected = np.unique(selected[(selected > x_min) & (selected < x_max)]) if min_interior is None and max_interior is None: # Fixed (non-adaptive) selector path: force exactly ``n_interior`` # interior knots so the width stays ``output_dim``. min_interior = max_interior = n_interior + # Search exactly this window, so the selector's own importance ranking + # picks the knots instead of a positional trim afterwards. + selected = selector.get_knot_locations( + x.reshape(-1, 1), y, task=task, min_knots=min_interior, max_knots=max_interior + ) + selected = np.asarray(selected, dtype=float) + x_min, x_max = x.min(), x.max() + selected = np.unique(selected[(selected > x_min) & (selected < x_max)]) selected = self._clamp_interior_knots(x, selected, min_interior, max_interior, strategy) return selected return supplement_interior_knots(x, generate_internal_knots(x, n_interior, strategy), n_interior) @@ -134,26 +139,16 @@ def _clamp_interior_knots(self, x, knots, min_count, max_count, strategy): """Clamp a data-driven set of interior knots into ``[min_count, max_count]``. Down-samples with :func:`pretab.core.knots.select_knots` when there are too - many knots and supplements with quantile / uniform interior candidates when - there are too few. Endpoints are never added -- the result stays strictly - interior. + many knots and keeps every knot while topping up with quantile / uniform + interior candidates when there are too few. Endpoints are never added -- + the result stays strictly interior. """ x = np.asarray(x) knots = np.unique(np.sort(np.asarray(knots, dtype=float))) if max_count is not None and len(knots) > max_count: knots = select_knots(knots, max_count) if min_count is not None and len(knots) < min_count: - x_min, x_max = x.min(), x.max() - candidates = [ - knots, - generate_internal_knots(x, min_count, "quantile"), - generate_internal_knots(x, min_count, "uniform"), - ] - combined = np.unique(np.concatenate(candidates)) - combined = combined[(combined > x_min) & (combined < x_max)] - if len(combined) > min_count: - combined = select_knots(combined, min_count) - knots = combined + knots = supplement_interior_knots(x, knots, min_count) return knots def _adaptive_interior_bounds(self, output_dim, selector, *, floor, offset): diff --git a/pretab/placement/adapters.py b/pretab/placement/adapters.py index 15a0e65..0f2e816 100644 --- a/pretab/placement/adapters.py +++ b/pretab/placement/adapters.py @@ -33,9 +33,9 @@ "SplinePlacementAdapter", ] -# The spline knot selectors have always searched a fixed basis-function window, -# independent of the requested output_dim (the transformer clamps to output_dim -# afterwards). These reproduce ``CART/LightGBMKnotSelector``'s defaults. +# Default basis-function search window of a standalone adapter, reproducing the +# historical ``CART/LightGBMKnotSelector`` defaults. The spline transformers do +# not rely on it: they pass their own resolved knot window per call. _SPLINE_MIN_BASIS = 3 _SPLINE_MAX_BASIS = 15 # Historical default seed used by the spline knot selectors when random_state is @@ -101,14 +101,24 @@ def get_knot_locations( X: np.ndarray, y: np.ndarray | None = None, task: Task | None = "regression", + *, + min_knots: int | None = None, + max_knots: int | None = None, ) -> np.ndarray: - """Return sorted internal knot locations for a single feature.""" + """Return sorted internal knot locations for a single feature. + + ``min_knots`` / ``max_knots`` override the adapter's search window for this + call. The spline transformers pass the interior-knot window resolved from + their own ``output_dim`` (or adaptive bounds), so that the selector's + importance ranking -- not a positional down-sample afterwards -- decides + which knots are kept. + """ seed = self.random_state if self.random_state is not None else _SPLINE_DEFAULT_SEED strategy = create_placement_strategy( target_aware=True, placement_strategy=self.placement_strategy, - min_count=self.min_knots, - max_count=self.max_knots, + min_count=self.min_knots if min_knots is None else min_knots, + max_count=self.max_knots if max_knots is None else max_knots, task=task, random_state=seed, ) diff --git a/tests/expansion/spline/test_spline_expansions.py b/tests/expansion/spline/test_spline_expansions.py index 38ea023..3ab0b88 100644 --- a/tests/expansion/spline/test_spline_expansions.py +++ b/tests/expansion/spline/test_spline_expansions.py @@ -3,8 +3,10 @@ from pretab.transformers import ( BSplineTransformer, + CubicRegressionSplineTransformer, ISplineTransformer, MSplineTransformer, + NaturalCubicSplineTransformer, ) @@ -211,3 +213,43 @@ def test_spline_penalty_matrix_symmetric(data): P = transformer.get_penalty_matrix() assert P.shape[0] == P.shape[1] assert np.allclose(P, P.T, atol=1e-9) + + +# --- target-aware splines search their own knot window (issue #69) ------------ + + +def _step_data(seed=1): + rng = np.random.default_rng(seed) + x = rng.uniform(0, 10, 1000) + y = np.where(x > 5.0, 1.0, 0.0) + 0.3 * rng.normal(size=1000) + return x.reshape(-1, 1), y + + +@pytest.mark.parametrize("cls", [BSplineTransformer, MSplineTransformer, ISplineTransformer]) +@pytest.mark.parametrize("output_dim", [5, 6, 7]) +def test_target_aware_bmi_spline_keeps_the_dominant_split(cls, output_dim): + """Regression guard for issue #69: selector knots were trimmed by position, so + the split where the target changes was routinely discarded.""" + X, y = _step_data() + transformer = cls(output_dim=output_dim, target_aware=True, placement_strategy="cart").fit(X, y) + interior = transformer.knots_[0][transformer.degree + 1 : -(transformer.degree + 1)] + assert len(interior) == output_dim - transformer.degree - 1 + assert np.abs(interior - 5.0).min() < 0.1, interior + + +@pytest.mark.parametrize("cls", [CubicRegressionSplineTransformer, NaturalCubicSplineTransformer]) +def test_target_aware_cubic_families_keep_the_dominant_split(cls): + X, y = _step_data() + knots = cls(output_dim=4, target_aware=True, placement_strategy="cart").fit(X, y).knots_[0] + assert np.abs(knots - 5.0).min() < 0.1, knots + + +@pytest.mark.parametrize("cls", [BSplineTransformer, CubicRegressionSplineTransformer, NaturalCubicSplineTransformer]) +def test_adaptive_spline_width_is_not_capped_at_the_legacy_window(cls): + rng = np.random.default_rng(0) + x = rng.uniform(0, 10, 3000) + y = np.sin(2 * x) * 3 + rng.normal(0, 0.3, 3000) + transformer = cls( + adaptive=True, min_output_dim=5, max_output_dim=40, target_aware=True, placement_strategy="cart" + ).fit(x.reshape(-1, 1), y) + assert 15 < transformer.total_output_dim_ <= 40 diff --git a/tests/placement/test_spline_placement_adapter.py b/tests/placement/test_spline_placement_adapter.py index 5d24fa1..52a1350 100644 --- a/tests/placement/test_spline_placement_adapter.py +++ b/tests/placement/test_spline_placement_adapter.py @@ -99,3 +99,18 @@ def test_lightgbm_adapter_reproducible(data): a = SplinePlacementAdapter(placement_strategy="lightgbm", degree=3).get_knot_locations(X, y) b = SplinePlacementAdapter(placement_strategy="lightgbm", degree=3).get_knot_locations(X, y) np.testing.assert_array_equal(a, b) + + +def _step_data(seed=1): + rng = np.random.default_rng(seed) + x = rng.uniform(0, 10, 1000) + y = np.where(x > 5.0, 1.0, 0.0) + 0.3 * rng.normal(size=1000) + return x.reshape(-1, 1), y + + +def test_get_knot_locations_window_override(): + X, y = _step_data() + adapter = SplinePlacementAdapter(placement_strategy="cart", degree=3) + knots = adapter.get_knot_locations(X, y, min_knots=1, max_knots=1) + assert len(knots) == 1 + assert abs(knots[0] - 5.0) < 0.1 From 87d396921ef53cbb8c718f41075d6eb704939062 Mon Sep 17 00:00:00 2001 From: ChrisW09 <50968720+ChrisW09@users.noreply.github.com> Date: Fri, 9 Oct 2026 13:25:14 +0200 Subject: [PATCH 08/18] fix(preprocessor): route columns by label and keep ColumnTransformer step names valid Preprocessor passed each column label to scikit-learn's ColumnTransformer both as the column selector and inside the step name. ColumnTransformer reads an integer selector as a position, so integer labels that differ from their positions silently swapped columns (or crashed), and labels starting with '_' or containing '__' produced step names scikit-learn rejects. The ColumnTransformer is now fitted and applied on a frame whose labels are strings (a shallow relabelled copy only when needed), and columns are selected by that string label. Steps keep the num_