ferrolearn_linear/

linear_regression.rs

//! Ordinary Least Squares linear regression.
//!
//! This module provides [`LinearRegression`], which fits a linear model by
//! solving the least squares problem via a single SVD (the LAPACK-`gelsd`
//! minimum-norm path, through `ferray::linalg::lstsq`):
//!
//! ```text
//! minimize ||X @ w - y||^2
//! ```
//!
//! ## REQ status (per `.design/linear/linear_regression.md`, mirrors `sklearn/linear_model/_base.py` @ 1.5.2)
//!
//! Mirrors `sklearn.linear_model.LinearRegression` (`_base.py:465`). Full-rank,
//! rank-deficient, and underdetermined OLS all match the live sklearn oracle to
//! 1e-8: the solve routes through `crate::linalg::solve_lstsq` →
//! `ferray::linalg::lstsq` (single-SVD, LAPACK-`gelsd`-equivalent min-norm),
//! mirroring sklearn's `linalg.lstsq(X, y)` (`_base.py:687`).
//!
//! | REQ | Status | Evidence |
//! |---|---|---|
//! | REQ-1 (full-rank OLS coef_/intercept_) | SHIPPED | `Fit for LinearRegression` (centering + `linalg::solve_lstsq` via `ferray::linalg::lstsq`); full-rank coef/intercept match oracle to 1e-8. Consumer: `RsLinearRegression` in `ferrolearn-python/src/regressors.rs`. Mirrors `_base.py:582`, intercept `_base.py:308`. |
//! | REQ-2 (predict = X·coef + intercept) | SHIPPED | `Predict for FittedLinearRegression`. Mirrors `_base.py:282`. |
//! | REQ-3 (fit_intercept incl. false) | SHIPPED | `with_fit_intercept`; `fit_intercept=false` forces intercept 0. Mirrors `_base.py:571`. |
//! | REQ-4 (HasCoefficients introspection) | SHIPPED | `HasCoefficients for FittedLinearRegression`. Mirrors fitted attrs `_base.py:499/511`. |
//! | REQ-5 (min-norm for rank-deficient / underdetermined X) | SHIPPED | `Fit for LinearRegression` calls `crate::linalg::solve_lstsq` → `ferray::linalg::lstsq` (`ferray-linalg/src/solve.rs:208`), the single-SVD gelsd-equivalent min-norm solver mirroring `_base.py:687`. Closes #376 (rank-deficient min-norm) + #377 (underdetermined accepted). Tests now passing (`#[ignore]` removed): `divergence_rank_deficient_no_intercept_min_norm`, `divergence_rank_deficient_with_intercept_min_norm`, `divergence_underdetermined_accepted_min_norm` in `tests/divergence_linreg_minnorm.rs`. |
//! | REQ-6 (positive=True / NNLS) | SHIPPED | `LinearRegression<F>` adds `pub positive: bool` (default `false`, `_base.py:574`) + `with_positive(bool)` builder. `fit_with_sample_weight`'s coefficient solve routes through `solve_coef`, which calls `crate::linalg::nnls` (Lawson-Hanson active-set NNLS solving the passive-set unconstrained LS via `solve_lstsq` on the passive columns) instead of `solve_lstsq` when `self.positive`, on the SAME centered-and-`√w`-rescaled design — mirroring sklearn's `self.coef_ = optimize.nnls(X, y)[0]` (`_base.py:647`) after `_preprocess_data`/`_rescale_data`. Intercept recovered identically (`y_off − x_off·coef` when fit_intercept, else 0; `_set_intercept`, `_base.py:692`). `rank_`/`singular_` are still taken from the `solve_lstsq` SVD of the design (sklearn leaves them unset on the positive path; ferrolearn reports the design's SVD as a documented analog). `positive=false` (default) is byte-identical to the unconstrained OLS path. Oracle tests: `linreg_positive_matches_sklearn` (coef `[2.03571429, 0.0]`, intercept `-1.46428571`, all ≥ 0, differs from unconstrained `[2.25, -0.75]`), `linreg_positive_no_intercept_matches_sklearn` (raw `nnls(X,y)` `[1.34210526, 0.0]`), `linreg_positive_false_unchanged` (byte-identical guard); `nnls_matches_scipy`/`nnls_equals_ols_when_unconstrained_nonneg` in `linalg.rs`. Closes #371. |
//! | REQ-7 (multi-output 2-D Y → 2-D coef_) | SHIPPED | Additive `Fit<Array2<F>, Array2<F>>` arm (does NOT touch the 1-D `Fit`/`FittedLinearRegression`/`Predict`) producing `FittedMultiOutputLinearRegression<F>` — `coefficients` shape `(n_targets, n_features)` (sklearn `coef_` orientation, `coef_.T` of the lstsq solution, `_base.py:688`), `intercepts` `(n_targets,)`, `rank_`/`singular_`. Solves all targets in one SVD via `linalg::solve_lstsq_multi` → `ferray::linalg::lstsq` with a 2-D `b` (mirrors `linalg.lstsq(X, Y)`, `_base.py:687`); shared X-centering + per-target y-offset, `intercepts = y_off − coefficients · x_off` (`_set_intercept`, `_base.py:322`); `fit_intercept=false` → raw solve, intercepts 0. `Predict<Array2<F>, Output=Array2<F>>` returns `X · coef_.T + intercepts` shape `(n_samples, n_targets)` (`_base.py:290`). Oracle tests `linreg_multioutput_coef_intercept_match_sklearn` (coef `[[2.06666667,-0.06666667],[0.86666667,0.23333333]]`, intercept `[-0.06666667,0.13333333]`), `linreg_multioutput_predict_shape_and_values` (`predict(X[:2]) = [[2.0,1.0],[4.0,2.1]]`), `linreg_multioutput_no_intercept` (coef `[[2.0195121951,-0.0097560976],[0.9609756098,0.1195121951]]`, intercepts 0), `linreg_single_output_unchanged` (1-D path byte-identical). Closes #372. |
//! | REQ-8 (sample_weight in fit) | SHIPPED | `LinearRegression::fit_with_sample_weight` solves WEIGHTED least squares `min Σᵢ wᵢ(yᵢ−xᵢ·w)²`: weighted offsets `x_off[j]=Σwᵢx[i,j]/Σwᵢ`, `y_off=Σwᵢyᵢ/Σwᵢ` (mirrors `_average(...,weights=sample_weight)`, `_base.py:193`/`:198`), centering, then `√wᵢ` row-rescaling (`_rescale_data`, `_base.py:641`), `linalg::solve_lstsq` on the rescaled design, `intercept = y_off − x_off·coef` (`_set_intercept`, `_base.py:320`); `fit_intercept=false` skips centering, intercept 0. `Fit::fit` delegates `fit_with_sample_weight(x, y, None)` (None path byte-identical to the historic OLS body). Oracle tests `linreg_fit_sample_weight_with_intercept_matches_sklearn` (coef 2.0935828877, intercept −0.2326203209), `linreg_fit_sample_weight_no_intercept_matches_sklearn` (coef 2.0350877193, intercept 0), `linreg_fit_none_sample_weight_equals_unweighted`. Mirrors `fit(..., sample_weight=None)` (`_base.py:582`). Closes #373. |
//! | REQ-9 (rank_/singular_/copy_X/n_jobs) | SHIPPED | `FittedLinearRegression` stores `rank_`/`singular_` (captured from `linalg::solve_lstsq` on the matrix actually solved — centered `X` when `fit_intercept`, raw `X` otherwise, matching sklearn `_base.py:687`), exposed via `rank()`/`singular_values()`; `LinearRegression` adds `copy_x` (default `true`) + `n_jobs` (default `None`) fields with `with_copy_x`/`with_n_jobs` builders, mirroring `_parameter_constraints` (`_base.py:561`) and the ctor (`_base.py:572-573`). `copy_x` is ABI-only (fit never mutates `x`); `n_jobs` stored-but-ignored (single-threaded). Oracle tests `linreg_rank_singular_match_sklearn_with_intercept` (rank 2, singular `[1.61803399, 0.61803399]` on centered X), `linreg_singular_no_intercept_matches_raw_x` (singular `[5.25371017, 0.63129192]` on raw X), `linreg_copy_x_default_and_builder`. Closes #374. |
//! | REQ-10 (ferray substrate) | NOT-STARTED | blocker #375 — OLS solve now on `ferray::linalg::lstsq`, but `LinearRegression`'s coef storage is still `ndarray` (coef return type tied to #359); fully on-substrate when the boundary `ndarray` types migrate. |
//! | REQ-11 (non-finite input rejected) | SHIPPED | `fit_with_sample_weight` (the shared entry `Fit::fit` delegates to) rejects any NaN/+/-inf in X or y BEFORE centering/solve with `FerroError::InvalidParameter`, mirroring sklearn's `_validate_data(force_all_finite=True)` (`_base.py:609`, default `force_all_finite=True` → `check_array` raises `ValueError("Input X contains NaN.")` / `"... contains infinity ..."`). `.iter().any(|v| !v.is_finite())` catches both NaN and Inf; the finite path is byte-identical (the guard never fires on finite input). Verified vs the live sklearn 1.5.2 oracle (R-CHAR-3): `LinearRegression().fit` raises `ValueError` for NaN/+inf/-inf in X and NaN/inf in y (`tests/divergence_linear_nonfinite.rs::linreg_*`). Non-test consumer: the existing `Fit::fit` / `RsLinearRegression` consumers. (#2256) |
//!
//! Two states only per goal.md R-DEFER-2. The OLS min-norm contract (#376/#377)
//! is fixed in `linalg.rs` via the ferray substrate.
//!
//! # Examples
//!
//! ```
//! use ferrolearn_linear::LinearRegression;
//! use ferrolearn_core::{Fit, Predict};
//! use ndarray::{array, Array1, Array2};
//!
//! let model = LinearRegression::<f64>::new();
//! let x = Array2::from_shape_vec((4, 1), vec![1.0, 2.0, 3.0, 4.0]).unwrap();
//! let y = array![2.0, 4.0, 6.0, 8.0];
//!
//! let fitted = model.fit(&x, &y).unwrap();
//! let preds = fitted.predict(&x).unwrap();
//! ```

