ferrolearn_preprocess/

polynomial_features.rs

//! Polynomial features: generate polynomial and interaction features.
//!
//! Given input features `[a, b]` and degree 2, this transformer generates:
//! - `[1, a, b, a², a·b, b²]` (default — `interaction_only = false`)
//! - `[1, a, b, a·b]` (with `interaction_only = true`)
//!
//! With `include_bias = false`, the constant column `1` is omitted.
//!
//! This transformer is stateless for value generation — call
//! [`Transform::transform`] directly. For scikit-learn API parity it ALSO
//! supports the stateful [`Fit`](ferrolearn_core::traits::Fit) →
//! [`FittedPolynomialFeatures`] path, which records `n_features_in_` /
//! `n_output_features_` / `powers_` and validates the input feature count in
//! `transform`; the fitted type's `transform` REUSES the very same value math as
//! the stateless path, so both paths are bit-identical to each other. That value
//! math (FIXED #2210) reproduces scikit-learn's INCREMENTAL degree-by-degree
//! column build (`_polynomial.py:525-564`) — each degree-`d` column is a stored
//! degree-`(d-1)` column times a single feature, in sklearn's `np.multiply` arg
//! order — so the polynomial VALUES are now bit-identical to the sklearn oracle
//! at EVERY degree, including degree≥3 (the prior independent per-column left-fold
//! diverged by 1 ULP on degree≥3 terms because FP multiplication is
//! non-associative).
//!
//! # `## REQ status`
//!
//! Binary (R-DEFER-2), translating `sklearn/preprocessing/_polynomial.py` (`class
//! PolynomialFeatures` `:99`). Design doc: `.design/preprocess/polynomial_features.md`. Expected
//! values from the live sklearn 1.5.2 oracle (R-CHAR-3). Consumers: `PipelineTransformer` impl +
//! crate re-export (`lib.rs`, grandfathered S5) + the `FittedPolynomialFeatures::transform` fitted
//! path. HONEST (R-HONEST-3): a dense int-degree transformer — the polynomial VALUES + column ORDER
//! match sklearn exactly, and the stateful `fit` → `FittedPolynomialFeatures` surface
//! (`n_features_in_`/`n_output_features_`/`powers_` + the transform feature-count check) SHIPS
//! (REQ-4/REQ-5), and the `ferrolearn-python` PyO3 binding (`_RsPolynomialFeatures` + the
//! `ferrolearn.PolynomialFeatures` stateful wrapper surfacing the fitted attributes) NOW SHIPS
//! (REQ-10); `get_feature_names_out`, degree-tuple, `order` (C/F memory layout), sparse, full-ctor,
//! ferray stay NOT-STARTED.
//!
//! | REQ | Status | Evidence |
//! |---|---|---|
//! | REQ-1 (int-degree dense values + column order) | SHIPPED | Column order: `fn feature_combinations` (bias=empty combo, DFS interaction_only strictly-increasing / else non-decreasing, sort by degree-then-lex) reproduces sklearn `_combinations` itertools order (`_polynomial.py:209-220`) column-for-column (this drives `powers_`/`n_output_features_`). Value math: `fn generate_poly_features` builds the columns DEGREE-BY-DEGREE via sklearn's incremental block recurrence (`_polynomial.py:525-564`) — each degree-`d` column = a stored degree-`(d-1)` column × a single feature in sklearn's `np.multiply` arg order — so values are bit-identical to the live oracle at EVERY degree. FIXED #2210: the prior independent per-column left-fold (`combo.iter().fold(1, |acc,j| acc*x[j])`) diverged by 1 ULP on degree≥3 terms (FP multiplication is non-associative); the incremental build now matches bit-for-bit (pin `divergence_polynomial_higher_degree_term_ulp_vs_sklearn`, x1*x2*x3 == `0x1.7b645a1cac082p+2`). Green guards (`tests/divergence_polynomial_features.rs`): 2-feat default `[1,2,3,4,6,9]`, 3-feat default `[1,2,3,5,4,6,10,9,15,25]`, 3-feat interaction `[1,2,3,5,6,10,15]`, deg-3 interaction no-bias `[2,3,5,6,10,15,30]`, multi-row. Consumers: `FittedPipelineTransformer` + re-export. |
//! | REQ-8 (input validation per check_array) | SHIPPED | FIXED #1180. `transform` guards (sklearn order) zero-samples → `InsufficientSamples`, zero-features → `InvalidParameter`, non-finite NaN/±inf → `InvalidParameter` — matching sklearn `transform` → `_validate_data` (`_polynomial.py:433-435`). Mirrors converged binarizer/normalizer. Critic two-round CLEAN: 13 tests incl. finite input (1e308) with overflowing product correctly ACCEPTED (input-only validation). |
//! | REQ-2 (degree tuple / min_degree) | NOT-STARTED | open prereq blocker #1181. Single `usize` degree, always starts at degree 1 (sklearn `:334-360`). |
//! | REQ-3 (order C/F param) | NOT-STARTED | open prereq blocker #1182. No `order` (sklearn `:201`,`:132-134`). |
//! | REQ-4 (stateful fit + n_features_in_/n_output_features_ + count check) | SHIPPED | FIXED #1183. `impl Fit<Array2<F>, ()> for PolynomialFeatures<F>` (`fit`): runs the SAME REQ-8 `validate_poly_input` guard (zero-samples → `InsufficientSamples`, zero-features/non-finite → `InvalidParameter`, sklearn `_validate_data` `_polynomial.py:323`), enumerates `feature_combinations(x.ncols())` (REUSED, not reimplemented), records `n_features_in_ = x.ncols()` and `n_output_features_ = combinations.len()` (== `_num_combinations`, sklearn `:362`), returns `FittedPolynomialFeatures`. `FittedPolynomialFeatures::transform` validates X (REQ-8) FIRST, generates via the SHARED `generate_poly_features` value math (delegated — byte-identical to the stateless path), THEN checks `x.ncols() != n_features_in_` → `ShapeMismatch` (sklearn `X has N features, but ... expecting M`, `:402-435`). The X-validation-before-n_features ORDER matches sklearn `_validate_data(reset=False)` → `check_array` before `_check_n_features` (#2207). Accessors `n_features_in()`/`n_output_features()`. Live-oracle tests: `fit_n_features_in_matches_ncols`, `fit_n_output_features_*`, `fitted_transform_matches_stateless_and_sklearn`, `fitted_transform_more/fewer_features_shape_mismatch`, `fitted_transform_nan_validation_before_n_features`. Consumers: `FittedPolynomialFeatures::transform` (the fitted path) + crate re-export `pub use polynomial_features::{PolynomialFeatures, FittedPolynomialFeatures}` (`lib.rs`). |
//! | REQ-5 (powers_ attribute) | SHIPPED | FIXED #1184. `FittedPolynomialFeatures<F>` carries `powers_: Array2<usize>` of shape `(n_output_features_, n_features_in_)` (matching `n_features_in_`'s `usize` storage idiom; sklearn `powers_` is int), built in `fit` from the SAME `feature_combinations` the value math uses (so row order == output column order): each combination → a `bincount` row counting occurrences of each feature index `0..n_features_in_`, the bias empty-combo → an all-zeros row (sklearn `np.vstack([np.bincount(c, minlength=n_features_in_) for c in combinations])`, `_polynomial.py:262-264`). Accessor `powers(&self) -> &Array2<usize>`. Live-oracle tests: `fit_powers_2feat/3feat_*` EXACT-match sklearn `.powers_` for default / `interaction_only` / `include_bias=False`, asserting shape + every entry. Consumer: the `powers_` accessor on the re-exported `FittedPolynomialFeatures`. |
//! | REQ-6 (get_feature_names_out) | NOT-STARTED | open prereq blocker #1185. None (sklearn `:266-303`). Depends on REQ-5. |
//! | REQ-7 (sparse CSR/CSC) | NOT-STARTED | open prereq blocker #1186. Dense-only (sklearn `:402-`,`:38-96`). |
//! | REQ-9 (full ctor + _parameter_constraints) | NOT-STARTED | open prereq blocker #1187. Positional `new`; validates only the `degree==0 && !include_bias` empty-output case (FIXED #2216: sklearn `_polynomial.py:326-330` rejects ONLY that case — `degree==0 && include_bias` is valid, emitting just the bias column). No keyword-only `*`-args, no `order`, no full `_parameter_constraints` (sklearn `:194-207`). |
//! | REQ-10 (PyO3 binding) | SHIPPED | FIXED #1188. `_RsPolynomialFeatures` (hand `#[pyclass]`, `ferrolearn-python/src/extras.rs`) over `ferrolearn_preprocess::{PolynomialFeatures, FittedPolynomialFeatures}<f64>` exposes `fit(x)` (ctor `PolynomialFeatures::new(degree, interaction_only, include_bias)`, `degree==0 && !include_bias`/NaN/±inf → `PyValueError`; FIXED #2216 `degree==0 && include_bias` is valid → bias-only column), `transform(x)` (`FittedPolynomialFeatures::transform` — wrong col count → `ShapeMismatch` → `PyValueError`), and `#[getter]`s `n_features_in_` (sklearn `_polynomial.py:323`), `n_output_features_` (`:362`), `powers_` (`:250-264`, `usize` `Array2` → numpy i64). Registered in `lib.rs`. The `_extras.py::PolynomialFeatures(_TransformerWrapper)` wrapper mirrors sklearn's ctor `(degree=2, *, interaction_only=False, include_bias=True, order="C")` (`_polynomial.py:201`, `degree` positional-or-keyword), is STATEFUL (overridden `fit` builds/fits `_rs`; pre-fit `transform`/attribute access → `NotFittedError` via `check_is_fitted`), does the FLOAT-only dtype cast-back (`check_array(dtype=FLOAT_DTYPES)`, `:432`: int→f64, float32→float32, float64→float64), and surfaces `n_features_in_`/`n_output_features_`/`powers_` as `@property`s reading from `_rs`. `order` is validated against sklearn's `StrOptions({"C","F"})` (bad string → `ValueError`); `'F'` is accepted-but-documented (the returned array is always C-contiguous — F memory-layout is REQ-3 #1182 NOT-STARTED — but the VALUES match sklearn's 'F' output, which is a layout-only knob). #2215-CLASS f64-ABI CAVEAT: for FLOAT32 input the polynomial products are computed in float64 (the f64 ABI) then cast back to float32 (sklearn computes them in float32), so the cast-back float32 VALUES can diverge from sklearn's by ~1e-6 on high-degree terms (the same reduced-precision caveat class as #2215/#2205) — float64 input is bit-exact; the float32 value pin uses a TOLERANT `atol=1e-5`. Exported in `__init__.py`. Non-test consumer: `ferrolearn.PolynomialFeatures(...).fit_transform(X)` + `lib.rs` `add_class` + `__init__.py` re-export. Live oracle: `tests/divergence_polynomial_features_py.py` (33 pass, sklearn 1.5.2) — fit_transform VALUES (d∈{2,3} × default/interaction_only/no-bias, float64 bit-exact), `powers_` int-exact, `n_output_features_`/`n_features_in_`, pre-fit `NotFittedError`, wrong-n_features `ValueError`, dtype int→f64/float32-tolerant/float64-bit-exact, get_params keys=={degree,interaction_only,include_bias,order}/defaults/clone/set_params, pipeline, order='F' values + bad-order `ValueError`. `get_feature_names_out` stays NOT-STARTED (REQ-6 #1185). |
//! | REQ-11 (ferray substrate) | NOT-STARTED | open prereq blocker #1189. `ndarray`+`num_traits`, not `ferray-core` (R-SUBSTRATE-1/2). |

