From 8074127a1f653e8b3fbf57a75b32cf0053729554 Mon Sep 17 00:00:00 2001 From: June Kim Date: Mon, 11 May 2026 12:19:48 -0700 Subject: [PATCH 1/3] Add include_bias parameter to Polynomial basis function Implements the feature requested in issue #193 to allow users to optionally include a bias (intercept) term in the Polynomial basis function, matching the interface already present in newer basis functions like Legendre. Changes: - Add `include_bias` parameter to Polynomial.__init__ (default=False for backward compatibility) - Modify fit() to prepend a column of ones when include_bias=True - Handle predefined_regressors correctly when bias is included - Support both NumPy and Array API backends Tests: - Add comprehensive test suite in test_polynomial.py - Verify backward compatibility with existing behavior - Test bias inclusion/exclusion - Test interaction with predefined_regressors - Test shape changes Maintains backward compatibility: default behavior unchanged. --- sysidentpy/basis_function/_polynomial.py | 39 ++++++++- .../basis_function/tests/test_polynomial.py | 85 +++++++++++++++++++ 2 files changed, 122 insertions(+), 2 deletions(-) create mode 100644 sysidentpy/basis_function/tests/test_polynomial.py diff --git a/sysidentpy/basis_function/_polynomial.py b/sysidentpy/basis_function/_polynomial.py index dbeff48e..03804fa5 100644 --- a/sysidentpy/basis_function/_polynomial.py +++ b/sysidentpy/basis_function/_polynomial.py @@ -36,6 +36,9 @@ class Polynomial(BaseBasisFunction): ---------- degree : int (max_degree), default=2 The maximum degree of the polynomial features. + include_bias : bool, default=False + Whether to include the bias (constant) term in the output feature matrix. + When set to True, a column of ones is prepended to the feature matrix. Notes ----- @@ -47,8 +50,10 @@ class Polynomial(BaseBasisFunction): def __init__( self, degree: int = 2, + include_bias: bool = False, ): self.degree = degree + self.include_bias = include_bias # Cache combination indices per (n_features, degree) to avoid rebuilding self._combination_cache: Dict[Tuple[int, int], np.ndarray] = {} @@ -156,9 +161,39 @@ def fit( The lagged matrix built in respect with each lag and column. """ + xp = get_namespace(data) + # Create combinations of all columns based on its index - psi = self._evaluate_terms(data, predefined_regressors) - return psi[max_lag:, :] + # Note: predefined_regressors filtering happens inside _evaluate_terms + # only when include_bias is False. When include_bias is True, we need + # to handle it after adding the bias column. + if self.include_bias: + # Evaluate all terms first, then add bias + psi = self._evaluate_terms(data, predefined_regressors=None) + psi = psi[max_lag:, :] + + # Add bias column at the beginning + target_device = _device(data) if not _is_numpy_namespace(xp) else None + bias_column = _ones( + xp, + (psi.shape[0], 1), + dtype=psi.dtype, + target_device=target_device, + ) + psi = _column_stack(xp, [bias_column, psi]) + + # Apply predefined_regressors filtering after adding bias + if predefined_regressors is not None: + predefined_regressors = self._normalize_predefined_regressors( + predefined_regressors + ) + psi = psi[:, predefined_regressors] + else: + # Original behavior: no bias term + psi = self._evaluate_terms(data, predefined_regressors) + psi = psi[max_lag:, :] + + return psi def transform( self, diff --git a/sysidentpy/basis_function/tests/test_polynomial.py b/sysidentpy/basis_function/tests/test_polynomial.py new file mode 100644 index 00000000..8c40c84d --- /dev/null +++ b/sysidentpy/basis_function/tests/test_polynomial.py @@ -0,0 +1,85 @@ +import numpy as np +from sysidentpy.basis_function import Polynomial + + +def test_polynomial_init(): + """Test that Polynomial class initializes correctly.""" + p = Polynomial(degree=3, include_bias=False) + assert p.degree == 3 + assert not p.include_bias + + p_default = Polynomial() + assert p_default.degree == 2 + assert not p_default.include_bias # Default should be False for backward compatibility + + +def test_polynomial_fit(): + """Test that fit correctly generates the polynomial feature matrix.""" + p = Polynomial(degree=2) + data = np.random.rand(10, 3) # 10 samples, 3 features + max_lag = 2 + + transformed = p.fit(data, max_lag=max_lag) + + assert transformed.shape[0] == data.shape[0] - max_lag + assert transformed.shape[1] > 0 # Ensure non-empty feature matrix + + +def test_polynomial_transform(): + """Test that transform behaves identically to fit.""" + p = Polynomial(degree=2) + data = np.random.rand(10, 3) + + transformed_fit = p.fit(data) + transformed_transform = p.transform(data) + + np.testing.assert_array_equal(transformed_fit, transformed_transform) + + +def test_polynomial_include_bias(): + """Test that bias is included when requested.""" + p = Polynomial(degree=2, include_bias=True) + data = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) + max_lag = 1 + + transformed = p.fit(data, max_lag=max_lag) + # First column should be all ones (bias term) + assert np.all(transformed[:, 0] == 1) + + p_no_bias = Polynomial(degree=2, include_bias=False) + transformed_no_bias = p_no_bias.fit(data, max_lag=max_lag) + # Should not have a bias column (no column of all ones) + has_bias_column = any(np.all(transformed_no_bias[:, i] == 1) + for i in range(transformed_no_bias.shape[1])) + assert not has_bias_column + + +def test_polynomial_include_bias_shapes(): + """Test that include_bias affects the output shape correctly.""" + data = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12]]) + max_lag = 1 + + p_with_bias = Polynomial(degree=2, include_bias=True) + p_without_bias = Polynomial(degree=2, include_bias=False) + + transformed_with = p_with_bias.fit(data, max_lag=max_lag) + transformed_without = p_without_bias.fit(data, max_lag=max_lag) + + # With bias should have one more column + assert transformed_with.shape[1] == transformed_without.shape[1] + 1 + # Same number of rows + assert transformed_with.shape[0] == transformed_without.shape[0] + + +def test_polynomial_fit_predefined_regressors_with_bias(): + """Test fit with predefined regressors when include_bias=True.""" + p = Polynomial(degree=2, include_bias=True) + data = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) + max_lag = 1 + predefined_regressors = np.array([0, 2]) # Selecting bias and one other regressor + + transformed = p.fit(data, max_lag=max_lag, predefined_regressors=predefined_regressors) + + assert transformed.shape[1] == len(predefined_regressors) + # First selected regressor should be bias (all ones) + assert np.all(transformed[:, 0] == 1) From 76c4e7a8ba25e02ec3ddab72dcbdf04613e3c606 Mon Sep 17 00:00:00 2001 From: June Kim Date: Mon, 11 May 2026 12:57:24 -0700 Subject: [PATCH 2/3] fix: address review feedback on include_bias implementation Fixes four critical bugs identified in code review: 1. Double bias bug: build_lagged_matrix already prepends a bias column at index 0. The previous implementation added another bias column, creating rank deficiency. Now we strip the existing bias column before polynomial expansion, then optionally re-add it. 2. include_bias=False now actually excludes bias: Previously the (0,0,...,0) combination still generated a constant term even when include_bias=False. Now we remove the bias column before expansion, so all-zero combinations cannot occur. 3. regressor_space index mapping: By stripping the bias column first and optionally re-adding it after expansion, regressor indices now map correctly whether include_bias is True or False. 4. Default consistency: Changed default from False to True to match Legendre basis function behavior. Implementation follows Legendre's pattern: - Remove bias column (index 0) from build_lagged_matrix output - Generate polynomial terms from remaining features - Re-add bias column if include_bias=True - Apply predefined_regressors filtering after bias handling --- sysidentpy/basis_function/_polynomial.py | 43 +++++++++---------- .../basis_function/tests/test_polynomial.py | 2 +- 2 files changed, 22 insertions(+), 23 deletions(-) diff --git a/sysidentpy/basis_function/_polynomial.py b/sysidentpy/basis_function/_polynomial.py index 03804fa5..06735b21 100644 --- a/sysidentpy/basis_function/_polynomial.py +++ b/sysidentpy/basis_function/_polynomial.py @@ -36,9 +36,11 @@ class Polynomial(BaseBasisFunction): ---------- degree : int (max_degree), default=2 The maximum degree of the polynomial features. - include_bias : bool, default=False + include_bias : bool, default=True Whether to include the bias (constant) term in the output feature matrix. - When set to True, a column of ones is prepended to the feature matrix. + When set to True, the bias column from build_lagged_matrix is preserved. + When set to False, the bias column is removed and all-zero combinations + are excluded from the polynomial expansion. Notes ----- @@ -50,7 +52,7 @@ class Polynomial(BaseBasisFunction): def __init__( self, degree: int = 2, - include_bias: bool = False, + include_bias: bool = True, ): self.degree = degree self.include_bias = include_bias @@ -143,6 +145,7 @@ def fit( ---------- data : ndarray of floats The lagged matrix built with respect to each lag and column. + Expected to include a bias column at index 0 (from build_lagged_matrix). max_lag : int Target data used on training phase. ylag : ndarray of int @@ -163,16 +166,16 @@ def fit( """ xp = get_namespace(data) - # Create combinations of all columns based on its index - # Note: predefined_regressors filtering happens inside _evaluate_terms - # only when include_bias is False. When include_bias is True, we need - # to handle it after adding the bias column. - if self.include_bias: - # Evaluate all terms first, then add bias - psi = self._evaluate_terms(data, predefined_regressors=None) - psi = psi[max_lag:, :] + # Remove the bias column that build_lagged_matrix already added (column 0) + # and trim max_lag rows, following Legendre's pattern + data_no_bias = data[max_lag:, 1:] + + # Generate polynomial terms from the non-bias columns + psi = self._evaluate_terms(data_no_bias, predefined_regressors=None) - # Add bias column at the beginning + # If include_bias=True, prepend a bias column + # If include_bias=False, we're done (no bias, no all-zeros combinations) + if self.include_bias: target_device = _device(data) if not _is_numpy_namespace(xp) else None bias_column = _ones( xp, @@ -182,16 +185,12 @@ def fit( ) psi = _column_stack(xp, [bias_column, psi]) - # Apply predefined_regressors filtering after adding bias - if predefined_regressors is not None: - predefined_regressors = self._normalize_predefined_regressors( - predefined_regressors - ) - psi = psi[:, predefined_regressors] - else: - # Original behavior: no bias term - psi = self._evaluate_terms(data, predefined_regressors) - psi = psi[max_lag:, :] + # Apply predefined_regressors filtering after bias handling + if predefined_regressors is not None: + predefined_regressors = self._normalize_predefined_regressors( + predefined_regressors + ) + psi = psi[:, predefined_regressors] return psi diff --git a/sysidentpy/basis_function/tests/test_polynomial.py b/sysidentpy/basis_function/tests/test_polynomial.py index 8c40c84d..4a6f527e 100644 --- a/sysidentpy/basis_function/tests/test_polynomial.py +++ b/sysidentpy/basis_function/tests/test_polynomial.py @@ -10,7 +10,7 @@ def test_polynomial_init(): p_default = Polynomial() assert p_default.degree == 2 - assert not p_default.include_bias # Default should be False for backward compatibility + assert p_default.include_bias # Default should be True to match Legendre def test_polynomial_fit(): From 9241e253c243ec3449016c9f1ca0a03b5a5bc914 Mon Sep 17 00:00:00 2001 From: June Kim Date: Mon, 18 May 2026 10:31:50 -0700 Subject: [PATCH 3/3] fix: preserve default behavior, drop bias row from regressor_code Prior commit broke backward compatibility: with include_bias=True the output lost the cross-bias polynomial terms (combinations like (0,k)) because fit() sliced column 0 from the input before polynomial expansion. test_fit_polynomial in test_basis_functions.py caught this. Restore the original fit() path on default. When include_bias=False, drop the (0,0,...,0) pure-bias column emitted first by combinations_with_replacement, and drop the matching row from RegressorDictionary.regressor_space so regressor_code stays aligned with psi. Per maintainer guidance on PR #197. --- sysidentpy/basis_function/_polynomial.py | 50 +++----- .../basis_function/tests/test_polynomial.py | 114 ++++++++---------- sysidentpy/narmax_base.py | 8 ++ 3 files changed, 77 insertions(+), 95 deletions(-) diff --git a/sysidentpy/basis_function/_polynomial.py b/sysidentpy/basis_function/_polynomial.py index 06735b21..e9dee620 100644 --- a/sysidentpy/basis_function/_polynomial.py +++ b/sysidentpy/basis_function/_polynomial.py @@ -37,10 +37,11 @@ class Polynomial(BaseBasisFunction): degree : int (max_degree), default=2 The maximum degree of the polynomial features. include_bias : bool, default=True - Whether to include the bias (constant) term in the output feature matrix. - When set to True, the bias column from build_lagged_matrix is preserved. - When set to False, the bias column is removed and all-zero combinations - are excluded from the polynomial expansion. + Whether the constant (pure-bias) regressor is part of the candidate set. + When True, behavior matches prior releases: the (0,0,...,0) combination + produces a column of ones at index 0 of the output. When False, that + column is dropped from the candidate set; ``RegressorDictionary.regressor_space`` + drops the matching row so ``regressor_code`` and ``psi`` stay aligned. Notes ----- @@ -145,7 +146,7 @@ def fit( ---------- data : ndarray of floats The lagged matrix built with respect to each lag and column. - Expected to include a bias column at index 0 (from build_lagged_matrix). + Column 0 is the bias column added by ``build_input_output_matrix``. max_lag : int Target data used on training phase. ylag : ndarray of int @@ -156,7 +157,8 @@ def fit( The type of the model (NARMAX, NAR or NFIR). predefined_regressors : ndarray of int The index of the selected regressors by the Model Structure - Selection algorithm. + Selection algorithm. Indices address the candidate set after + ``include_bias`` has been applied. Returns ------- @@ -164,35 +166,21 @@ def fit( The lagged matrix built in respect with each lag and column. """ - xp = get_namespace(data) - - # Remove the bias column that build_lagged_matrix already added (column 0) - # and trim max_lag rows, following Legendre's pattern - data_no_bias = data[max_lag:, 1:] - - # Generate polynomial terms from the non-bias columns - psi = self._evaluate_terms(data_no_bias, predefined_regressors=None) - - # If include_bias=True, prepend a bias column - # If include_bias=False, we're done (no bias, no all-zeros combinations) if self.include_bias: - target_device = _device(data) if not _is_numpy_namespace(xp) else None - bias_column = _ones( - xp, - (psi.shape[0], 1), - dtype=psi.dtype, - target_device=target_device, - ) - psi = _column_stack(xp, [bias_column, psi]) + psi = self._evaluate_terms(data, predefined_regressors) + return psi[max_lag:, :] - # Apply predefined_regressors filtering after bias handling + # include_bias=False: drop the pure-bias combination (0,0,...,0), which + # combinations_with_replacement always emits at index 0. predefined_regressors + # addresses the bias-dropped index space, so shift by +1 before slicing combos. if predefined_regressors is not None: - predefined_regressors = self._normalize_predefined_regressors( - predefined_regressors - ) - psi = psi[:, predefined_regressors] + shifted = self._normalize_predefined_regressors(predefined_regressors) + 1 + psi = self._evaluate_terms(data, predefined_regressors=shifted) + else: + psi = self._evaluate_terms(data, predefined_regressors=None) + psi = psi[:, 1:] - return psi + return psi[max_lag:, :] def transform( self, diff --git a/sysidentpy/basis_function/tests/test_polynomial.py b/sysidentpy/basis_function/tests/test_polynomial.py index 4a6f527e..6418972c 100644 --- a/sysidentpy/basis_function/tests/test_polynomial.py +++ b/sysidentpy/basis_function/tests/test_polynomial.py @@ -1,85 +1,71 @@ import numpy as np -from sysidentpy.basis_function import Polynomial +from numpy.testing import assert_array_equal +from sysidentpy.basis_function import Polynomial -def test_polynomial_init(): - """Test that Polynomial class initializes correctly.""" - p = Polynomial(degree=3, include_bias=False) - assert p.degree == 3 - assert not p.include_bias - p_default = Polynomial() - assert p_default.degree == 2 - assert p_default.include_bias # Default should be True to match Legendre +def _lagged_with_bias(n_rows, n_features): + """Build a synthetic lagged matrix with the bias column at index 0.""" + rng = np.random.default_rng(0) + data = rng.uniform(size=(n_rows, n_features)) + data[:, 0] = 1.0 + return data -def test_polynomial_fit(): - """Test that fit correctly generates the polynomial feature matrix.""" - p = Polynomial(degree=2) - data = np.random.rand(10, 3) # 10 samples, 3 features - max_lag = 2 +def test_polynomial_init_defaults(): + p = Polynomial() + assert p.degree == 2 + assert p.include_bias is True - transformed = p.fit(data, max_lag=max_lag) - assert transformed.shape[0] == data.shape[0] - max_lag - assert transformed.shape[1] > 0 # Ensure non-empty feature matrix +def test_polynomial_init_explicit(): + p = Polynomial(degree=3, include_bias=False) + assert p.degree == 3 + assert p.include_bias is False -def test_polynomial_transform(): - """Test that transform behaves identically to fit.""" +def test_polynomial_default_matches_original_behavior(): + """include_bias=True must preserve pre-PR fit output exactly.""" p = Polynomial(degree=2) - data = np.random.rand(10, 3) - - transformed_fit = p.fit(data) - transformed_transform = p.transform(data) - - np.testing.assert_array_equal(transformed_fit, transformed_transform) - - -def test_polynomial_include_bias(): - """Test that bias is included when requested.""" - p = Polynomial(degree=2, include_bias=True) - data = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) - max_lag = 1 - - transformed = p.fit(data, max_lag=max_lag) - # First column should be all ones (bias term) - assert np.all(transformed[:, 0] == 1) + data = np.array([[1, 1, 1], [2, 3, 4], [3, 3, 3]]) + expected = np.array([[4, 6, 8, 9, 12, 16], [9, 9, 9, 9, 9, 9]]) + assert_array_equal(p.fit(data=data, max_lag=1), expected) - p_no_bias = Polynomial(degree=2, include_bias=False) - transformed_no_bias = p_no_bias.fit(data, max_lag=max_lag) - # Should not have a bias column (no column of all ones) - has_bias_column = any(np.all(transformed_no_bias[:, i] == 1) - for i in range(transformed_no_bias.shape[1])) - assert not has_bias_column +def test_polynomial_include_bias_false_drops_constant(): + """include_bias=False drops the (0,0,...,0) pure-bias column.""" + p_with = Polynomial(degree=2, include_bias=True) + p_without = Polynomial(degree=2, include_bias=False) + data = np.array([[1, 1, 1], [2, 3, 4], [3, 3, 3]]) + psi_with = p_with.fit(data=data, max_lag=1) + psi_without = p_without.fit(data=data, max_lag=1) -def test_polynomial_include_bias_shapes(): - """Test that include_bias affects the output shape correctly.""" - data = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12]]) - max_lag = 1 + assert psi_with.shape[1] == psi_without.shape[1] + 1 + assert_array_equal(psi_without, psi_with[:, 1:]) - p_with_bias = Polynomial(degree=2, include_bias=True) - p_without_bias = Polynomial(degree=2, include_bias=False) - transformed_with = p_with_bias.fit(data, max_lag=max_lag) - transformed_without = p_without_bias.fit(data, max_lag=max_lag) +def test_polynomial_include_bias_false_no_constant_column(): + """The remaining columns under include_bias=False are not pure constants.""" + p = Polynomial(degree=2, include_bias=False) + data = _lagged_with_bias(20, 4) + psi = p.fit(data=data, max_lag=2) - # With bias should have one more column - assert transformed_with.shape[1] == transformed_without.shape[1] + 1 - # Same number of rows - assert transformed_with.shape[0] == transformed_without.shape[0] + constant_columns = [ + i for i in range(psi.shape[1]) if np.allclose(psi[:, i], psi[0, i]) + ] + assert constant_columns == [] -def test_polynomial_fit_predefined_regressors_with_bias(): - """Test fit with predefined regressors when include_bias=True.""" - p = Polynomial(degree=2, include_bias=True) - data = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) - max_lag = 1 - predefined_regressors = np.array([0, 2]) # Selecting bias and one other regressor +def test_polynomial_transform_matches_fit(): + p = Polynomial(degree=2, include_bias=False) + data = _lagged_with_bias(10, 3) + assert_array_equal(p.fit(data=data, max_lag=1), p.transform(data=data, max_lag=1)) - transformed = p.fit(data, max_lag=max_lag, predefined_regressors=predefined_regressors) - assert transformed.shape[1] == len(predefined_regressors) - # First selected regressor should be bias (all ones) - assert np.all(transformed[:, 0] == 1) +def test_polynomial_predefined_regressors_index_into_candidate_set(): + """predefined_regressors indexes the candidate set after include_bias is applied.""" + p = Polynomial(degree=2, include_bias=False) + data = np.array([[1, 1, 1], [2, 3, 4], [3, 3, 3]]) + full = p.fit(data=data, max_lag=1) + picked = p.fit(data=data, max_lag=1, predefined_regressors=np.array([0, 2])) + assert_array_equal(picked, full[:, [0, 2]]) diff --git a/sysidentpy/narmax_base.py b/sysidentpy/narmax_base.py index f773c827..951a153e 100644 --- a/sysidentpy/narmax_base.py +++ b/sysidentpy/narmax_base.py @@ -217,6 +217,14 @@ def regressor_space(self, n_inputs: int) -> np.ndarray: regressor_code = np.array(regressor_code) regressor_code = regressor_code[:, regressor_code.shape[1] :: -1] + if ( + isinstance(self.basis_function, Polynomial) + and not getattr(self.basis_function, "include_bias", True) + ): + # combinations_with_replacement emits the (0,0,...,0) pure-bias tuple + # first; Polynomial.fit drops the matching psi column, so drop the row + # here to keep regressor_code aligned with psi. + regressor_code = regressor_code[1:] return regressor_code def _get_max_lag(self):