use ferrolearn_core::error::FerroError;
use ferrolearn_core::introspection::HasCoefficients;
use ferrolearn_core::pipeline::{FittedPipelineEstimator, PipelineEstimator};
use ferrolearn_core::traits::{Fit, Predict};
use ndarray::{Array1, Array2, Axis, ScalarOperand};
use num_traits::{Float, FromPrimitive};

use crate::linalg;

/// Ordinary least squares linear regression.
///
/// Solves the least-squares problem via a single SVD (minimum-norm,
/// LAPACK-`gelsd`-equivalent, through `ferray::linalg::lstsq`). The
/// `fit_intercept` option controls whether a bias (intercept) term is
/// included.
///
/// # Type Parameters
///
/// - `F`: The floating-point type (`f32` or `f64`).
#[derive(Debug, Clone)]
pub struct LinearRegression<F> {
    /// Whether to fit an intercept (bias) term.
    pub fit_intercept: bool,
    /// Whether `X` may be overwritten during fit (sklearn `copy_X`,
    /// `_base.py:480`). ferrolearn's `fit` never mutates `x` (it reads via
    /// `.iter()`/`.mean_axis()`), so the observable non-mutation contract
    /// holds for either value; the field is exposed for ABI parity. Default
    /// `true`, matching sklearn (`_base.py:572`).
    pub copy_x: bool,
    /// Number of jobs for the computation (sklearn `n_jobs`, `_base.py:483`).
    /// ferrolearn's dense OLS solve is single-threaded, so this is stored but
    /// ignored — parallelism is a no-op here and behaviour matches sklearn's
    /// `n_jobs=None` single-job default. Default `None` (`_base.py:573`).
    pub n_jobs: Option<usize>,
    /// When `true`, constrains the fitted coefficients to be non-negative via
    /// non-negative least squares (sklearn `positive`, `_base.py:574`). sklearn
    /// solves the (centered, `√w`-rescaled) coefficient system with
    /// `scipy.optimize.nnls` instead of `linalg.lstsq` when `positive=True`
    /// (`_base.py:645-647`). Default `false`, matching sklearn's
    /// `positive=False` (`_base.py:574`); when `false`, the fit is
    /// byte-identical to the unconstrained OLS path.
    pub positive: bool,
    _marker: std::marker::PhantomData<F>,
}

impl<
    F: Float
        + Send
        + Sync
        + ScalarOperand
        + num_traits::FromPrimitive
        + ferray::linalg::LinalgFloat
        + 'static,