use ferrolearn_core::error::FerroError;
use ferrolearn_core::pipeline::{FittedPipelineTransformer, PipelineTransformer};
use ferrolearn_core::traits::{Fit, Transform};
use ndarray::{Array1, Array2};
use num_traits::Float;

// ---------------------------------------------------------------------------
// PolynomialFeatures
// ---------------------------------------------------------------------------

/// A stateless polynomial feature generator.
///
/// Generates all polynomial combinations of the input features up to the
/// specified degree.
///
/// # Configuration
///
/// - `degree`: maximum polynomial degree (default `2`).
/// - `interaction_only`: if `true`, only cross-product terms are generated
///   (no pure powers like `a²`). Default `false`.
/// - `include_bias`: if `true`, a constant column of ones is prepended.
///   Default `true`.
///
/// # Examples
///
/// ```
/// use ferrolearn_preprocess::polynomial_features::PolynomialFeatures;
/// use ferrolearn_core::traits::Transform;
/// use ndarray::array;
///
/// let poly = PolynomialFeatures::<f64>::new(2, false, true).unwrap();
/// let x = array![[2.0, 3.0]];
/// let out = poly.transform(&x).unwrap();
/// // out = [[1, 2, 3, 4, 6, 9]]
/// ```
#[derive(Debug, Clone)]
pub struct PolynomialFeatures<F> {
    /// Maximum polynomial degree.
    pub(crate) degree: usize,
    /// If `true`, only interaction terms are produced (no pure powers).
    pub(crate) interaction_only: bool,
    /// If `true`, prepend a bias (constant ones) column.
    pub(crate) include_bias: bool,
    _marker: std::marker::PhantomData<F>,
}

