Add include_bias option to Polynomial basis function

sysidentpy/basis_function/_polynomial.py

@@ -36,6 +36,12 @@ class Polynomial(BaseBasisFunction):
     ----------
     degree : int (max_degree), default=2
         The maximum degree of the polynomial features.
+    include_bias : bool, default=True
+        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
     -----
@@ -47,8 +53,10 @@ class Polynomial(BaseBasisFunction):
     def __init__(
         self,
         degree: int = 2,
+        include_bias: bool = True,
     ):
         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] = {}
 
@@ -138,6 +146,7 @@ def fit(
         ----------
         data : ndarray of floats
             The lagged matrix built with respect to each lag and column.
+            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
@@ -148,16 +157,29 @@ 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
         -------
         psi = ndarray of floats
             The lagged matrix built in respect with each lag and column.
 
         """
-        # Create combinations of all columns based on its index
-        psi = self._evaluate_terms(data, predefined_regressors)
+        if self.include_bias:
+            psi = self._evaluate_terms(data, predefined_regressors)
+            return psi[max_lag:, :]
+
+        # 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:
+            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[max_lag:, :]
 
     def transform(

sysidentpy/basis_function/tests/test_polynomial.py

@@ -0,0 +1,71 @@
+import numpy as np
+from numpy.testing import assert_array_equal
+
+from sysidentpy.basis_function import Polynomial
+
+
+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_init_defaults():
+    p = Polynomial()
+    assert p.degree == 2
+    assert p.include_bias is True
+
+
+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_default_matches_original_behavior():
+    """include_bias=True must preserve pre-PR fit output exactly."""
+    p = Polynomial(degree=2)
+    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)
+
+
+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)
+
+    assert psi_with.shape[1] == psi_without.shape[1] + 1
+    assert_array_equal(psi_without, psi_with[:, 1:])
+
+
+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)
+
+    constant_columns = [
+        i for i in range(psi.shape[1]) if np.allclose(psi[:, i], psi[0, i])
+    ]
+    assert constant_columns == []
+
+
+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))
+
+
+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]])

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