> LinearRegression<F>
{
    /// Fit the linear regression model with optional per-sample weights.
    ///
    /// Mirrors scikit-learn's `LinearRegression.fit(X, y, sample_weight=None)`
    /// (`sklearn/linear_model/_base.py:582`). When `sample_weight` is `Some(w)`,
    /// this solves the WEIGHTED least-squares problem `min Σᵢ wᵢ (yᵢ − xᵢ·w)²`:
    ///
    /// - `fit_intercept=true`: offsets are the WEIGHTED means
    ///   `x_off[j] = Σᵢ wᵢ·x[i,j] / Σwᵢ`, `y_off = Σᵢ wᵢ·yᵢ / Σwᵢ`
    ///   (sklearn `_preprocess_data` → `_average(..., weights=sample_weight)`,
    ///   `_base.py:193`/`:198`). `X` and `y` are centered by those offsets, each
    ///   row is then rescaled by `√wᵢ` (sklearn `_rescale_data`, `_base.py:641`),
    ///   `linalg.lstsq` solves for `coef`, and
    ///   `intercept = y_off − x_off · coef` (`_set_intercept`, `_base.py:320`).
    /// - `fit_intercept=false`: no centering; each row is rescaled by `√wᵢ`, the
    ///   solve runs on the rescaled `X`, and `intercept = 0`.
    ///
    /// `sample_weight=None` is BYTE-IDENTICAL to [`Fit::fit`] (the unweighted
    /// centering + `solve_lstsq` path), which delegates here.
    ///
    /// `rank_`/`singular_` are captured from `solve_lstsq` on the matrix actually
    /// solved (centered-and-rescaled `X` when `fit_intercept`, rescaled `X`
    /// otherwise), matching sklearn's `linalg.lstsq` operands (`_base.py:687`).
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if the number of samples in `x` and
    /// `y` (or, when provided, `sample_weight`) differ.
    /// Returns [`FerroError::InsufficientSamples`] if there are no samples.
    /// Returns [`FerroError::NumericalInstability`] if the system is singular or
    /// the weighted-offset denominator (`Σwᵢ`) cannot be formed.
    /// Solve the coefficient system on the (already centered / `√w`-rescaled)
    /// design `a` and target `b`, returning `(coef, rank_, singular_)`.
    ///
    /// `rank_`/`singular_` always come from the unconstrained `linalg.lstsq`
    /// SVD of the design `a` (matching sklearn's `linalg.lstsq(X, y)` operands,
    /// `_base.py:687`). When `self.positive`, the COEFFICIENTS are overridden by
    /// the non-negative least-squares solution (`scipy.optimize.nnls`,
    /// `_base.py:647`) on the same design; otherwise the lstsq coefficients are
    /// returned unchanged, keeping the `positive=false` path byte-identical to
    /// the unconstrained OLS solve.
    fn solve_coef(
        &self,
        a: &Array2<F>,
        b: &Array1<F>,
    ) -> Result<(Array1<F>, usize, Array1<F>), FerroError> {
        let (coef, rank, singular) = linalg::solve_lstsq(a, b)?;
        if self.positive {
            let coef_pos = linalg::nnls(a, b)?;
            Ok((coef_pos, rank, singular))
        } else {
            Ok((coef, rank, singular))
        }
    }

    pub fn fit_with_sample_weight(
        &self,
        x: &Array2<F>,
        y: &Array1<F>,
        sample_weight: Option<&Array1<F>>,
    ) -> Result<FittedLinearRegression<F>, FerroError> {
        let (n_samples, n_features) = x.dim();

        // Validate input shapes.
        if n_samples != y.len() {
            return Err(FerroError::ShapeMismatch {
                expected: vec![n_samples],
                actual: vec![y.len()],
                context: "y length must match number of samples in X".into(),
            });
        }

        if n_samples == 0 {
            return Err(FerroError::InsufficientSamples {
                required: 1,
                actual: 0,
                context: "LinearRegression requires at least one sample".into(),
            });
        }

        if let Some(w) = sample_weight
            && w.len() != n_samples
        {
            return Err(FerroError::ShapeMismatch {
                expected: vec![n_samples],
                actual: vec![w.len()],
                context: "sample_weight length must match number of samples in X".into(),
            });
        }

        // sklearn `LinearRegression.fit` -> `self._validate_data(X, y, ...)`
        // (`_base.py:609`); the call keeps the default `force_all_finite=True`,
        // so `check_array` rejects any NaN or +/-inf in X OR y with a
        // `ValueError` BEFORE the solve. `.iter().any(|v| !v.is_finite())`
        // rejects both NaN and Inf (bounds-safe, no panic, R-CODE-2), matching
        // the crate idiom (`multi_task_lasso.rs`). (#2256)
        if x.iter().any(|v| !v.is_finite()) {
            return Err(FerroError::InvalidParameter {
                name: "X".into(),
                reason: "Input X contains NaN or infinity.".into(),
            });
        }
        if y.iter().any(|v| !v.is_finite()) {
            return Err(FerroError::InvalidParameter {
                name: "y".into(),
                reason: "Input y contains NaN or infinity.".into(),
            });
        }

        // sklearn validates `sample_weight` via `_check_sample_weight` ->
        // `check_array(..., input_name="sample_weight")`
        // (`sklearn/utils/validation.py:2043-2050`), keeping the default
        // `force_all_finite=True`, so any NaN or +/-inf weight raises a
        // `ValueError` BEFORE the weighted centering / √w rescaling. Mirror it
        // with the same NaN+Inf-rejecting idiom as X/y above (#2258).
        if let Some(w) = sample_weight
            && w.iter().any(|v| !v.is_finite())
        {
            return Err(FerroError::InvalidParameter {
                name: "sample_weight".into(),
                reason: "Input sample_weight contains NaN or infinity.".into(),
            });
        }

        match sample_weight {
            None => {
                // Unweighted path — identical to the original `Fit::fit` body.
                if self.fit_intercept {
                    // Centering trick: center X and y, solve the (uncentered)
                    // OLS problem on the centered design, then recover the
                    // intercept as y_mean - x_mean . w. sklearn centers
                    // identically before its `linalg.lstsq` call (`_base.py`
                    // `_preprocess_data` + `:687`).
                    let n = <F as num_traits::NumCast>::from(n_samples).ok_or_else(|| {
                        FerroError::NumericalInstability {
                            message: "could not represent n_samples as the float type".into(),
                        }
                    })?;
                    let x_mean =
                        x.mean_axis(Axis(0))
                            .ok_or_else(|| FerroError::InsufficientSamples {
                                required: 1,
                                actual: 0,
                                context: "cannot compute feature means of an empty design".into(),
                            })?;
                    let y_mean = y.sum() / n;

                    let x_centered = x - &x_mean;
                    let y_centered = y - y_mean;

                    let (w, rank, singular) = self.solve_coef(&x_centered, &y_centered)?;

                    let intercept = y_mean - x_mean.dot(&w);

                    Ok(FittedLinearRegression {
                        coefficients: w,
                        intercept,
                        rank_: rank,
                        singular_: singular,
                    })
                } else {
                    let (w, rank, singular) = self.solve_coef(x, y)?;

                    Ok(FittedLinearRegression {
                        coefficients: w,
                        intercept: <F as num_traits::Zero>::zero(),
                        rank_: rank,
                        singular_: singular,
                    })
                }
            }
            Some(w) => {
                // Per-row √w factor (sklearn `_rescale_data`, `_base.py:641`).
                let w_sqrt = w.mapv(<F as Float>::sqrt);

                if self.fit_intercept {
                    // WEIGHTED centering: offsets are the weighted means
                    // x_off[j] = Σ wᵢ x[i,j] / Σ wᵢ, y_off = Σ wᵢ yᵢ / Σ wᵢ
                    // (sklearn `_average(..., weights=sample_weight)`,
                    // `_base.py:193`/`:198`).
                    let w_sum = w.sum();
                    if w_sum <= <F as num_traits::Zero>::zero() {
                        return Err(FerroError::NumericalInstability {
                            message: "sum of sample_weight must be positive to center".into(),
                        });
                    }

                    let mut x_off = Array1::<F>::zeros(n_features);
                    for (i, row) in x.outer_iter().enumerate() {
                        let wi = w[i];
                        x_off = &x_off + &row.mapv(|v| v * wi);
                    }
                    x_off.mapv_inplace(|v| v / w_sum);

                    let y_off = y
                        .iter()
                        .zip(w.iter())
                        .fold(<F as num_traits::Zero>::zero(), |acc, (&yi, &wi)| {
                            acc + wi * yi
                        })
                        / w_sum;

                    // Center, then row-rescale by √w.
                    let x_centered = x - &x_off;
                    let y_centered = y - y_off;
                    let x_scaled = &x_centered * &w_sqrt.view().insert_axis(Axis(1));
                    let y_scaled = &y_centered * &w_sqrt;

                    let (coef, rank, singular) = self.solve_coef(&x_scaled, &y_scaled)?;

                    let intercept = y_off - x_off.dot(&coef);

                    Ok(FittedLinearRegression {
                        coefficients: coef,
                        intercept,
                        rank_: rank,
                        singular_: singular,
                    })
                } else {
                    // No centering; just √w row-rescaling, intercept 0.
                    let x_scaled = x * &w_sqrt.view().insert_axis(Axis(1));
                    let y_scaled = y * &w_sqrt;

                    let (coef, rank, singular) = self.solve_coef(&x_scaled, &y_scaled)?;

                    Ok(FittedLinearRegression {
                        coefficients: coef,
                        intercept: <F as num_traits::Zero>::zero(),
                        rank_: rank,
                        singular_: singular,
                    })
                }
            }
        }
    }
}

impl<F: Float> LinearRegression<F> {
    /// Create a new `LinearRegression` with default settings.
    ///
    /// Defaults: `fit_intercept = true`, `copy_x = true`, `n_jobs = None`,
    /// `positive = false` (mirroring sklearn's ctor defaults,
    /// `_base.py:571-574`).
    #[must_use]
    pub fn new() -> Self {
        Self {
            fit_intercept: true,
            copy_x: true,
            n_jobs: None,
            positive: false,
            _marker: std::marker::PhantomData,
        }
    }

    /// Set whether to fit an intercept term.
    #[must_use]
    pub fn with_fit_intercept(mut self, fit_intercept: bool) -> Self {
        self.fit_intercept = fit_intercept;
        self
    }

    /// Set the `copy_X` flag (sklearn `copy_X`, `_base.py:480`).
    ///
    /// ferrolearn's fit never mutates `x`, so this is exposed for ABI parity
    /// with sklearn and does not change the result.
    #[must_use]
    pub fn with_copy_x(mut self, copy_x: bool) -> Self {
        self.copy_x = copy_x;
        self
    }