impl<F: Float + Send + Sync + 'static> PolynomialFeatures<F> {
    /// Create a new `PolynomialFeatures` transformer.
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::InvalidParameter`] if `degree == 0 && !include_bias`.
    /// `degree == 0` is ALLOWED when `include_bias == true` — the output is a
    /// single all-ones bias column (`powers_ == [[0,..,0]]`,
    /// `n_output_features_ == 1`), mirroring scikit-learn, which only rejects the
    /// `degree == 0 && include_bias == False` case
    /// (`sklearn/preprocessing/_polynomial.py:326-330`: "Setting degree to zero
    /// and include_bias to False would result in an empty output array.").
    pub fn new(
        degree: usize,
        interaction_only: bool,
        include_bias: bool,
    ) -> Result<Self, FerroError> {
        // sklearn `_polynomial.py:326-330`: ONLY reject degree==0 when
        // include_bias is False (an empty output array); degree==0 with
        // include_bias=True is valid and emits just the bias column.
        if degree == 0 && !include_bias {
            return Err(FerroError::InvalidParameter {
                name: "degree".into(),
                reason: "Setting degree to zero and include_bias to False would \
                         result in an empty output array."
                    .into(),
            });
        }
        Ok(Self {
            degree,
            interaction_only,
            include_bias,
            _marker: std::marker::PhantomData,
        })
    }

    /// Create a `PolynomialFeatures` with default settings:
    /// degree=2, interaction_only=false, include_bias=true.
    #[must_use]
    pub fn default_config() -> Self {
        Self {
            degree: 2,
            interaction_only: false,
            include_bias: true,
            _marker: std::marker::PhantomData,
        }
    }

    /// Return the configured degree.
    #[must_use]
    pub fn degree(&self) -> usize {
        self.degree
    }

    /// Return whether only interaction terms are generated.
    #[must_use]
    pub fn interaction_only(&self) -> bool {
        self.interaction_only
    }

    /// Return whether a bias column is included.
    #[must_use]
    pub fn include_bias(&self) -> bool {
        self.include_bias
    }

    /// Generate all combinations (with repetition unless `interaction_only`)
    /// of feature indices up to `degree`.
    ///
    /// Returns a list of index-tuples, where each tuple specifies which
    /// feature indices to multiply together to produce one output column.
    fn feature_combinations(&self, n_features: usize) -> Vec<Vec<usize>> {
        let mut combos: Vec<Vec<usize>> = Vec::new();

        // Bias term: empty product = 1
        if self.include_bias {
            combos.push(vec![]);
        }

        // Generate combinations of degrees 1..=self.degree
        let mut stack: Vec<(Vec<usize>, usize)> = Vec::new();

        // Start with each feature at degree 1. For `degree == 0` there are NO
        // degree>=1 combinations (sklearn `_combinations` ranges over
        // `range(max(1, min_degree), max_degree + 1)` = `range(1, 1)` = empty,
        // `_polynomial.py:215-217`), so only the bias empty-combo (if
        // `include_bias`) is emitted.
        if self.degree >= 1 {
            for i in 0..n_features {
                stack.push((vec![i], i));
            }
        }

        while let Some((combo, last_idx)) = stack.pop() {
            combos.push(combo.clone());

            if combo.len() < self.degree {
                // Extend with another feature
                let start = if self.interaction_only {
                    // Strictly increasing indices — no repeated features
                    last_idx + 1
                } else {
                    // Non-decreasing indices — allows repeated features (pure powers)
                    last_idx
                };
                for i in start..n_features {
                    let mut new_combo = combo.clone();
                    new_combo.push(i);
                    stack.push((new_combo, i));
                }
            }
        }

        // Sort: bias first, then by combo length, then lexicographically
        combos.sort_by(|a, b| a.len().cmp(&b.len()).then_with(|| a.cmp(b)));

        combos
    }
}

/// Run the shared `check_array` input validation (REQ-8) used by the stateless
/// [`Transform::transform`], the stateful [`Fit::fit`], and the
/// [`FittedPolynomialFeatures`] transform path, in sklearn's `check_array`
/// order: zero-samples → zero-features → non-finite
/// (`sklearn/utils/validation.py:1084`, `:1093`, `:1063`). Mirrors sklearn
/// `PolynomialFeatures.fit`/`.transform` → `_validate_data`
/// (`_polynomial.py:323`, `:433-435`), whose default `force_all_finite=True`
/// REJECTS NaN/±inf.
///
/// `context` names the calling site for diagnostics (e.g.
/// `"PolynomialFeatures::transform"` vs `"PolynomialFeatures::fit"`).
fn validate_poly_input<F: Float>(x: &Array2<F>, context: &str) -> Result<(), FerroError> {
    if x.nrows() == 0 {
        return Err(FerroError::InsufficientSamples {
            required: 1,
            actual: 0,
            context: context.into(),
        });
    }
    if x.ncols() == 0 {
        return Err(FerroError::InvalidParameter {
            name: "X".to_string(),
            reason: "Found array with 0 feature(s); a minimum of 1 is required \
                     by PolynomialFeatures"
                .to_string(),
        });
    }
    if x.iter().any(|v| !v.is_finite()) {
        return Err(FerroError::InvalidParameter {
            name: "X".to_string(),
            reason: "Input X contains non-finite values (NaN or infinity); \
                     PolynomialFeatures requires all-finite input"
                .to_string(),
        });
    }
    Ok(())
}

/// Generate the dense polynomial feature matrix (REQ-1 value math), shared by
/// the stateless [`Transform::transform`] and the stateful
/// [`FittedPolynomialFeatures`] transform path so the two are byte-identical.
///
/// Builds the output columns DEGREE-BY-DEGREE using scikit-learn's incremental
/// block recurrence (`_polynomial.py:525-564`), NOT an independent per-column
/// left-fold: the degree-1 block is `X` (copied in feature order); each degree-`d`
/// block is built feature-by-feature, where for `feature_idx` the new columns are
/// the just-built degree-`(d-1)` columns in slice `[start, end)` each multiplied
/// by `X[:, feature_idx]` — i.e. `previous_degree_column * X[:, feature_idx]` in
/// THAT multiply order, matching `np.multiply(XP[:, start:end], X[:, feature_idx],
/// casting="no")` bit-for-bit (FP multiplication is non-associative, so the
/// association order is load-bearing for the last ULP on degree≥3 terms). The
/// resulting column order is exactly sklearn's (bias, then degree-1, then each
/// degree-`d` block by `feature_idx`), which equals the [`feature_combinations`]
/// order the `powers_`/REQ-1 guards pin. The bias column (if `include_bias`) is
/// `1`. The input is assumed already validated by [`validate_poly_input`].
fn generate_poly_features<F: Float>(
    x: &Array2<F>,
    degree: usize,
    interaction_only: bool,
    include_bias: bool,
) -> Array2<F> {
    let n_samples = x.nrows();
    let n_features = x.ncols();
    let n_out = num_output_columns(n_features, degree, interaction_only, include_bias);
    let mut out = Array2::zeros((n_samples, n_out));

    // degree 0 (bias) term
    let mut current_col = 0usize;
    if include_bias {
        for i in 0..n_samples {
            out[[i, 0]] = F::one();
        }
        current_col = 1;
    }

    if n_features == 0 {
        return out;
    }

    if degree == 0 {
        // sklearn `_polynomial.py:532-533`: for `_max_degree == 0` the build
        // returns right after the bias column — NO degree-1 block, NO higher
        // blocks. degree==0 here implies include_bias (new() rejects
        // degree==0 && !include_bias), so `out` already holds the single
        // all-ones bias column.
        return out;
    }

    // degree 1 block: copy X column-for-column.
    // `index[j]` is the output column where degree-(d-1) block for feature `j`
    // begins; the trailing `index[n_features]` marks one past the block end.
    let mut index: Vec<usize> = Vec::with_capacity(n_features + 1);
    for j in 0..n_features {
        index.push(current_col);
        for i in 0..n_samples {
            out[[i, current_col]] = x[[i, j]];
        }
        current_col += 1;
    }
    index.push(current_col);

    // degree >= 2 blocks, built incrementally from the previous degree block.
    for _ in 2..=degree {
        let mut new_index: Vec<usize> = Vec::with_capacity(n_features + 1);
        // sklearn `end = index[-1]` (`_polynomial.py:544`): the LAST element of
        // `index`, NOT `index[n_features]`. When `interaction_only` runs out of
        // combinations at a high degree the inner loop `break`s early, leaving
        // `index` (the previous `new_index`) with FEWER than `n_features+1`
        // entries — so indexing `[n_features]` would panic (#2211, R-CODE-2). The
        // last element is the correct degree-(d-1) block end either way.
        let end = index.last().copied().unwrap_or(0);
        for feature_idx in 0..n_features {
            let mut start = index[feature_idx];
            new_index.push(current_col);
            if interaction_only {
                start += index[feature_idx + 1] - index[feature_idx];
            }
            if end <= start {
                // next_col <= current_col → no new columns for this feature.
                break;
            }
            // new column = (stored degree-(d-1) column) * X[:, feature_idx],
            // in that np.multiply arg order (prev FIRST, single feature SECOND).
            for src_col in start..end {
                for i in 0..n_samples {
                    out[[i, current_col]] = out[[i, src_col]] * x[[i, feature_idx]];
                }
                current_col += 1;
            }
        }
        new_index.push(current_col);
        index = new_index;
    }

    out
}

/// Number of output columns for the incremental build (== number of
/// [`feature_combinations`] entries / `n_output_features_`): the count of
/// degree-`1..=degree` combinations (with repetition unless `interaction_only`)
/// over `n_features` features, plus one for the bias if `include_bias`.
fn num_output_columns(
    n_features: usize,
    degree: usize,
    interaction_only: bool,
    include_bias: bool,
) -> usize {
    let mut total = usize::from(include_bias);
    if n_features == 0 {
        return total;
    }
    if degree == 0 {
        // sklearn `_combinations` emits no degree>=1 terms for degree==0
        // (`range(max(1,0), 0+1)` = `range(1,1)` = empty); the only output
        // column is the bias (`_polynomial.py:215-219`). degree==0 reaching here
        // implies include_bias (new() rejects degree==0 && !include_bias).
        return total;
    }
    // Count combinations per degree by mirroring the incremental block sizes.
    // degree-1 block has `n_features` columns; subsequent block sizes follow the
    // same start/end recurrence as the value build.
    let mut index: Vec<usize> = (0..=n_features).collect();
    total += n_features;
    let mut current = n_features;
    for _ in 2..=degree {
        let mut new_index: Vec<usize> = Vec::with_capacity(n_features + 1);
        // sklearn `index[-1]` (the last element), not `index[n_features]` —
        // robust to the `interaction_only` high-degree early break (#2211).
        let end = index.last().copied().unwrap_or(0);
        for feature_idx in 0..n_features {
            let mut start = index[feature_idx];
            new_index.push(current);
            if interaction_only {
                start += index[feature_idx + 1] - index[feature_idx];
            }
            if end <= start {
                break;
            }
            current += end - start;
            total += end - start;
        }
        new_index.push(current);
        index = new_index;
    }
    total
}

/// Build the `powers_` exponent matrix from a precomputed combination list
/// (REQ-5). Returns an `Array2<usize>` of shape `(combos.len(), n_features)`
/// where `powers_[k, j]` is the number of occurrences of input feature `j` in
/// combination `k` (its exponent) — sklearn's
/// `np.vstack([np.bincount(c, minlength=n_features_in_) for c in combinations])`
/// (`_polynomial.py:262-264`). The bias empty-combo → an all-zeros row. Built
/// from the SAME combinations the value math uses, so the row order equals the
/// output column order.
fn build_powers(combos: &[Vec<usize>], n_features: usize) -> Array2<usize> {
    let mut powers = Array2::<usize>::zeros((combos.len(), n_features));
    for (k, combo) in combos.iter().enumerate() {
        for &j in combo {
            powers[[k, j]] += 1;
        }
    }
    powers
}

impl<F: Float + Send + Sync + 'static> Default for PolynomialFeatures<F> {
    fn default() -> Self {
        Self::default_config()
    }
}