    /// Set the `n_jobs` parameter (sklearn `n_jobs`, `_base.py:483`).
    ///
    /// The dense OLS solve is single-threaded; this is stored but ignored.
    #[must_use]
    pub fn with_n_jobs(mut self, n_jobs: Option<usize>) -> Self {
        self.n_jobs = n_jobs;
        self
    }

    /// Set the `positive` flag (sklearn `positive`, `_base.py:574`).
    ///
    /// When `true`, the fitted coefficients are constrained to be non-negative,
    /// solved via non-negative least squares (`scipy.optimize.nnls`,
    /// `_base.py:647`) instead of unconstrained OLS. Default `false`.
    #[must_use]
    pub fn with_positive(mut self, positive: bool) -> Self {
        self.positive = positive;
        self
    }
}

impl<F: Float> Default for LinearRegression<F> {
    fn default() -> Self {
        Self::new()
    }
}

/// Fitted ordinary least squares linear regression model.
///
/// Stores the learned coefficients and intercept. Implements [`Predict`]
/// to generate predictions and [`HasCoefficients`] for introspection.
#[derive(Debug, Clone)]
pub struct FittedLinearRegression<F> {
    /// Learned coefficient vector (one per feature).
    coefficients: Array1<F>,
    /// Learned intercept (bias) term.
    intercept: F,
    /// Effective rank of the design matrix actually solved (sklearn `rank_`,
    /// `_base.py:505`/`:687`) — the centered `X` when `fit_intercept`, the
    /// raw `X` otherwise.
    rank_: usize,
    /// Singular values of the design matrix actually solved (sklearn
    /// `singular_`, `_base.py:508`/`:687`).
    singular_: Array1<F>,
}

impl<
    F: Float
        + Send
        + Sync
        + ScalarOperand
        + num_traits::FromPrimitive
        + ferray::linalg::LinalgFloat
        + 'static,
> Fit<Array2<F>, Array1<F>> for LinearRegression<F>
{
    type Fitted = FittedLinearRegression<F>;
    type Error = FerroError;

    /// Fit the linear regression model.
    ///
    /// Solves the OLS least-squares problem via the SVD-based
    /// minimum-norm solver [`crate::linalg::solve_lstsq`] (routed through
    /// [`ferray::linalg::lstsq`], LAPACK-`gelsd`-equivalent), matching
    /// scikit-learn's dense path `linalg.lstsq(X, y)`
    /// (`sklearn/linear_model/_base.py:687`). When `fit_intercept` is true,
    /// `X` and `y` are centered first and the intercept is recovered as
    /// `y_mean - x_mean . w`.
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if the number of samples in `x`
    /// and `y` differ.
    /// Returns [`FerroError::InsufficientSamples`] if there are fewer samples
    /// than features.
    /// Returns [`FerroError::NumericalInstability`] if the system is singular.
    fn fit(&self, x: &Array2<F>, y: &Array1<F>) -> Result<FittedLinearRegression<F>, FerroError> {
        // Unweighted OLS is the `sample_weight=None` arm of the weighted fit;
        // delegating keeps the None path byte-identical to the historic body
        // (centering + `solve_lstsq`), mirroring sklearn's single `fit` entry
        // (`_base.py:582`, `sample_weight=None` default).
        self.fit_with_sample_weight(x, y, None)
    }
}

/// Fitted multi-output ordinary least squares linear regression model.
///
/// The 2-D-target companion to [`FittedLinearRegression`]: produced by
/// `Fit<Array2<F>, Array2<F>>` when fitting a 2-D `Y` of shape
/// `(n_samples, n_targets)`. Mirrors scikit-learn's multi-output
/// `LinearRegression` (`MultiOutputMixin`, `_base.py:465`), whose `coef_` is a
/// 2-D array of shape `(n_targets, n_features)` and `intercept_` an array of
/// shape `(n_targets,)` (`_base.py:499`/`:511`). Stored in sklearn's `coef_`
/// orientation (target rows), so `coefficients()` maps directly onto
/// `sklearn.coef_`.
#[derive(Debug, Clone)]
pub struct FittedMultiOutputLinearRegression<F> {
    /// Learned coefficient matrix in sklearn `coef_` orientation: shape
    /// `(n_targets, n_features)`, row `t` the coefficients for target `t`
    /// (`_base.py:499`).
    coefficients: Array2<F>,
    /// Learned per-target intercepts, shape `(n_targets,)` (sklearn
    /// `intercept_`, `_base.py:511`).
    intercepts: Array1<F>,
    /// Effective rank of the design matrix actually solved (sklearn `rank_`,
    /// `_base.py:505`/`:687`) — the centered `X` when `fit_intercept`, the raw
    /// `X` otherwise.
    rank_: usize,
    /// Singular values of the design matrix actually solved (sklearn
    /// `singular_`, `_base.py:508`/`:687`).
    singular_: Array1<F>,
}

impl<F: Float> FittedMultiOutputLinearRegression<F> {
    /// Learned coefficient matrix, shape `(n_targets, n_features)` (sklearn
    /// `coef_`, `_base.py:499`).
    #[must_use]
    pub fn coefficients(&self) -> &Array2<F> {
        &self.coefficients
    }

    /// Learned per-target intercepts, shape `(n_targets,)` (sklearn
    /// `intercept_`, `_base.py:511`).
    #[must_use]
    pub fn intercepts(&self) -> &Array1<F> {
        &self.intercepts
    }

    /// Effective rank of the design matrix (sklearn `rank_`, `_base.py:505`).
    #[must_use]
    pub fn rank(&self) -> usize {
        self.rank_
    }

    /// Singular values of the design matrix (sklearn `singular_`,
    /// `_base.py:508`).
    #[must_use]
    pub fn singular_values(&self) -> &Array1<F> {
        &self.singular_
    }
}

impl<
    F: Float
        + Send
        + Sync
        + ScalarOperand
        + num_traits::FromPrimitive
        + ferray::linalg::LinalgFloat
        + 'static,