// ---------------------------------------------------------------------------
// Trait implementations
// ---------------------------------------------------------------------------

impl<F: Float + Send + Sync + 'static> Transform<Array2<F>> for PolynomialFeatures<F> {
    type Output = Array2<F>;
    type Error = FerroError;

    /// Generate polynomial and interaction features.
    ///
    /// # Errors
    ///
    /// Mirrors scikit-learn's `PolynomialFeatures.transform`
    /// (`sklearn/preprocessing/_polynomial.py:433`, `self._validate_data(...)` ->
    /// `check_array`), validating in `check_array` order
    /// (samples -> features -> finite):
    ///
    /// - Returns [`FerroError::InsufficientSamples`] if `x` has zero rows
    ///   (`check_array` `ensure_min_samples=1`; sklearn raises `ValueError:
    ///   Found array with 0 sample(s) ... a minimum of 1 is required`).
    /// - Returns [`FerroError::InvalidParameter`] if `x` has zero columns
    ///   (`check_array` `ensure_min_features=1`; sklearn raises `ValueError:
    ///   Found array with 0 feature(s) ... a minimum of 1 is required`).
    /// - Returns [`FerroError::InvalidParameter`] if `x` contains any non-finite
    ///   value (NaN, +inf, or -inf; `check_array` `force_all_finite=True`;
    ///   sklearn raises `ValueError: Input X contains NaN` / `infinity ...`).
    fn transform(&self, x: &Array2<F>) -> Result<Array2<F>, FerroError> {
        validate_poly_input(x, "PolynomialFeatures::transform")?;
        Ok(generate_poly_features(
            x,
            self.degree,
            self.interaction_only,
            self.include_bias,
        ))
    }
}

// ---------------------------------------------------------------------------
// FittedPolynomialFeatures (sklearn stateful `fit` -> fitted estimator path)
// ---------------------------------------------------------------------------

/// A fitted [`PolynomialFeatures`].
///
/// `PolynomialFeatures.fit` (sklearn `_polynomial.py:306-400`, "Compute number
/// of output features") learns NO numeric statistics — it records only the
/// shape metadata derivable from `n_features_in_` and the configured params:
/// `n_features_in_`, `n_output_features_` (the number of output columns), and
/// the `powers_` exponent matrix. The fitted type's [`Transform::transform`]
/// REUSES the very same combination enumeration ([`PolynomialFeatures::
/// feature_combinations`]) and value math ([`generate_poly_features`]) as the
/// stateless [`PolynomialFeatures`] path, so the two paths are bit-identical;
/// it additionally enforces the transform-time feature-count check against
/// `n_features_in_` (sklearn `X has N features, but ... expecting M`).
#[derive(Debug, Clone)]
pub struct FittedPolynomialFeatures<F> {
    /// Maximum polynomial degree (carried from the unfitted [`PolynomialFeatures`]).
    pub(crate) degree: usize,
    /// If `true`, only interaction terms are produced (no pure powers).
    pub(crate) interaction_only: bool,
    /// If `true`, a bias (constant ones) column is prepended.
    pub(crate) include_bias: bool,
    /// Number of features (columns) seen during [`Fit::fit`] — sklearn's
    /// `n_features_in_` (`_polynomial.py:323`, set by `_validate_data`).
    pub(crate) n_features_in_: usize,
    /// Total number of polynomial output columns — sklearn's
    /// `n_output_features_` (`_polynomial.py:362`, == `_num_combinations`).
    pub(crate) n_output_features_: usize,
    /// The exponent matrix of shape `(n_output_features_, n_features_in_)`:
    /// `powers_[i, j]` is the exponent of input feature `j` in output feature
    /// `i` — sklearn's `powers_` (`_polynomial.py:250-264`). The bias empty-combo
    /// row is all zeros.
    pub(crate) powers_: Array2<usize>,
    _marker: std::marker::PhantomData<F>,
}