> Fit<Array2<F>, Array2<F>> for LinearRegression<F>
{
    type Fitted = FittedMultiOutputLinearRegression<F>;
    type Error = FerroError;

    /// Fit the multi-output linear regression model on a 2-D target `Y`.
    ///
    /// Mirrors scikit-learn's multi-output dense path: `linalg.lstsq(X, Y)`
    /// with `Y` of shape `(n_samples, n_targets)` solves all targets in one
    /// SVD, yielding `coef_` of shape `(n_targets, n_features)` and a per-target
    /// `intercept_` of shape `(n_targets,)` (`sklearn/linear_model/_base.py:687`,
    /// `coef_.T`; intercept `_set_intercept`, `_base.py:308`/`:322`). When
    /// `fit_intercept` is true, `X` and each column of `Y` are centered by their
    /// column means and the intercept is recovered as
    /// `y_off − coefficients · x_off` per target; when false, the solve runs on
    /// raw `X`/`Y` and the intercepts are all `0`.
    ///
    /// The 1-D `Fit<Array2<F>, Array1<F>>` impl is unaffected — this is an
    /// additive 2-D arm.
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if the number of samples in `x`
    /// and `y` differ.
    /// Returns [`FerroError::InsufficientSamples`] if there are no samples.
    /// Returns [`FerroError::NumericalInstability`] if the system is singular.
    fn fit(
        &self,
        x: &Array2<F>,
        y: &Array2<F>,
    ) -> Result<FittedMultiOutputLinearRegression<F>, FerroError> {
        let n_samples = x.nrows();
        let n_targets = y.ncols();

        if n_samples != y.nrows() {
            return Err(FerroError::ShapeMismatch {
                expected: vec![n_samples],
                actual: vec![y.nrows()],
                context: "Y rows must match number of samples in X".into(),
            });
        }

        if n_samples == 0 {
            return Err(FerroError::InsufficientSamples {
                required: 1,
                actual: 0,
                context: "LinearRegression requires at least one sample".into(),
            });
        }

        // sklearn `LinearRegression.fit` -> `self._validate_data(X, y, ...,
        // multi_output=True, ...)` (`_base.py:609`) keeps the default
        // `force_all_finite=True`, so `check_array` rejects any NaN or +/-inf in
        // X OR the 2-D Y with a `ValueError` BEFORE the solve — regardless of
        // output dimensionality. This separate multi-output arm does NOT route
        // through `fit_with_sample_weight`, so it needs the SAME finite-check.
        // `.iter().any(|v| !v.is_finite())` (Array2's element iterator) rejects
        // both NaN and Inf (bounds-safe, no panic, R-CODE-2). (#2257)
        if x.iter().any(|v| !v.is_finite()) {
            return Err(FerroError::InvalidParameter {
                name: "X".into(),
                reason: "Input X contains NaN or infinity.".into(),
            });
        }
        if y.iter().any(|v| !v.is_finite()) {
            return Err(FerroError::InvalidParameter {
                name: "y".into(),
                reason: "Input y contains NaN or infinity.".into(),
            });
        }

        if self.fit_intercept {
            // Same column-centering as the 1-D fit, generalized to Y's columns
            // (sklearn `_preprocess_data` centers X and every column of Y by
            // their per-column means, `_base.py:193`/`:198`).
            let x_off = x
                .mean_axis(Axis(0))
                .ok_or_else(|| FerroError::InsufficientSamples {
                    required: 1,
                    actual: 0,
                    context: "cannot compute feature means of an empty design".into(),
                })?;
            let y_off = y
                .mean_axis(Axis(0))
                .ok_or_else(|| FerroError::InsufficientSamples {
                    required: 1,
                    actual: 0,
                    context: "cannot compute target means of an empty Y".into(),
                })?;

            let x_centered = x - &x_off;
            let y_centered = y - &y_off;

            // coef_ft is (n_features, n_targets); store in sklearn `coef_`
            // orientation (n_targets, n_features) via transpose.
            let (coef_ft, rank, singular) = linalg::solve_lstsq_multi(&x_centered, &y_centered)?;
            let coefficients = coef_ft.t().to_owned();

            // intercept_[t] = y_off[t] − coefficients[t] · x_off
            // (sklearn `_set_intercept`: `y_offset − X_offset @ coef_.T`,
            // `_base.py:322`).
            let intercepts = &y_off - &coefficients.dot(&x_off);

            Ok(FittedMultiOutputLinearRegression {
                coefficients,
                intercepts,
                rank_: rank,
                singular_: singular,
            })
        } else {
            let (coef_ft, rank, singular) = linalg::solve_lstsq_multi(x, y)?;
            let coefficients = coef_ft.t().to_owned();
            let intercepts = Array1::<F>::zeros(n_targets);

            Ok(FittedMultiOutputLinearRegression {
                coefficients,
                intercepts,
                rank_: rank,
                singular_: singular,
            })
        }
    }
}

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

    /// Predict 2-D target values for the given feature matrix.
    ///
    /// Computes `X @ coefficients.T + intercepts` (broadcasting the per-target
    /// intercepts over rows), shape `(n_samples, n_targets)`, mirroring sklearn's
    /// 2-D `_decision_function` arm `X @ coef_.T + self.intercept_`
    /// (`_base.py:290`).
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if the number of features does not
    /// match the fitted model.
    fn predict(&self, x: &Array2<F>) -> Result<Array2<F>, FerroError> {
        let n_features = x.ncols();
        if n_features != self.coefficients.ncols() {
            return Err(FerroError::ShapeMismatch {
                expected: vec![self.coefficients.ncols()],
                actual: vec![n_features],
                context: "number of features must match fitted model".into(),
            });
        }

        // X (n_samples, n_features) @ coef_.T (n_features, n_targets) -> (n_samples, n_targets)
        let preds = x.dot(&self.coefficients.t());
        Ok(preds + &self.intercepts)
    }
}

impl<F: Float> FittedLinearRegression<F> {
    /// Effective rank of the design matrix (sklearn `rank_`, `_base.py:505`).
    ///
    /// The rank of the matrix actually solved by `linalg.lstsq` — the
    /// centered `X` when `fit_intercept` is true, the raw `X` otherwise
    /// (`_base.py:687`).
    #[must_use]
    pub fn rank(&self) -> usize {
        self.rank_
    }

    /// Singular values of the design matrix (sklearn `singular_`,
    /// `_base.py:508`).
    ///
    /// The singular values of the matrix actually solved by `linalg.lstsq`
    /// — the centered `X` when `fit_intercept` is true, the raw `X`
    /// otherwise (`_base.py:687`).
    #[must_use]
    pub fn singular_values(&self) -> &Array1<F> {
        &self.singular_
    }
}

impl<F: Float + Send + Sync + ScalarOperand + 'static> Predict<Array2<F>>
    for FittedLinearRegression<F>
{
    type Output = Array1<F>;
    type Error = FerroError;

    /// Predict target values for the given feature matrix.
    ///
    /// Computes `X @ coefficients + intercept`.
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if the number of features
    /// does not match the fitted model.
    fn predict(&self, x: &Array2<F>) -> Result<Array1<F>, FerroError> {
        let n_features = x.ncols();
        if n_features != self.coefficients.len() {
            return Err(FerroError::ShapeMismatch {
                expected: vec![self.coefficients.len()],
                actual: vec![n_features],
                context: "number of features must match fitted model".into(),
            });
        }

        let preds = x.dot(&self.coefficients) + self.intercept;
        Ok(preds)
    }
}

impl<F: Float + Send + Sync + ScalarOperand + 'static> HasCoefficients<F>
    for FittedLinearRegression<F>
{
    fn coefficients(&self) -> &Array1<F> {
        &self.coefficients
    }

    fn intercept(&self) -> F {
        self.intercept
    }
}

// Pipeline integration.
impl<F> PipelineEstimator<F> for LinearRegression<F>
where
    F: Float + FromPrimitive + ScalarOperand + ferray::linalg::LinalgFloat + Send + Sync + 'static,
{
    fn fit_pipeline(
        &self,
        x: &Array2<F>,
        y: &Array1<F>,
    ) -> Result<Box<dyn FittedPipelineEstimator<F>>, FerroError> {
        let fitted = self.fit(x, y)?;
        Ok(Box::new(fitted))
    }
}