impl<F: Float + Send + Sync + 'static> FittedPolynomialFeatures<F> {
    /// Return the number of features (columns) seen during [`Fit::fit`].
    ///
    /// Mirrors scikit-learn's `PolynomialFeatures.n_features_in_`
    /// (`_polynomial.py:323`).
    #[must_use]
    pub fn n_features_in(&self) -> usize {
        self.n_features_in_
    }

    /// Return the total number of polynomial output columns.
    ///
    /// Mirrors scikit-learn's `PolynomialFeatures.n_output_features_`
    /// (`_polynomial.py:362`).
    #[must_use]
    pub fn n_output_features(&self) -> usize {
        self.n_output_features_
    }

    /// Return the `powers_` exponent matrix of shape
    /// `(n_output_features_, n_features_in_)`.
    ///
    /// `powers_[i, j]` is the exponent of input feature `j` in output feature
    /// `i`, mirroring scikit-learn's `PolynomialFeatures.powers_`
    /// (`_polynomial.py:250-264`).
    #[must_use]
    pub fn powers(&self) -> &Array2<usize> {
        &self.powers_
    }

    /// Return the configured degree.
    #[must_use]
    pub fn degree(&self) -> usize {
        self.degree
    }

    /// Return whether only interaction terms are generated.
    #[must_use]
    pub fn interaction_only(&self) -> bool {
        self.interaction_only
    }

    /// Return whether a bias column is included.
    #[must_use]
    pub fn include_bias(&self) -> bool {
        self.include_bias
    }
}

impl<F: Float + Send + Sync + 'static> Fit<Array2<F>, ()> for PolynomialFeatures<F> {
    type Fitted = FittedPolynomialFeatures<F>;
    type Error = FerroError;

    /// Validate the input and record `n_features_in_` / `n_output_features_` /
    /// `powers_`, returning a [`FittedPolynomialFeatures`].
    ///
    /// `PolynomialFeatures.fit` (sklearn `_polynomial.py:306-400`, "Compute
    /// number of output features") learns NO numeric statistics. This runs the
    /// SAME REQ-8 `check_array` validation as [`Transform::transform`] (via the
    /// shared [`validate_poly_input`] helper; sklearn's `_validate_data` default
    /// `force_all_finite=True` REJECTS NaN/±inf), enumerates the polynomial
    /// combinations by REUSING [`PolynomialFeatures::feature_combinations`] (no
    /// reimplementation), and records `n_features_in_ = x.ncols()`,
    /// `n_output_features_ = combinations.len()` (== sklearn's
    /// `_num_combinations`, `:362`), and the `powers_` exponent matrix
    /// (`build_powers` = sklearn's `np.bincount` over the combinations, `:262`).
    /// `powers_` is built from the SAME combinations the transform value math
    /// uses, so its row order equals the output column order.
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::InsufficientSamples`] for zero rows and
    /// [`FerroError::InvalidParameter`] for zero features or any non-finite
    /// value (NaN, +inf, -inf) — matching `check_array`
    /// (`sklearn/utils/validation.py:1084`, `:1093`, `:1063`) as routed through
    /// `PolynomialFeatures.fit` → `_validate_data` (`_polynomial.py:323`).
    fn fit(&self, x: &Array2<F>, _y: &()) -> Result<FittedPolynomialFeatures<F>, FerroError> {
        validate_poly_input(x, "PolynomialFeatures::fit")?;
        let n_features_in_ = x.ncols();
        // REUSE the SHIPPED combination enumeration (REQ-1) — do NOT reimplement.
        let combos = self.feature_combinations(n_features_in_);
        let n_output_features_ = combos.len();
        let powers_ = build_powers(&combos, n_features_in_);
        Ok(FittedPolynomialFeatures {
            degree: self.degree,
            interaction_only: self.interaction_only,
            include_bias: self.include_bias,
            n_features_in_,
            n_output_features_,
            powers_,
            _marker: std::marker::PhantomData,
        })
    }
}

impl<F: Float + Send + Sync + 'static> Transform<Array2<F>> for FittedPolynomialFeatures<F> {
    type Output = Array2<F>;
    type Error = FerroError;

    /// Generate polynomial features, delegating to the SAME combination
    /// enumeration + value math as the stateless [`PolynomialFeatures`] path.
    ///
    /// First applies the REQ-8 `check_array` guards (min-samples / min-features
    /// / finite) and generates the polynomial matrix, THEN validates that `x`
    /// has the same number of columns recorded during [`Fit::fit`]. This ORDER
    /// matches sklearn's `_validate_data(reset=False)`, which runs `check_array`
    /// BEFORE the `n_features_in_` consistency check (`_polynomial.py:433-435`,
    /// #2207): a NaN / ±inf / zero-sample / zero-feature input raises its
    /// `check_array` error EVEN when the column count is also wrong. Only after
    /// that does the feature-count comparison fire. Because the combinations are
    /// re-enumerated from the SAME params (degree / interaction_only /
    /// include_bias) and the value math is the shared [`generate_poly_features`],
    /// the output is byte-identical to `PolynomialFeatures::transform`.
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if the column count differs from
    /// `n_features_in_` (sklearn `ValueError: X has N features, but
    /// PolynomialFeatures is expecting M features as input.`,
    /// `_polynomial.py:402-435`). Returns [`FerroError::InsufficientSamples`]
    /// for zero rows and [`FerroError::InvalidParameter`] for zero features or
    /// any non-finite value (REQ-8, via [`validate_poly_input`]).
    fn transform(&self, x: &Array2<F>) -> Result<Array2<F>, FerroError> {
        // sklearn `_validate_data(reset=False)` runs `check_array` (finite /
        // min-samples / min-features) BEFORE the `n_features_in_` check (#2207).
        // So validate + generate FIRST; a NaN / ±inf / zero-sample / zero-feature
        // input must raise its check_array error EVEN when the column count is
        // also wrong. Only then does the feature-count comparison fire.
        validate_poly_input(x, "FittedPolynomialFeatures::transform")?;
        // REUSE the SHARED incremental value math — build from the configured
        // params so the column order matches the stateless path (and `powers_`)
        // exactly.
        let out = generate_poly_features(x, self.degree, self.interaction_only, self.include_bias);
        if x.ncols() != self.n_features_in_ {
            return Err(FerroError::ShapeMismatch {
                expected: vec![x.nrows(), self.n_features_in_],
                actual: vec![x.nrows(), x.ncols()],
                context: "FittedPolynomialFeatures::transform".into(),
            });
        }
        Ok(out)
    }
}

// ---------------------------------------------------------------------------
// Pipeline integration (generic)
// ---------------------------------------------------------------------------

impl<F: Float + Send + Sync + 'static> PipelineTransformer<F> for PolynomialFeatures<F> {
    /// Fit the polynomial features transformer using the pipeline interface.
    ///
    /// Because `PolynomialFeatures` is stateless, this simply boxes `self`
    /// as a [`FittedPipelineTransformer`].
    ///
    /// # Errors
    ///
    /// This implementation never returns an error.
    fn fit_pipeline(
        &self,
        _x: &Array2<F>,
        _y: &Array1<F>,
    ) -> Result<Box<dyn FittedPipelineTransformer<F>>, FerroError> {
        Ok(Box::new(self.clone()))
    }
}

impl<F: Float + Send + Sync + 'static> FittedPipelineTransformer<F> for PolynomialFeatures<F> {
    /// Transform data using the pipeline interface.
    ///
    /// # Errors
    ///
    /// Propagates errors from [`Transform::transform`].
    fn transform_pipeline(&self, x: &Array2<F>) -> Result<Array2<F>, FerroError> {
        self.transform(x)
    }
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

#[cfg(test)]
mod tests {
    use super::*;
    use approx::assert_abs_diff_eq;
    use ndarray::array;

    #[test]
    fn test_degree2_two_features_with_bias() {
        // degree=2, interaction_only=false, include_bias=true
        // Expected: [1, a, b, a², a·b, b²]
        let poly = PolynomialFeatures::<f64>::new(2, false, true).unwrap();
        let x = array![[2.0, 3.0]];
        let out = poly.transform(&x).unwrap();
        assert_eq!(out.shape()[0], 1);
        assert_eq!(out.shape()[1], 6); // 1 + 2 + 3 combinations
        assert_abs_diff_eq!(out[[0, 0]], 1.0, epsilon = 1e-10); // bias
        assert_abs_diff_eq!(out[[0, 1]], 2.0, epsilon = 1e-10); // a
        assert_abs_diff_eq!(out[[0, 2]], 3.0, epsilon = 1e-10); // b
        assert_abs_diff_eq!(out[[0, 3]], 4.0, epsilon = 1e-10); // a²
        assert_abs_diff_eq!(out[[0, 4]], 6.0, epsilon = 1e-10); // a·b
        assert_abs_diff_eq!(out[[0, 5]], 9.0, epsilon = 1e-10); // b²
    }

    #[test]
    fn test_degree2_interaction_only() {
        // degree=2, interaction_only=true, include_bias=true
        // Expected: [1, a, b, a·b]
        let poly = PolynomialFeatures::<f64>::new(2, true, true).unwrap();
        let x = array![[2.0, 3.0]];
        let out = poly.transform(&x).unwrap();
        assert_eq!(out.shape()[1], 4);
        assert_abs_diff_eq!(out[[0, 0]], 1.0, epsilon = 1e-10); // bias
        assert_abs_diff_eq!(out[[0, 1]], 2.0, epsilon = 1e-10); // a
        assert_abs_diff_eq!(out[[0, 2]], 3.0, epsilon = 1e-10); // b
        assert_abs_diff_eq!(out[[0, 3]], 6.0, epsilon = 1e-10); // a·b
    }

    #[test]
    fn test_no_bias() {
        // degree=2, interaction_only=false, include_bias=false
        // Expected: [a, b, a², a·b, b²]
        let poly = PolynomialFeatures::<f64>::new(2, false, false).unwrap();
        let x = array![[2.0, 3.0]];
        let out = poly.transform(&x).unwrap();
        assert_eq!(out.shape()[1], 5);
        assert_abs_diff_eq!(out[[0, 0]], 2.0, epsilon = 1e-10); // a
    }