impl<F> FittedPipelineEstimator<F> for FittedLinearRegression<F>
where
    F: Float + ScalarOperand + Send + Sync + 'static,
{
    fn predict_pipeline(&self, x: &Array2<F>) -> Result<Array1<F>, FerroError> {
        self.predict(x)
    }
}

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

    #[test]
    fn test_simple_linear_regression() {
        // y = 2*x + 1
        let x = Array2::from_shape_vec((5, 1), vec![1.0, 2.0, 3.0, 4.0, 5.0]).unwrap();
        let y = array![3.0, 5.0, 7.0, 9.0, 11.0];

        let model = LinearRegression::<f64>::new();
        let fitted = model.fit(&x, &y).unwrap();

        assert_relative_eq!(fitted.coefficients()[0], 2.0, epsilon = 1e-10);
        assert_relative_eq!(fitted.intercept(), 1.0, epsilon = 1e-10);

        let preds = fitted.predict(&x).unwrap();
        for (p, &actual) in preds.iter().zip(y.iter()) {
            assert_relative_eq!(*p, actual, epsilon = 1e-10);
        }
    }

    #[test]
    fn test_multiple_linear_regression() {
        // y = 1*x1 + 2*x2 + 3
        let x =
            Array2::from_shape_vec((4, 2), vec![1.0, 1.0, 2.0, 1.0, 3.0, 2.0, 4.0, 2.0]).unwrap();
        let y = array![6.0, 7.0, 10.0, 11.0];

        let model = LinearRegression::<f64>::new();
        let fitted = model.fit(&x, &y).unwrap();

        assert_relative_eq!(fitted.coefficients()[0], 1.0, epsilon = 1e-10);
        assert_relative_eq!(fitted.coefficients()[1], 2.0, epsilon = 1e-10);
        assert_relative_eq!(fitted.intercept(), 3.0, epsilon = 1e-10);
    }

    #[test]
    fn test_no_intercept() {
        // y = 2*x (through origin)
        let x = Array2::from_shape_vec((4, 1), vec![1.0, 2.0, 3.0, 4.0]).unwrap();
        let y = array![2.0, 4.0, 6.0, 8.0];

        let model = LinearRegression::<f64>::new().with_fit_intercept(false);
        let fitted = model.fit(&x, &y).unwrap();

        assert_relative_eq!(fitted.coefficients()[0], 2.0, epsilon = 1e-10);
        assert_relative_eq!(fitted.intercept(), 0.0, epsilon = 1e-10);
    }

    #[test]
    fn test_shape_mismatch_fit() {
        let x = Array2::from_shape_vec((3, 1), vec![1.0, 2.0, 3.0]).unwrap();
        let y = array![1.0, 2.0]; // Wrong length

        let model = LinearRegression::<f64>::new();
        let result = model.fit(&x, &y);
        assert!(result.is_err());
    }

    #[test]
    fn test_shape_mismatch_predict() {
        let x = Array2::from_shape_vec((3, 2), vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]).unwrap();
        let y = array![1.0, 2.0, 3.0];

        let model = LinearRegression::<f64>::new();
        let fitted = model.fit(&x, &y).unwrap();

        // Wrong number of features
        let x_bad = Array2::from_shape_vec((2, 3), vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]).unwrap();
        let result = fitted.predict(&x_bad);
        assert!(result.is_err());
    }

    #[test]
    fn test_has_coefficients() {
        let x = Array2::from_shape_vec((3, 1), vec![1.0, 2.0, 3.0]).unwrap();
        let y = array![2.0, 4.0, 6.0];

        let model = LinearRegression::<f64>::new();
        let fitted = model.fit(&x, &y).unwrap();

        assert_eq!(fitted.coefficients().len(), 1);
    }

    #[test]
    fn test_pipeline_integration() {
        let x = Array2::from_shape_vec((4, 1), vec![1.0, 2.0, 3.0, 4.0]).unwrap();
        let y = array![3.0, 5.0, 7.0, 9.0];

        let model = LinearRegression::<f64>::new();
        let fitted = model.fit_pipeline(&x, &y).unwrap();
        let preds = fitted.predict_pipeline(&x).unwrap();
        assert_eq!(preds.len(), 4);
    }

    #[test]
    fn linreg_rank_singular_match_sklearn_with_intercept() {
        // Live sklearn 1.5.2 oracle (fit_intercept=True centers X before
        // linalg.lstsq, so singular_ are the singular values of CENTERED X):
        //   cd /tmp && python3 -c "import numpy as np; \
        //     from sklearn.linear_model import LinearRegression; \
        //     X=np.array([[1.,1.],[1.,2.],[2.,2.],[2.,3.]]); \
        //     y=np.array([6.,8.,9.,11.]); m=LinearRegression().fit(X,y); \
        //     print(m.rank_, [round(s,8) for s in m.singular_])"
        //   -> 2 [1.61803399, 0.61803399]
        let x =
            Array2::from_shape_vec((4, 2), vec![1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 2.0, 3.0]).unwrap();
        let y = array![6.0, 8.0, 9.0, 11.0];

        let model = LinearRegression::<f64>::new();
        let fitted = model.fit(&x, &y).unwrap();

        assert_eq!(fitted.rank(), 2);
        let sv = fitted.singular_values();
        assert_eq!(sv.len(), 2);
        assert_relative_eq!(sv[0], 1.618_033_99, epsilon = 1e-6);
        assert_relative_eq!(sv[1], 0.618_033_99, epsilon = 1e-6);
    }

    #[test]
    fn linreg_singular_no_intercept_matches_raw_x() {
        // Live sklearn 1.5.2 oracle (fit_intercept=False → singular_ are the
        // singular values of the RAW X):
        //   cd /tmp && python3 -c "import numpy as np; \
        //     from sklearn.linear_model import LinearRegression; \
        //     X=np.array([[1.,1.],[1.,2.],[2.,2.],[2.,3.]]); \
        //     y=np.array([6.,8.,9.,11.]); \
        //     m=LinearRegression(fit_intercept=False).fit(X,y); \
        //     print(m.rank_, [round(s,8) for s in m.singular_])"
        //   -> 2 [5.25371017, 0.63129192]
        let x =
            Array2::from_shape_vec((4, 2), vec![1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 2.0, 3.0]).unwrap();
        let y = array![6.0, 8.0, 9.0, 11.0];

        let model = LinearRegression::<f64>::new().with_fit_intercept(false);
        let fitted = model.fit(&x, &y).unwrap();

        assert_eq!(fitted.rank(), 2);
        let sv = fitted.singular_values();
        assert_eq!(sv.len(), 2);
        assert_relative_eq!(sv[0], 5.253_710_17, epsilon = 1e-6);
        assert_relative_eq!(sv[1], 0.631_291_92, epsilon = 1e-6);
    }

    #[test]
    fn linreg_copy_x_default_and_builder() {
        // copy_X default is true (sklearn `_base.py:572`); the builder flips
        // it; n_jobs builder stores Some(4); and fit produces identical coef_
        // regardless of copy_x (no behaviour change — fit never mutates X).
        assert!(LinearRegression::<f64>::new().copy_x);
        assert!(!LinearRegression::<f64>::new().with_copy_x(false).copy_x);
        assert_eq!(
            LinearRegression::<f64>::new().with_n_jobs(Some(4)).n_jobs,
            Some(4)
        );

        let x = Array2::from_shape_vec((4, 1), vec![1.0, 2.0, 3.0, 4.0]).unwrap();
        let y = array![3.0, 5.0, 7.0, 9.0];

        let fitted_copy = LinearRegression::<f64>::new()
            .with_copy_x(true)
            .fit(&x, &y)
            .unwrap();
        let fitted_nocopy = LinearRegression::<f64>::new()
            .with_copy_x(false)
            .fit(&x, &y)
            .unwrap();

        assert_relative_eq!(
            fitted_copy.coefficients()[0],
            fitted_nocopy.coefficients()[0],
            epsilon = 1e-12
        );
        assert_relative_eq!(
            fitted_copy.intercept(),
            fitted_nocopy.intercept(),
            epsilon = 1e-12
        );
    }

    #[test]
    fn linreg_fit_sample_weight_with_intercept_matches_sklearn() {
        // Live sklearn 1.5.2 oracle (WEIGHTED OLS, fit_intercept=True):
        //   cd /tmp && python3 -c "import numpy as np; \
        //     from sklearn.linear_model import LinearRegression; \
        //     X=np.array([[1.],[2.],[3.],[4.],[5.]]); \
        //     y=np.array([2.1,3.9,6.2,7.7,10.3]); w=np.array([1.,5.,1.,1.,5.]); \
        //     m=LinearRegression().fit(X,y,sample_weight=w); \
        //     print(round(m.coef_[0],10), round(m.intercept_,10))"
        //   -> 2.0935828877 -0.2326203209
        let x = Array2::from_shape_vec((5, 1), vec![1.0, 2.0, 3.0, 4.0, 5.0]).unwrap();
        let y = array![2.1, 3.9, 6.2, 7.7, 10.3];
        let w = array![1.0, 5.0, 1.0, 1.0, 5.0];

        let model = LinearRegression::<f64>::new();
        let fitted = model.fit_with_sample_weight(&x, &y, Some(&w)).unwrap();

        assert_relative_eq!(fitted.coefficients()[0], 2.093_582_887_7, epsilon = 1e-7);
        assert_relative_eq!(fitted.intercept(), -0.232_620_320_9, epsilon = 1e-7);

        // Non-tautological: the weighted result MUST differ from the unweighted
        // fit (oracle unweighted coef_ 2.02, intercept_ -0.02).
        let unweighted = model.fit(&x, &y).unwrap();
        assert_relative_eq!(unweighted.coefficients()[0], 2.02, epsilon = 1e-7);
        assert!((fitted.coefficients()[0] - unweighted.coefficients()[0]).abs() > 1e-3);
        assert!((fitted.intercept() - unweighted.intercept()).abs() > 1e-3);
    }

    #[test]
    fn linreg_fit_sample_weight_no_intercept_matches_sklearn() {
        // Live sklearn 1.5.2 oracle (WEIGHTED OLS, fit_intercept=False):
        //   cd /tmp && python3 -c "import numpy as np; \
        //     from sklearn.linear_model import LinearRegression; \
        //     X=np.array([[1.],[2.],[3.],[4.],[5.]]); \
        //     y=np.array([2.1,3.9,6.2,7.7,10.3]); w=np.array([1.,5.,1.,1.,5.]); \
        //     m=LinearRegression(fit_intercept=False).fit(X,y,sample_weight=w); \
        //     print(round(m.coef_[0],10))"
        //   -> 2.0350877193
        let x = Array2::from_shape_vec((5, 1), vec![1.0, 2.0, 3.0, 4.0, 5.0]).unwrap();
        let y = array![2.1, 3.9, 6.2, 7.7, 10.3];
        let w = array![1.0, 5.0, 1.0, 1.0, 5.0];

        let model = LinearRegression::<f64>::new().with_fit_intercept(false);
        let fitted = model.fit_with_sample_weight(&x, &y, Some(&w)).unwrap();

        assert_relative_eq!(fitted.coefficients()[0], 2.035_087_719_3, epsilon = 1e-7);
        assert_eq!(fitted.intercept(), 0.0);
    }

    #[test]
    fn linreg_fit_none_sample_weight_equals_unweighted() {
        // Regression guard: the `None` path is BYTE-IDENTICAL to `fit`.
        let x = Array2::from_shape_vec((5, 1), vec![1.0, 2.0, 3.0, 4.0, 5.0]).unwrap();
        let y = array![2.1, 3.9, 6.2, 7.7, 10.3];

        let model = LinearRegression::<f64>::new();
        let via_fit = model.fit(&x, &y).unwrap();
        let via_none = model.fit_with_sample_weight(&x, &y, None).unwrap();

        assert_eq!(
            via_fit.coefficients()[0].to_bits(),
            via_none.coefficients()[0].to_bits()
        );
        assert_eq!(
            via_fit.intercept().to_bits(),
            via_none.intercept().to_bits()
        );

        // Same for fit_intercept=false.
        let model_ni = LinearRegression::<f64>::new().with_fit_intercept(false);
        let via_fit_ni = model_ni.fit(&x, &y).unwrap();
        let via_none_ni = model_ni.fit_with_sample_weight(&x, &y, None).unwrap();
        assert_eq!(
            via_fit_ni.coefficients()[0].to_bits(),
            via_none_ni.coefficients()[0].to_bits()
        );
        assert_eq!(
            via_fit_ni.intercept().to_bits(),
            via_none_ni.intercept().to_bits()
        );
    }

    #[test]
    fn linreg_multioutput_coef_intercept_match_sklearn() {
        // Live sklearn 1.5.2 oracle (multi-output, fit_intercept=True):
        //   cd /tmp && python3 -c "import numpy as np; \
        //     from sklearn.linear_model import LinearRegression; \
        //     X=np.array([[1.,0.],[2.,1.],[3.,1.],[4.,2.],[5.,3.]]); \
        //     Y=np.array([[2.1,1.0],[3.9,2.1],[6.2,2.9],[7.7,4.2],[10.3,5.1]]); \
        //     m=LinearRegression().fit(X,Y); print(m.coef_.shape); \
        //     print([[round(v,8) for v in r] for r in m.coef_]); \
        //     print([round(v,8) for v in m.intercept_])"
        //   -> (2, 2)
        //      [[2.06666667, -0.06666667], [0.86666667, 0.23333333]]
        //      [-0.06666667, 0.13333333]
        let x = Array2::from_shape_vec(
            (5, 2),
            vec![1.0, 0.0, 2.0, 1.0, 3.0, 1.0, 4.0, 2.0, 5.0, 3.0],
        )
        .unwrap();
        let y = Array2::from_shape_vec(
            (5, 2),
            vec![2.1, 1.0, 3.9, 2.1, 6.2, 2.9, 7.7, 4.2, 10.3, 5.1],
        )
        .unwrap();

        let model = LinearRegression::<f64>::new();
        let fitted = Fit::<Array2<f64>, Array2<f64>>::fit(&model, &x, &y).unwrap();

        assert_eq!(fitted.coefficients().dim(), (2, 2));
        let c = fitted.coefficients();
        assert_relative_eq!(c[[0, 0]], 2.066_666_67, epsilon = 1e-7);
        assert_relative_eq!(c[[0, 1]], -0.066_666_67, epsilon = 1e-7);
        assert_relative_eq!(c[[1, 0]], 0.866_666_67, epsilon = 1e-7);
        assert_relative_eq!(c[[1, 1]], 0.233_333_33, epsilon = 1e-7);

        let b = fitted.intercepts();
        assert_eq!(b.len(), 2);
        assert_relative_eq!(b[0], -0.066_666_67, epsilon = 1e-7);
        assert_relative_eq!(b[1], 0.133_333_33, epsilon = 1e-7);
    }

    #[test]
    fn linreg_multioutput_predict_shape_and_values() {
        // Oracle (same model as above): predict(X).shape == (5, 2);
        //   m.predict(X[:2]) -> [[2.0, 1.0], [4.0, 2.1]]
        let x = Array2::from_shape_vec(
            (5, 2),
            vec![1.0, 0.0, 2.0, 1.0, 3.0, 1.0, 4.0, 2.0, 5.0, 3.0],
        )
        .unwrap();
        let y = Array2::from_shape_vec(
            (5, 2),
            vec![2.1, 1.0, 3.9, 2.1, 6.2, 2.9, 7.7, 4.2, 10.3, 5.1],
        )
        .unwrap();

        let model = LinearRegression::<f64>::new();
        let fitted = Fit::<Array2<f64>, Array2<f64>>::fit(&model, &x, &y).unwrap();

        let preds = fitted.predict(&x).unwrap();
        assert_eq!(preds.dim(), (5, 2));

        let x2 = x.slice(ndarray::s![0..2, ..]).to_owned();
        let preds2 = fitted.predict(&x2).unwrap();
        assert_eq!(preds2.dim(), (2, 2));
        assert_relative_eq!(preds2[[0, 0]], 2.0, epsilon = 1e-6);
        assert_relative_eq!(preds2[[0, 1]], 1.0, epsilon = 1e-6);
        assert_relative_eq!(preds2[[1, 0]], 4.0, epsilon = 1e-6);
        assert_relative_eq!(preds2[[1, 1]], 2.1, epsilon = 1e-6);
    }

    #[test]
    fn linreg_multioutput_no_intercept() {
        // Live sklearn 1.5.2 oracle (multi-output, fit_intercept=False):
        //   cd /tmp && python3 -c "import numpy as np; \
        //     from sklearn.linear_model import LinearRegression; \
        //     X=np.array([[1.,0.],[2.,1.],[3.,1.],[4.,2.],[5.,3.]]); \
        //     Y=np.array([[2.1,1.0],[3.9,2.1],[6.2,2.9],[7.7,4.2],[10.3,5.1]]); \
        //     m=LinearRegression(fit_intercept=False).fit(X,Y); \
        //     print([[round(v,10) for v in r] for r in m.coef_]); print(m.intercept_)"
        //   -> [[2.0195121951, -0.0097560976], [0.9609756098, 0.1195121951]]
        //      0.0
        let x = Array2::from_shape_vec(
            (5, 2),
            vec![1.0, 0.0, 2.0, 1.0, 3.0, 1.0, 4.0, 2.0, 5.0, 3.0],
        )
        .unwrap();
        let y = Array2::from_shape_vec(
            (5, 2),
            vec![2.1, 1.0, 3.9, 2.1, 6.2, 2.9, 7.7, 4.2, 10.3, 5.1],
        )
        .unwrap();

        let model = LinearRegression::<f64>::new().with_fit_intercept(false);
        let fitted = Fit::<Array2<f64>, Array2<f64>>::fit(&model, &x, &y).unwrap();

        let c = fitted.coefficients();
        assert_eq!(c.dim(), (2, 2));
        assert_relative_eq!(c[[0, 0]], 2.019_512_195_1, epsilon = 1e-7);
        assert_relative_eq!(c[[0, 1]], -0.009_756_097_6, epsilon = 1e-7);
        assert_relative_eq!(c[[1, 0]], 0.960_975_609_8, epsilon = 1e-7);
        assert_relative_eq!(c[[1, 1]], 0.119_512_195_1, epsilon = 1e-7);

        let b = fitted.intercepts();
        assert_eq!(b.len(), 2);
        assert_eq!(b[0], 0.0);
        assert_eq!(b[1], 0.0);
    }

    #[test]
    fn linreg_single_output_unchanged() {
        // Regression guard: the additive 2-D arm must not disturb the 1-D path.
        // y = 2*x + 1 (same fixture as `test_simple_linear_regression`).
        let x = Array2::from_shape_vec((5, 1), vec![1.0, 2.0, 3.0, 4.0, 5.0]).unwrap();
        let y1 = array![3.0, 5.0, 7.0, 9.0, 11.0];

        let model = LinearRegression::<f64>::new();
        let fitted = model.fit(&x, &y1).unwrap();

        assert_relative_eq!(fitted.coefficients()[0], 2.0, epsilon = 1e-10);
        assert_relative_eq!(fitted.intercept(), 1.0, epsilon = 1e-10);
    }

    #[test]
    fn linreg_positive_matches_sklearn() {
        // Live sklearn 1.5.2 oracle (positive=True, fit_intercept=True →
        // centers, runs scipy.optimize.nnls on the centered design, recovers
        // intercept = y_off − X_off·coef):
        //   cd /tmp && python3 -c "import numpy as np; \
        //     from sklearn.linear_model import LinearRegression; \
        //     X=np.array([[1.,1.],[1.,2.],[2.,1.],[3.,2.],[2.,3.]]); \
        //     y=np.array([1.,0.5,3.,5.,1.5]); \
        //     m=LinearRegression(positive=True).fit(X,y); \
        //     print([round(c,8) for c in m.coef_], round(m.intercept_,8))"
        //   -> [2.03571429, 0.0] -1.46428571
        // The 2nd coef CLAMPS to 0; the unconstrained fit is [2.25, -0.75].
        let x = array![[1.0, 1.0], [1.0, 2.0], [2.0, 1.0], [3.0, 2.0], [2.0, 3.0]];
        let y = array![1.0, 0.5, 3.0, 5.0, 1.5];

        let model = LinearRegression::<f64>::new().with_positive(true);
        let res = model.fit(&x, &y);
        let unc_res = LinearRegression::<f64>::new().fit(&x, &y);
        assert!(res.is_ok());
        assert!(unc_res.is_ok());
        if let (Ok(fitted), Ok(unconstrained)) = (res, unc_res) {
            assert_relative_eq!(fitted.coefficients()[0], 2.035_714_29, epsilon = 1e-6);
            assert_relative_eq!(fitted.coefficients()[1], 0.0, epsilon = 1e-6);
            assert_relative_eq!(fitted.intercept(), -1.464_285_71, epsilon = 1e-6);

            // Non-negativity contract.
            assert!(fitted.coefficients().iter().all(|&c| c >= 0.0));

            // Non-tautological: the constrained result MUST differ from the
            // unconstrained OLS fit (oracle coef_ [2.25, -0.75]).
            assert_relative_eq!(unconstrained.coefficients()[0], 2.25, epsilon = 1e-6);
            assert_relative_eq!(unconstrained.coefficients()[1], -0.75, epsilon = 1e-6);
            assert!((fitted.coefficients()[1] - unconstrained.coefficients()[1]).abs() > 0.5);
        }
    }

    #[test]
    fn linreg_positive_false_unchanged() {
        // Regression guard: with_positive(false) (the default) is
        // BYTE-IDENTICAL to the historic unconstrained fit.
        let x = array![[1.0, 1.0], [1.0, 2.0], [2.0, 1.0], [3.0, 2.0], [2.0, 3.0]];
        let y = array![1.0, 0.5, 3.0, 5.0, 1.5];

        let d_res = LinearRegression::<f64>::new().fit(&x, &y);
        let e_res = LinearRegression::<f64>::new()
            .with_positive(false)
            .fit(&x, &y);
        assert!(d_res.is_ok());
        assert!(e_res.is_ok());
        if let (Ok(default), Ok(explicit)) = (d_res, e_res) {
            assert_eq!(
                default.coefficients()[0].to_bits(),
                explicit.coefficients()[0].to_bits()
            );
            assert_eq!(
                default.coefficients()[1].to_bits(),
                explicit.coefficients()[1].to_bits()
            );
            assert_eq!(
                default.intercept().to_bits(),
                explicit.intercept().to_bits()
            );

            // And matches the unconstrained oracle (coef_ [2.25, -0.75]).
            assert_relative_eq!(default.coefficients()[0], 2.25, epsilon = 1e-6);
            assert_relative_eq!(default.coefficients()[1], -0.75, epsilon = 1e-6);
        }
    }

    #[test]
    fn linreg_positive_no_intercept_matches_sklearn() {
        // Live sklearn 1.5.2 oracle (positive=True, fit_intercept=False →
        // raw nnls(X, y), intercept 0):
        //   cd /tmp && python3 -c "import numpy as np; \
        //     from sklearn.linear_model import LinearRegression; \
        //     X=np.array([[1.,1.],[1.,2.],[2.,1.],[3.,2.],[2.,3.]]); \
        //     y=np.array([1.,0.5,3.,5.,1.5]); \
        //     m=LinearRegression(positive=True,fit_intercept=False).fit(X,y); \
        //     print([round(c,8) for c in m.coef_], round(m.intercept_,8))"
        //   -> [1.34210526, 0.0] 0.0  (== raw nnls(X, y))
        let x = array![[1.0, 1.0], [1.0, 2.0], [2.0, 1.0], [3.0, 2.0], [2.0, 3.0]];
        let y = array![1.0, 0.5, 3.0, 5.0, 1.5];

        let res = LinearRegression::<f64>::new()
            .with_positive(true)
            .with_fit_intercept(false)
            .fit(&x, &y);
        assert!(res.is_ok());
        if let Ok(fitted) = res {
            assert_relative_eq!(fitted.coefficients()[0], 1.342_105_26, epsilon = 1e-6);
            assert_relative_eq!(fitted.coefficients()[1], 0.0, epsilon = 1e-6);
            assert_eq!(fitted.intercept(), 0.0);
            assert!(fitted.coefficients().iter().all(|&c| c >= 0.0));
        }
    }

    #[test]
    fn test_f32_support() {
        let x = Array2::from_shape_vec((4, 1), vec![1.0f32, 2.0, 3.0, 4.0]).unwrap();
        let y = Array1::from_vec(vec![2.0f32, 4.0, 6.0, 8.0]);

        let model = LinearRegression::<f32>::new();
        let fitted = model.fit(&x, &y).unwrap();
        let preds = fitted.predict(&x).unwrap();
        assert_eq!(preds.len(), 4);
    }
}