    #[test]
    fn test_degree1_only_linear() {
        let poly = PolynomialFeatures::<f64>::new(1, false, true).unwrap();
        let x = array![[2.0, 3.0]];
        let out = poly.transform(&x).unwrap();
        // [1, a, b]
        assert_eq!(out.shape()[1], 3);
        assert_abs_diff_eq!(out[[0, 0]], 1.0, epsilon = 1e-10);
        assert_abs_diff_eq!(out[[0, 1]], 2.0, epsilon = 1e-10);
        assert_abs_diff_eq!(out[[0, 2]], 3.0, epsilon = 1e-10);
    }

    #[test]
    fn test_multiple_rows() {
        let poly = PolynomialFeatures::<f64>::new(2, false, false).unwrap();
        let x = array![[1.0, 2.0], [3.0, 4.0]];
        let out = poly.transform(&x).unwrap();
        assert_eq!(out.shape(), &[2, 5]);
        // Row 0: a=1, b=2 → [1, 2, 1, 2, 4]
        assert_abs_diff_eq!(out[[0, 0]], 1.0, epsilon = 1e-10);
        assert_abs_diff_eq!(out[[0, 1]], 2.0, epsilon = 1e-10);
        assert_abs_diff_eq!(out[[0, 2]], 1.0, epsilon = 1e-10);
        assert_abs_diff_eq!(out[[0, 3]], 2.0, epsilon = 1e-10);
        assert_abs_diff_eq!(out[[0, 4]], 4.0, epsilon = 1e-10);
    }

    #[test]
    fn test_single_feature_degree2() {
        // [a] → [1, a, a²]
        let poly = PolynomialFeatures::<f64>::new(2, false, true).unwrap();
        let x = array![[3.0]];
        let out = poly.transform(&x).unwrap();
        assert_eq!(out.shape()[1], 3);
        assert_abs_diff_eq!(out[[0, 0]], 1.0, epsilon = 1e-10);
        assert_abs_diff_eq!(out[[0, 1]], 3.0, epsilon = 1e-10);
        assert_abs_diff_eq!(out[[0, 2]], 9.0, epsilon = 1e-10);
    }

    #[test]
    fn test_degree0_include_bias_false_is_err() {
        // sklearn `_polynomial.py:326-330`: degree==0 with include_bias=False is
        // the ONLY degree==0 case that raises ("Setting degree to zero and
        // include_bias to False would result in an empty output array.").
        assert!(PolynomialFeatures::<f64>::new(0, false, false).is_err());
        assert!(PolynomialFeatures::<f64>::new(0, true, false).is_err());
    }

    #[test]
    fn test_degree0_include_bias_true_bias_only() -> Result<(), FerroError> {
        // sklearn ORACLE (degree=0, include_bias=True): a single all-ones bias
        // column, powers_ == [[0, 0]] (one all-zeros row), n_output_features_ == 1.
        //   python3 -c "import numpy as np; from sklearn.preprocessing import \
        //   PolynomialFeatures as P; m=P(degree=0, include_bias=True).fit( \
        //   np.array([[2.,3.]])); print(m.transform(np.array([[2.,3.]])).tolist(), \
        //   m.powers_.tolist(), m.n_output_features_)"
        //   -> [[1.0]] [[0, 0]] 1
        let poly = PolynomialFeatures::<f64>::new(0, false, true)?;
        let x = array![[2.0, 3.0]];

        // Stateless transform → single all-ones bias column.
        let out = poly.transform(&x)?;
        assert_eq!(out.shape(), &[1, 1]);
        assert_abs_diff_eq!(out[[0, 0]], 1.0, epsilon = 1e-10);

        // Fitted attributes: powers_ == [[0, 0]], n_output_features_ == 1.
        let fitted = poly.fit(&x, &())?;
        assert_eq!(fitted.n_output_features(), 1);
        assert_eq!(fitted.powers().shape(), &[1, 2]);
        assert_eq!(fitted.powers()[[0, 0]], 0);
        assert_eq!(fitted.powers()[[0, 1]], 0);

        // Fitted transform is byte-identical to the stateless path.
        let out_fitted = fitted.transform(&x)?;
        assert_eq!(out_fitted, out);
        Ok(())
    }

    #[test]
    fn test_default_config() {
        let poly = PolynomialFeatures::<f64>::default();
        assert_eq!(poly.degree(), 2);
        assert!(!poly.interaction_only());
        assert!(poly.include_bias());
    }

    #[test]
    fn test_pipeline_integration() {
        use ferrolearn_core::pipeline::PipelineTransformer;
        let poly = PolynomialFeatures::<f64>::new(2, false, true).unwrap();
        let x = array![[1.0, 2.0], [3.0, 4.0]];
        let y = Array1::zeros(2);
        let fitted = poly.fit_pipeline(&x, &y).unwrap();
        let result = fitted.transform_pipeline(&x).unwrap();
        assert_eq!(result.shape(), &[2, 6]);
    }

    #[test]
    fn test_degree3_single_feature() {
        // [a] with degree=3, no bias → [a, a², a³]
        let poly = PolynomialFeatures::<f64>::new(3, false, false).unwrap();
        let x = array![[2.0]];
        let out = poly.transform(&x).unwrap();
        assert_eq!(out.shape()[1], 3);
        assert_abs_diff_eq!(out[[0, 0]], 2.0, epsilon = 1e-10);
        assert_abs_diff_eq!(out[[0, 1]], 4.0, epsilon = 1e-10);
        assert_abs_diff_eq!(out[[0, 2]], 8.0, epsilon = 1e-10);
    }
}