ferrolearn_linear/

lda.rs

//! Linear Discriminant Analysis (LDA).
//!
//! LDA is both a supervised dimensionality-reduction technique and a linear
//! classifier. This module mirrors scikit-learn's **default** `solver="svd"`
//! path (`sklearn/discriminant_analysis.py:487-559`, commit 156ef14): rather
//! than forming a covariance and solving the classical `Sw⁻¹·Sb` Fisher
//! eigenproblem, it whitens the within-class data with two SVDs, derives the
//! whitened projection `scalings_` and the weighted overall mean `xbar_`, then
//! forms the **affine** classifier `coef_`/`intercept_` (whose `intercept_`
//! embeds `log(priors_)`).
//!
//! The [`Solver::Lsqr`] least-squares path (`_solve_lstsq`,
//! `discriminant_analysis.py:365-419`) is also available (`LDA::with_solver`):
//! it forms the prior-weighted within-class covariance `covariance_` and solves
//! `coef_ = lstsq(covariance_, means_.T).T`, supporting covariance
//! [`Shrinkage`] (`None`/`Auto` Ledoit-Wolf/`Fixed`); it does NOT do
//! dimensionality reduction (no `transform`). The [`Solver::Eigen`]
//! generalized-eigenvalue path (`_solve_eigen`,
//! `discriminant_analysis.py:421-485`) is also available (`LDA::with_solver`):
//! it forms the within-class scatter `Sw = covariance_` and total scatter
//! `St = cov(X)`, then solves the generalized symmetric-definite eigenproblem
//! `eigh(Sb, Sw)` (with `Sb = St - Sw`) — reduced to a STANDARD symmetric
//! eigenproblem via the Cholesky factor of `Sw` (`Sw = L·Lᵀ`,
//! `M = L⁻¹·Sb·L⁻ᵀ`, `eigh(M)`, `evecs = L⁻ᵀ·W`) since ferray exposes only the
//! standard solver. It supports `shrinkage` and dimensionality reduction
//! (`transform` is the un-centered `X @ scalings_`, since eigen has no `xbar_`).
//!
//! # Algorithm (`_solve_svd`, `discriminant_analysis.py:487-559`)
//!
//! With `n = n_samples`, `c = n_classes`:
//! 1. `priors_`: empirical `n_k / n` when the constructor `priors` is `None`
//!    (sklearn's default), else the provided `priors` used VERBATIM
//!    (`discriminant_analysis.py:601-605`).
//! 2. `means_` = per-class mean; `xbar_ = priors_ @ means_`.
//! 3. `Xc` = each sample minus its class mean (stacked); `std = std(Xc, axis=0)`
//!    (population, `ddof=0`), zeros replaced by `1`.
//! 4. `Xw = sqrt(1/(n-c)) · (Xc / std)`; thin SVD `Xw = U·diag(S)·Vt`;
//!    `rank = Σ(S > tol)`; `scalings = (Vt[:rank]/std).T / S[:rank]`.
//! 5. Between-class scaled centers `Xb = (sqrt(n·priors_·1/(c-1)) ⊙
//!    (means_-xbar_).T).T @ scalings`; thin SVD `Xb = U2·diag(S2)·Vt2`;
//!    `explained_variance_ratio_ = (S2²/ΣS2²)[:max_components]`;
//!    `rank2 = Σ(S2 > tol·S2[0])`; `scalings_ = scalings @ Vt2.T[:, :rank2]`.
//! 6. `coef = (means_-xbar_) @ scalings_`;
//!    `intercept_ = -½·Σ(coef²) + log(priors_)`;
//!    `coef_ = coef @ scalings_.T`; `intercept_ -= xbar_ @ coef_.T`.
//!
//! Inference (the `LinearClassifierMixin`, `discriminant_analysis.py:739`):
//! - `transform(X) = ((X - xbar_) @ scalings_)[:, :max_components]`
//!   (`discriminant_analysis.py:684-689`).
//! - `decision_function(X) = X @ coef_.T + intercept_`
//!   (`discriminant_analysis.py:739`).
//! - `predict(X)` = `classes_[argmax(decision_function)]`.
//! - `predict_proba(X)` = `softmax(decision_function)`
//!   (`discriminant_analysis.py:706-711`).
//!
//! The number of discriminant directions is at most `min(n_classes - 1,
//! n_features)`.
//!
//! ## REQ status (per `.design/linear/lda.md`, mirrors `sklearn/discriminant_analysis.py` @ 1.5.2)
//!
//! | REQ | Status | Evidence |
//! |---|---|---|
//! | REQ-1 (svd fit + decision_function parity) | SHIPPED | `_solve_svd` in `fn fit` (`fn svd_s_vt` → `ferray::linalg::svd`) builds `coef_`/`intercept_`/`xbar_`/`scalings_` (`discriminant_analysis.py:556-559`); `fn decision_function` = `X @ coef_.T + intercept_` (`:739`). Consumer: `Predict for FittedLDA` + crate-root `pub use`. Test `lda_decision_function_parity` <1e-6 vs live oracle. #588. |
//! | REQ-2 (predict argmax) | SHIPPED | `Predict::predict` = `classes_[argmax(decision_function)]`; the affine decision carries `log(priors_)` via `intercept_`. Test `lda_imbalanced_priors_predict` (prior shifts the boundary, label-for-label vs live oracle). #589. |
//! | REQ-3 (predict_proba prior-aware) | SHIPPED | `FittedLDA::predict_proba` = `softmax(decision_function)` (`discriminant_analysis.py:711`); rows sum to 1. Test `lda_imbalanced_priors_predict` proba block <1e-6 vs live oracle. #590 (partial: multiclass softmax; binary `expit` collapse pends #600). |
//! | REQ-5 (transform parity) | SHIPPED | `fn transform` = `((X - xbar_) @ scalings_)[:, :max_components]` (`discriminant_analysis.py:684-689`). Test `lda_transform_parity` <1e-6 (per-column sign) vs live oracle. #592. |
//! | REQ-6 (n_components bound) | SHIPPED | `fn fit` computes `max_components = min(n_classes-1, n_features)`, defaults `None` to it, errors `Some(0)`/`Some(k>max)` (`discriminant_analysis.py:614-625`). Tests `test_lda_default_n_components`, `test_lda_error_zero_n_components`, `test_lda_error_n_components_too_large`. |
//! | REQ-7 (priors: None=empirical + provided) | SHIPPED | `fn fit` resolves `priors_`: empirical `n_k/n` when `priors` is `None` (`discriminant_analysis.py:601-603`), else the provided `LDA::with_priors` array, now VALIDATED like sklearn LDA (`:607-612`, unlike QDA): R-DEV-4 length check (`p.len() != n_classes` → `ShapeMismatch`, sklearn would mis-index `:540,557`); negative entries → `InvalidParameter` (`:607-608`, `raise ValueError("priors must be non-negative")`); renormalized `p / p.sum()` with an `eprintln!` warning (the crate's warning channel, cf. qda.rs) when `|Σ-1| > 1e-5` (`:610-612`). The resolved priors flow into `xbar_ = priors_ @ means_` (`:517`), the between-class scaling `sqrt(n·priors_·fac)` (`:540`), and `intercept_ += log(priors_)` (`:557`). `FittedLDA::priors` exposes `priors_`. Consumer: the resolved `priors` is read by `fn fit` (xbar_/scaling/intercept_); `Predict for FittedLDA` consumes the prior-shifted decision. Tests: `lda_imbalanced_priors_predict` (empirical `[0.9091,0.0909]` flips the label), `lda_provided_priors` (`with_priors([0.9,0.1])` `predict_proba` <1e-6 vs live oracle; empirical default differs), `lda_priors_negative_rejected` (`[-0.1,1.1]` → `Err`), `lda_priors_renormalized` (`[0.5,0.6]` → `priors_=[0.4545…,0.5454…]`, `predict_proba` <1e-6 vs the live oracle which renormalizes internally). #593, #603. |
//! | REQ-8 (coef_/intercept_/xbar_) | SHIPPED | `FittedLDA::{coef, intercept, xbar}` accessors expose the `_solve_svd` arrays (`discriminant_analysis.py:556-559,517`). Consumer: `fn decision_function` reads `coef_`/`intercept_`; `fn transform` reads `xbar_`. Verified via `lda_decision_function_parity` (decision = `X@coef_.T+intercept_`) + `lda_transform_parity` (uses `xbar_`). |
//! | REQ-13 (explained_variance_ratio_) | SHIPPED | `fn fit` sets `explained_variance_ratio_ = (S2²/ΣS2²)[:max_components]` from the SECOND (between-class) SVD (`discriminant_analysis.py:550-552`). Test `test_lda_explained_variance_ratio_oracle` <1e-9 vs live `L().explained_variance_ratio_`. #599. |
//! | REQ-4 (predict_log_proba smallest_normal floor) | SHIPPED | `FittedLDA::predict_log_proba` mirrors sklearn exactly (`discriminant_analysis.py:713-737`): `predict_proba` then `prediction[prediction == 0.0] += smallest_normal` (`F::min_positive_value()` = numpy `finfo.smallest_normal`, `:729-736`) before `log`, so nonzero probas keep their true `ln` and exact zeros become `log(MIN_POSITIVE)` (not `-inf`). Consumer: shares `FittedLDA::predict_proba` (the `Predict` path). Test `lda_predict_log_proba` (overlapping 3-class, all-finite log-probas) <1e-6 vs live `LinearDiscriminantAnalysis().predict_log_proba`. #591. |
//! | REQ-9 (lsqr solver) | SHIPPED | `Solver::Lsqr` (`LDA::with_solver`) dispatches `fn fit` to `fn solve_lstsq` (sklearn `_solve_lstsq`, `discriminant_analysis.py:365-419`): `covariance_ = Σ_k priors_[k] · cov(X_k)` (`_class_cov` `:128-172`, ALWAYS populated for lsqr, `:413`); `coef_ = lstsq(covariance_, means_.T)[0].T` (`:416`) via `fn lstsq_multi` → `ferray::linalg::lstsq` (multi-RHS, `ferray-linalg/src/solve.rs:208`, R-SUBSTRATE-4 bridge); `intercept_ = -½·diag(means_ @ coef_.T) + log(priors_)` (`:417-418`). No `scalings_`/`xbar_`/`explained_variance_ratio_` / `transform` (sklearn raises for lsqr `transform`, `:676-679`; `max_components=0` ⇒ empty projection). Binary collapse `coef_[1]-coef_[0]` deferred to #600 (coef_ stays `(n_classes, n_features)`, matching the svd path). Consumer: `fn fit` reads `self.solver` and dispatches; `Predict`/`predict_proba` for `FittedLDA` consume the lsqr `coef_`/`intercept_`. Test `lda_lsqr_solver` (collapsed `coef_[1]-coef_[0]` = `[14.7368…, 14.7368…]`, predict/predict_proba) <1e-6 vs live `LinearDiscriminantAnalysis(solver='lsqr').fit(X,y)`. #595. |
//! | REQ-10 (eigen solver) | SHIPPED | `Solver::Eigen` (`LDA::with_solver`) dispatches `fn fit` to `fn solve_eigen` (sklearn `_solve_eigen`, `discriminant_analysis.py:421-485`): `Sw = covariance_ = Σ_k priors_[k]·cov(X_k)` (`_class_cov` `:467-471`); `St = _cov(WHOLE X, shrinkage)` (total scatter, `:472`); `Sb = St - Sw` (`:473`); the GENERALIZED symmetric-definite eigenproblem `eigh(Sb, Sw)` (`:475`) reduced to STANDARD form via the Cholesky factor of `Sw` (ferray has `eigh`/`cholesky` but no generalized solver): `Sw = L·Lᵀ` (`fn cholesky_lower` → `ferray::linalg::cholesky`, `ferray-linalg/src/decomp/cholesky.rs:22`), `M = L⁻¹·Sb·L⁻ᵀ` (`fn matrix_inverse` → `ferray::linalg::inv`, `ferray-linalg/src/solve.rs:367`) SYMMETRIZED `M = (M+Mᵀ)/2`, `(evals, W) = eigh(M)` (`fn eigh_sym` → `ferray::linalg::eigh`, ascending, `ferray-linalg/src/decomp/eigen.rs:105`), generalized `evecs = L⁻ᵀ·W` sorted DESCENDING by eigenvalue (`:479`); `explained_variance_ratio_ = sort(evals/Σevals)[::-1][:max_components]` (`:476-478`); `scalings_ = evecs` (ALL columns, `:481`); `coef_ = (means_@evecs)@evecs.T` (`:482`, SIGN/ORDER-INVARIANT so it matches sklearn despite the eigenvector ambiguity); `intercept_ = -½·diag(means_@coef_.T) + log(priors_)` (`:483-485`). Supports `shrinkage` (like lsqr); `transform` is the un-centered `X @ scalings_[:, :max_components]` (eigen has NO `xbar_`, `:687`). Consumer: `fn fit` reads `self.solver` and dispatches; `Predict`/`predict_proba` for `FittedLDA` consume the eigen `coef_`/`intercept_`; `Transform` consumes `scalings_`. Tests `lda_eigen_solver` (collapsed `coef_[1]-coef_[0]` = `[14.7368…, 14.7368…]`, `explained_variance_ratio_` = `[1.0]`, predict/predict_proba) and `lda_eigen_shrinkage` (`Fixed(0.5)` collapsed coef = `[12.043…, 12.043…]`) <1e-6 vs live `LinearDiscriminantAnalysis(solver='eigen').fit(X,y)`. #596. |
//! | REQ-11 (shrinkage None/auto/float) | SHIPPED | `Shrinkage::{None, Auto, Fixed(F)}` (`LDA::with_shrinkage`) drives `fn cov_shrunk` (sklearn `_cov`, `discriminant_analysis.py:36-93`) inside `fn solve_lstsq`: `None` → maximum-likelihood empirical covariance (`fn empirical_covariance`, `np.cov(...,bias=1)`, `:76-77`); `Fixed(s)` → `(1-s)·emp + s·(trace(emp)/p)·I` (`shrunk_covariance`, `covariance/_shrunk_covariance.py:153-156`), validated `0 ≤ s ≤ 1` (`Interval(Real,0,1,closed=both)`, `:339`) else `InvalidParameter`; `Auto` → analytical Ledoit-Wolf (`fn ledoit_wolf_shrinkage`, transcribed from `covariance/_shrunk_covariance.py:365-401` unblocked case) on StandardScaler-standardized data then rescaled (`_cov` `:70-75`). `Solver::Svd` + non-`None` shrinkage → `InvalidParameter("shrinkage not supported with svd solver")` mirroring sklearn `NotImplementedError` (`:628-629`). Consumer: `fn fit`/`fn solve_lstsq` read `self.shrinkage`. Tests `lda_shrinkage_fixed` (`Fixed(0.5)` coef = `[12.043…, 12.043…]`), `lda_shrinkage_auto` (`Auto` coef = `[11.3706…, 11.3706…]`, validates the Ledoit-Wolf transcription), `lda_svd_shrinkage_rejected` (svd+shrinkage → `Err`) <1e-6 vs the live oracle. #597. |
//! | REQ-12 (store_covariance / covariance_) | SHIPPED | `LDA::with_store_covariance` sets the flag (sklearn default `false`, `discriminant_analysis.py:353`); when `true`, `fn fit` computes the shared within-class covariance `covariance_ = Σ_k priors_[k] · cov(X_k)` (`:509-510`, `_class_cov` `:128-172`) with the maximum-likelihood (`bias=1`, ÷`n_k`) per-class empirical covariance (`empirical_covariance`, `np.cov(...,bias=1)`), stored on `FittedLDA::covariance` (`None` when the flag is unset, matching sklearn). Consumer: `fn fit` reads `self.store_covariance`/`priors`/`means` and populates the field; `FittedLDA::covariance` exposes it. Test `lda_store_covariance` matches the live oracle `LinearDiscriminantAnalysis(store_covariance=True).fit(X,y).covariance_` to 1e-9 and asserts `None` for the default/`false` path. #598. |
//! | REQ-14 (binary decision_function shape `(n,)`) | NOT-STARTED | open prereq blocker #600. `fn decision_function` always returns `(n, n_classes)`; sklearn collapses binary to `(n,)` (`discriminant_analysis.py:651-657,739`). Binding-ABI layer (parallel to QDA #581). |
//! | REQ-15 (tol parameter) | SHIPPED | `LDA::with_tol` sets the svd-solver rank threshold (sklearn default `1e-4`, `discriminant_analysis.py:354,362`); `fn fit` reads `self.tol` into BOTH rank cutoffs `rank = Σ(S > tol)` (`:532`) and `rank2 = Σ(S2 > tol·S2[0])` (`:554`), REPLACING the prior hardcoded `1e-4`. Default `1e-4` ⇒ byte-identical to prior behavior (all existing svd-fit oracle tests stay green). Consumer: `fn fit` reads `self.tol` in both rank thresholds. Test `lda_tol_param` (field default `1e-4` + `with_tol` plumb-through). #601. |
//! | REQ-16 (ferray array-type substrate) | NOT-STARTED | open prereq blocker #602. The two SVDs run on `ferray::linalg::svd`; the owned array type is still `ndarray` (crate-wide deferral, cf. qda.rs REQ-12 #585). |
//! | REQ-17 (non-finite input rejected) | SHIPPED | The shared `fn fit` entry rejects any NaN/+/-inf in X BEFORE the solver dispatch (svd/lsqr/eigen) with `FerroError::InvalidParameter`, mirroring sklearn's `_validate_data(force_all_finite=True)` (`discriminant_analysis.py:589`) → `ValueError("Input X contains NaN.")` / `"... contains infinity ..."`. `y` is `Array1<usize>` (integer labels), finite by type; LDA's `fit` takes no `sample_weight`, so X is the only runtime check. All three solvers dispatch downstream of the guard, so it covers every solver path. `.iter().any(|v| !v.is_finite())` catches NaN and Inf; the finite path is byte-identical. Verified vs the live sklearn 1.5.2 oracle (R-CHAR-3): `LinearDiscriminantAnalysis().fit` raises `ValueError` for NaN/+inf/-inf in X (`tests/divergence_linear_nonfinite_batch4.rs::lda_*`). Non-test consumer: the existing `Fit for LDA` consumers (`Predict for FittedLDA`, crate-root `pub use`). (#2263) |
//!
//! # Examples
//!
//! ```
//! use ferrolearn_linear::lda::LDA;
//! use ferrolearn_core::{Fit, Predict};
//! use ndarray::{array, Array1, Array2};
//!
//! let lda = LDA::new(Some(1));
//! let x = Array2::from_shape_vec(
//!     (6, 2),
//!     vec![1.0, 1.0, 1.5, 1.2, 1.2, 0.8, 5.0, 5.0, 5.5, 4.8, 4.8, 5.2],
//! ).unwrap();
//! let y = array![0usize, 0, 0, 1, 1, 1];
//! let fitted = lda.fit(&x, &y).unwrap();
//! let preds = fitted.predict(&x).unwrap();
//! assert_eq!(preds.len(), 6);
//! ```

use ferray::linalg::{LinalgFloat, cholesky, eigh, inv, svd};
use ferray::{Array as FerrayArray, Ix2 as FerrayIx2};
use ferrolearn_core::error::FerroError;
use ferrolearn_core::introspection::HasClasses;
use ferrolearn_core::pipeline::{FittedPipelineEstimator, PipelineEstimator};
use ferrolearn_core::traits::{Fit, Predict, Transform};
use ndarray::{Array1, Array2, ScalarOperand};
use num_traits::{Float, NumCast};

// ---------------------------------------------------------------------------
// Solver / Shrinkage enums
// ---------------------------------------------------------------------------

/// LDA solver selector (sklearn's `solver`, `discriminant_analysis.py:204-216`,
/// `_parameter_constraints` `StrOptions({svd, lsqr, eigen})` `:338`).
///
/// - [`Solver::Svd`] (default) — the singular-value-decomposition path
///   (`_solve_svd`, `discriminant_analysis.py:487-559`); supports `transform`
///   (dimensionality reduction) but NOT `shrinkage`.
/// - [`Solver::Lsqr`] — the least-squares path (`_solve_lstsq`,
///   `discriminant_analysis.py:365-419`): `coef_ = lstsq(covariance_,
///   means_.T).T`, `intercept_ = -½·diag(means_ @ coef_.T) + log(priors_)`.
///   Supports `shrinkage`; does NOT support `transform` (sklearn raises
///   `NotImplementedError`, `:676-679`).
/// - [`Solver::Eigen`] — the generalized-eigenvalue path
///   (`_solve_eigen`, `discriminant_analysis.py:421-485`): forms `Sw`/`St`,
///   solves the generalized `eigh(Sb, Sw)` (reduced to a STANDARD symmetric
///   eigenproblem via the Cholesky factor of `Sw`), and yields `scalings_`,
///   `coef_`, `intercept_`, `explained_variance_ratio_`. Supports `shrinkage`
///   and `transform` (the un-centered `X @ scalings_`).
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum Solver {
    /// Singular-value-decomposition solver (sklearn default).
    #[default]
    Svd,
    /// Least-squares solver (`_solve_lstsq`).
    Lsqr,
    /// Generalized-eigenvalue solver (`_solve_eigen`), via Cholesky reduction.
    Eigen,
}

/// LDA covariance-shrinkage selector (sklearn's `shrinkage`,
/// `discriminant_analysis.py:218-225`, `_parameter_constraints`
/// `[StrOptions({auto}), Interval(Real, 0, 1, closed=both), None]` `:339`).
///
/// Drives the per-class covariance estimate `_cov`
/// (`discriminant_analysis.py:36-93`) inside the `lsqr`/`eigen` solvers
/// (sklearn note `:225`: shrinkage works only with the `lsqr` and `eigen`
/// solvers):
/// - [`Shrinkage::None`] — no shrinkage; the maximum-likelihood empirical
///   covariance (`np.cov(..., bias=1)`, `:76-77`).
/// - [`Shrinkage::Auto`] — automatic Ledoit-Wolf shrinkage (`:70-75`):
///   standardize features, run the Ledoit-Wolf lemma, then rescale.
/// - [`Shrinkage::Fixed`]`(s)` — fixed shrinkage `s ∈ [0, 1]`
///   (`shrunk_covariance`, `_shrunk_covariance.py:111-158`):
///   `(1 - s)·emp_cov + s·(trace(emp_cov)/p)·I`.
///
/// # Type Parameters
///
/// - `F`: the floating-point scalar type (`f32` or `f64`).
#[derive(Debug, Clone, Copy, PartialEq, Default)]
pub enum Shrinkage<F> {
    /// No shrinkage (sklearn `None`/`'empirical'`). Default.
    #[default]
    None,
    /// Automatic Ledoit-Wolf shrinkage (sklearn `'auto'`).
    Auto,
    /// Fixed shrinkage coefficient `s ∈ [0, 1]` (sklearn `float`).
    Fixed(F),
}

// ---------------------------------------------------------------------------
// LDA (unfitted)
// ---------------------------------------------------------------------------

/// Linear Discriminant Analysis configuration.
///
/// Holds hyperparameters. Calling [`Fit::fit`] runs sklearn's default
/// `solver="svd"` path (`discriminant_analysis.py:487-559`) and returns a
/// [`FittedLDA`].
///
/// Use [`LDA::with_priors`] to supply class priors (sklearn's `priors`,
/// `discriminant_analysis.py:359`); the default `None` infers the empirical
/// priors `n_k / n` at fit time.
///
/// # Type Parameters
///
/// - `F`: The floating-point scalar type (`f32` or `f64`).
#[derive(Debug, Clone)]
pub struct LDA<F> {
    /// Number of discriminant components to retain.
    ///
    /// If `None`, defaults to `min(n_classes - 1, n_features)` at fit time.
    n_components: Option<usize>,

    /// Class prior probabilities (sklearn's `priors`,
    /// `discriminant_analysis.py:351,359`).
    ///
    /// `None` (sklearn's default) ⇒ the empirical priors `n_k / n` are inferred
    /// from the training data at fit time. `Some(p)` ⇒ `p` is used VERBATIM as
    /// `priors_` (matching sklearn `:605`, `self.priors_ = xp.asarray(self.priors)`).
    priors: Option<Array1<F>>,

    /// Solver selector (sklearn's `solver`, default `"svd"`,
    /// `discriminant_analysis.py:204,349`). See [`Solver`].
    solver: Solver,

    /// Covariance-shrinkage selector (sklearn's `shrinkage`, default `None`,
    /// `discriminant_analysis.py:218,350`). See [`Shrinkage`]. Only honored by
    /// the `lsqr` solver here; combined with `Solver::Svd` it is rejected at fit
    /// (sklearn `NotImplementedError`, `:628-629`).
    shrinkage: Shrinkage<F>,

    /// Whether to compute and store the shared within-class covariance matrix
    /// `covariance_` during fit (sklearn's `store_covariance`, default `false`,
    /// `discriminant_analysis.py:353,361`). When `true`, the svd-solver `fit`
    /// computes `covariance_ = Σ_k priors_[k] · cov(X_k)` (`:509-510`,
    /// `_class_cov` `:128-172`).
    store_covariance: bool,

    /// Singular-value rank threshold used by the svd-solver (sklearn's `tol`,
    /// default `1e-4`, `discriminant_analysis.py:354,362`). It drives the two
    /// rank cutoffs `rank = Σ(S > tol)` (`:532`) and
    /// `rank2 = Σ(S2 > tol·S2[0])` (`:554`).
    tol: F,

    _marker: std::marker::PhantomData<F>,
}

impl<F: Float + Send + Sync + 'static> LDA<F> {
    /// Create a new `LDA`.
    ///
    /// - `n_components`: number of discriminant directions to retain.
    ///   Pass `None` to use `min(n_classes - 1, n_features)`.
    #[must_use]
    pub fn new(n_components: Option<usize>) -> Self {
        Self {
            n_components,
            priors: None,
            solver: Solver::Svd,
            shrinkage: Shrinkage::None,
            store_covariance: false,
            // sklearn default `tol=1e-4` (`discriminant_analysis.py:354`).
            // `1e-4` is exactly representable in f32/f64; the fallback to
            // `F::epsilon()` is unreachable for those but keeps `new`
            // panic-free for any conforming `Float`.
            tol: F::from(1e-4).unwrap_or_else(F::epsilon),
            _marker: std::marker::PhantomData,
        }
    }

    /// Return the configured number of components (may be `None`).
    #[must_use]
    pub fn n_components(&self) -> Option<usize> {
        self.n_components
    }

    /// Set the class prior probabilities (sklearn's `priors`,
    /// `discriminant_analysis.py:351,359`).
    ///
    /// The provided vector is used VERBATIM as `priors_` (sklearn does not
    /// normalize it here when it already sums to 1, `:605`). Its length must
    /// equal the number of classes seen at fit time, or [`Fit::fit`] returns
    /// [`FerroError::ShapeMismatch`]. Pass nothing (the `None` default) to infer
    /// the empirical priors `n_k / n` from the training data (`:601-603`).
    #[must_use]
    pub fn with_priors(mut self, priors: Array1<F>) -> Self {
        self.priors = Some(priors);
        self
    }

    /// Return the configured class priors (`None` ⇒ empirical at fit time).
    /// Mirrors sklearn's constructor `priors` (`discriminant_analysis.py:359`).
    #[must_use]
    pub fn priors(&self) -> Option<&Array1<F>> {
        self.priors.as_ref()
    }

    /// Set the solver (sklearn's `solver`, `discriminant_analysis.py:204,349`).
    /// Default [`Solver::Svd`]. See [`Solver`].
    ///
    /// [`Solver::Lsqr`] enables the least-squares path (and `shrinkage`);
    /// [`Solver::Eigen`] enables the generalized-eigenvalue path (also supports
    /// `shrinkage` and `transform`).
    #[must_use]
    pub fn with_solver(mut self, solver: Solver) -> Self {
        self.solver = solver;
        self
    }

    /// Return the configured solver. Mirrors sklearn's constructor `solver`
    /// (`discriminant_analysis.py:349`, default `"svd"`).
    #[must_use]
    pub fn solver(&self) -> Solver {
        self.solver
    }

    /// Set the covariance shrinkage (sklearn's `shrinkage`,
    /// `discriminant_analysis.py:218,350`). Default [`Shrinkage::None`]. See
    /// [`Shrinkage`].
    ///
    /// Honored only by [`Solver::Lsqr`] here (sklearn note `:225`: shrinkage
    /// works only with the `lsqr`/`eigen` solvers). Combined with
    /// [`Solver::Svd`], [`Fit::fit`] returns a [`FerroError`] mirroring sklearn's
    /// `NotImplementedError("shrinkage not supported with 'svd' solver.")`
    /// (`:628-629`). [`Shrinkage::Fixed`]`(s)` requires `0 <= s <= 1` (sklearn
    /// `Interval(Real, 0, 1, closed="both")`, `:339`), else [`Fit::fit`] returns
    /// [`FerroError::InvalidParameter`].
    #[must_use]
    pub fn with_shrinkage(mut self, shrinkage: Shrinkage<F>) -> Self {
        self.shrinkage = shrinkage;
        self
    }

    /// Return the configured covariance shrinkage. Mirrors sklearn's constructor
    /// `shrinkage` (`discriminant_analysis.py:350`, default `None`).
    #[must_use]
    pub fn shrinkage(&self) -> Shrinkage<F> {
        self.shrinkage
    }

    /// Set whether to compute and store the shared within-class covariance
    /// matrix `covariance_` during fit (sklearn's `store_covariance`,
    /// `discriminant_analysis.py:353,361`). Default `false`.
    ///
    /// When `true`, [`Fit::fit`] computes `covariance_ = Σ_k priors_[k] ·
    /// cov(X_k)` (`:509-510`, `_class_cov` `:128-172`) and
    /// [`FittedLDA::covariance`] returns `Some`. When `false` it returns `None`
    /// (matching sklearn, where the attribute only exists when the flag is set).
    #[must_use]
    pub fn with_store_covariance(mut self, store_covariance: bool) -> Self {
        self.store_covariance = store_covariance;
        self
    }

    /// Return whether `covariance_` will be stored during fit (sklearn's
    /// `store_covariance`, `discriminant_analysis.py:353`).
    #[must_use]
    pub fn store_covariance(&self) -> bool {
        self.store_covariance
    }

    /// Set the singular-value rank threshold `tol` used by the svd-solver
    /// (sklearn's `tol`, `discriminant_analysis.py:354,362`). Default `1e-4`.
    ///
    /// It drives the two rank cutoffs `rank = Σ(S > tol)` (`:532`) and
    /// `rank2 = Σ(S2 > tol·S2[0])` (`:554`).
    #[must_use]
    pub fn with_tol(mut self, tol: F) -> Self {
        self.tol = tol;
        self
    }

    /// Return the configured svd-solver rank threshold `tol`. Mirrors sklearn's
    /// constructor `tol` (`discriminant_analysis.py:354`, default `1e-4`).
    #[must_use]
    pub fn tol(&self) -> F {
        self.tol
    }
}

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

// ---------------------------------------------------------------------------
// FittedLDA
// ---------------------------------------------------------------------------

/// A fitted LDA model (sklearn's `svd` solver).
///
/// Created by calling [`Fit::fit`] on an [`LDA`]. Implements:
/// - [`Transform<Array2<F>>`] — project data via `(X - xbar_) @ scalings_`.
/// - [`Predict<Array2<F>>`] — classify by argmax of the affine
///   `decision_function`.
#[derive(Debug, Clone)]
pub struct FittedLDA<F> {
    /// Whitened projection matrix `scalings_`, shape `(n_features, rank2)`.
    /// Mirrors sklearn's `scalings_` (`discriminant_analysis.py:555`).
    scalings: Array2<F>,

    /// Per-class means in the ORIGINAL feature space, shape
    /// `(n_classes, n_features)`. Mirrors sklearn's `means_`
    /// (`discriminant_analysis.py:508`).
    means: Array2<F>,

    /// Weighted overall mean `xbar_ = priors_ @ means_`, length `n_features`.
    /// Mirrors sklearn's `xbar_` (`discriminant_analysis.py:517`).
    xbar: Array1<F>,

    /// Resolved class priors `priors_`, length `n_classes`. Empirical `n_k / n`
    /// when the constructor `priors` was `None`, else the provided `priors`
    /// verbatim. Mirrors sklearn's `priors_` (`discriminant_analysis.py:601-605`).
    priors: Array1<F>,

    /// Affine classifier coefficients `coef_`, shape `(n_classes, n_features)`.
    /// Mirrors sklearn's `coef_` (`discriminant_analysis.py:558`). (Binary
    /// collapse to `(1, n_features)` pends #600.)
    coef: Array2<F>,

    /// Affine classifier intercepts `intercept_`, length `n_classes` (embeds
    /// `log(priors_)`). Mirrors sklearn's `intercept_`
    /// (`discriminant_analysis.py:557,559`).
    intercept: Array1<F>,

    /// Ratio of explained variance per discriminant direction, length
    /// `max_components`. Mirrors sklearn's `explained_variance_ratio_`
    /// (`discriminant_analysis.py:550-552`).
    explained_variance_ratio: Array1<F>,

    /// Shared within-class covariance matrix `covariance_`, shape
    /// `(n_features, n_features)`, present only when the model was configured
    /// with [`LDA::with_store_covariance`]`(true)`. Mirrors sklearn's
    /// `covariance_` (`discriminant_analysis.py:509-510`, `_class_cov`
    /// `:128-172`): `Σ_k priors_[k] · cov(X_k)`. `None` otherwise (matching
    /// sklearn, where the attribute only exists when `store_covariance=True`).
    covariance: Option<Array2<F>>,

    /// Class labels corresponding to rows of `means`/`coef`.
    classes: Vec<usize>,

    /// Number of components to keep on `transform` output (sklearn's
    /// `_max_components`, `discriminant_analysis.py:619/625`).
    max_components: usize,

    /// Number of features seen during fitting.
    n_features: usize,
}

impl<F: Float + Send + Sync + 'static> FittedLDA<F> {
    /// Whitened projection (`scalings_`) matrix, shape `(n_features, rank2)`.
    /// Mirrors sklearn's `scalings_` (`discriminant_analysis.py:555`).
    #[must_use]
    pub fn scalings(&self) -> &Array2<F> {
        &self.scalings
    }

    /// Per-class means in the original feature space, shape
    /// `(n_classes, n_features)`. Mirrors sklearn's `means_`
    /// (`discriminant_analysis.py:508`).
    #[must_use]
    pub fn means(&self) -> &Array2<F> {
        &self.means
    }

    /// Weighted overall mean `xbar_`, length `n_features`. Mirrors sklearn's
    /// `xbar_` (`discriminant_analysis.py:517`).
    #[must_use]
    pub fn xbar(&self) -> &Array1<F> {
        &self.xbar
    }

    /// Resolved class priors `priors_`, length `n_classes` (empirical `n_k / n`
    /// when the constructor `priors` was `None`, else the provided `priors`
    /// verbatim). Mirrors sklearn's `priors_` (`discriminant_analysis.py:601-605`).
    #[must_use]
    pub fn priors(&self) -> &Array1<F> {
        &self.priors
    }

    /// Affine classifier coefficients `coef_`, shape `(n_classes, n_features)`.
    /// Mirrors sklearn's `coef_` (`discriminant_analysis.py:558`).
    #[must_use]
    pub fn coef(&self) -> &Array2<F> {
        &self.coef
    }

    /// Affine classifier intercepts `intercept_`, length `n_classes`. Mirrors
    /// sklearn's `intercept_` (`discriminant_analysis.py:557,559`).
    #[must_use]
    pub fn intercept(&self) -> &Array1<F> {
        &self.intercept
    }

    /// Explained-variance ratio per discriminant direction. Mirrors sklearn's
    /// `explained_variance_ratio_` (`discriminant_analysis.py:550-552`).
    #[must_use]
    pub fn explained_variance_ratio(&self) -> &Array1<F> {
        &self.explained_variance_ratio
    }

    /// Shared within-class covariance matrix `covariance_`, shape
    /// `(n_features, n_features)`. Mirrors sklearn's `covariance_`
    /// (`discriminant_analysis.py:509-510`, `_class_cov` `:128-172`):
    /// `Σ_k priors_[k] · cov(X_k)` where `cov(X_k)` is the maximum-likelihood
    /// empirical covariance of class `k`'s samples (`np.cov(..., bias=1)`,
    /// normalized by `n_k`, via `empirical_covariance`,
    /// `covariance/_empirical_covariance.py:109`).
    ///
    /// Returns `Some` only when the model was configured with
    /// [`LDA::with_store_covariance`]`(true)`; `None` otherwise — matching
    /// sklearn, where the `covariance_` attribute only exists when
    /// `store_covariance=True`.
    #[must_use]
    pub fn covariance(&self) -> Option<&Array2<F>> {
        self.covariance.as_ref()
    }

    /// Sorted class labels as seen during fitting.
    #[must_use]
    pub fn classes(&self) -> &[usize] {
        &self.classes
    }

    /// Per-class discriminant scores. Mirrors sklearn
    /// `LinearDiscriminantAnalysis.decision_function` (the `LinearClassifierMixin`,
    /// `discriminant_analysis.py:739`): the affine map `X @ coef_.T + intercept_`.
    ///
    /// Returns shape `(n_samples, n_classes)`. (Binary collapse to `(n,)` pends
    /// REQ-14/#600.) argmax of each row agrees with [`Predict`].
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if the number of features does not
    /// match the fitted model.
    pub fn decision_function(&self, x: &Array2<F>) -> Result<Array2<F>, FerroError> {
        if x.ncols() != self.n_features {
            return Err(FerroError::ShapeMismatch {
                expected: vec![x.nrows(), self.n_features],
                actual: vec![x.nrows(), x.ncols()],
                context: "FittedLDA::decision_function".into(),
            });
        }
        // X @ coef_.T + intercept_  (coef_ is (n_classes, n_features)).
        let mut out = x.dot(&self.coef.t());
        let n_classes = self.intercept.len();
        for mut row in out.rows_mut() {
            for c in 0..n_classes {
                row[c] = row[c] + self.intercept[c];
            }
        }
        Ok(out)
    }

    /// Predict per-class probabilities. Mirrors sklearn
    /// `LinearDiscriminantAnalysis.predict_proba` (`discriminant_analysis.py:706-711`):
    /// the multiclass `softmax(decision_function)` (the row-max-shifted softmax
    /// of `sklearn.utils.extmath.softmax`, `extmath.py:949-985`).
    ///
    /// Returns shape `(n_samples, n_classes)`; rows sum to 1. (The binary
    /// `[1-expit(d), expit(d)]` collapse pends REQ-14/#600; the multiclass
    /// softmax here is correct for `n_classes >= 2` because `coef_`/`intercept_`
    /// are not yet collapsed to the binary single-row form.)
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if the number of features does not
    /// match the model.
    pub fn predict_proba(&self, x: &Array2<F>) -> Result<Array2<F>, FerroError> {
        let decision = self.decision_function(x)?;
        let n_samples = decision.nrows();
        let n_classes = decision.ncols();
        let mut proba = Array2::<F>::zeros((n_samples, n_classes));
        for i in 0..n_samples {
            let max_l = (0..n_classes)
                .map(|c| decision[[i, c]])
                .fold(F::neg_infinity(), |a, b| if b > a { b } else { a });
            let mut sum_exp = F::zero();
            for c in 0..n_classes {
                let e = (decision[[i, c]] - max_l).exp();
                proba[[i, c]] = e;
                sum_exp = sum_exp + e;
            }
            for c in 0..n_classes {
                proba[[i, c]] = proba[[i, c]] / sum_exp;
            }
        }
        Ok(proba)
    }

    /// Element-wise log of [`predict_proba`](Self::predict_proba). Mirrors
    /// sklearn `predict_log_proba` exactly (`discriminant_analysis.py:713-737`):
    /// entries that are EXACTLY `0.0` are bumped by the dtype's
    /// `smallest_normal` (`f32`/`f64::MIN_POSITIVE`) before taking `log`
    /// (`:729-736`), so `log(0)` becomes `log(MIN_POSITIVE)` rather than `-inf`;
    /// every nonzero probability keeps its true `ln`.
    ///
    /// # Errors
    ///
    /// Forwards any error from [`predict_proba`](Self::predict_proba).
    pub fn predict_log_proba(&self, x: &Array2<F>) -> Result<Array2<F>, FerroError> {
        let proba = self.predict_proba(x)?;
        // sklearn: prediction[prediction == 0.0] += smallest_normal; log(prediction).
        // `F::min_positive_value()` is numpy's `finfo(dtype).smallest_normal`
        // (`f64::MIN_POSITIVE` ≈ 2.2250738585072014e-308).
        let smallest_normal = F::min_positive_value();
        Ok(proba.mapv(|p| {
            if p == F::zero() {
                (p + smallest_normal).ln()
            } else {
                p.ln()
            }
        }))
    }
}

// ---------------------------------------------------------------------------
// Internal helpers
// ---------------------------------------------------------------------------

/// Convert a `usize` count to `F` without panicking. Returns
/// [`FerroError::NumericalInstability`] if the value is not representable.
#[inline]
fn usize_to_f<F: Float>(v: usize) -> Result<F, FerroError> {
    F::from(v).ok_or_else(|| FerroError::NumericalInstability {
        message: format!("could not represent count {v} as the float type"),
    })
}

/// `0.5` as `F`, built panic-free from `1 / (1 + 1)` (exact for binary floats).
#[inline]
fn half<F: Float>() -> F {
    F::one() / (F::one() + F::one())
}

/// Singular values `S` and right singular vectors transposed `Vt` of the thin
/// SVD `A = U·diag(S)·Vt` (`full_matrices=False`), on the ferray substrate
/// (`ferray::linalg::svd`, the analog of `scipy.linalg.svd(X,
/// full_matrices=False)`, `discriminant_analysis.py:530,545`). Mirrors the
/// bridging pattern in `qda.rs::svd_s_vt` / `bayesian_ridge.rs::svd_thin`
/// (R-SUBSTRATE-4): the caller keeps its `ndarray` signature and the
/// ndarray↔ferray conversion happens here.
///
/// Returns `(S, Vt)` with `S` of length `k = min(m, n)` (descending) and `Vt`
/// of shape `(k, n)`.
///
/// # Errors
///
/// Returns [`FerroError::NumericalInstability`] if the ferray array build or
/// the SVD itself fails.
fn svd_s_vt<F: LinalgFloat>(a: &Array2<F>) -> Result<(Array1<F>, Array2<F>), FerroError> {
    let (m, n) = a.dim();
    let a_flat: Vec<F> = a.iter().copied().collect();
    let fa =
        FerrayArray::<F, FerrayIx2>::from_vec(FerrayIx2::new([m, n]), a_flat).map_err(|e| {
            FerroError::NumericalInstability {
                message: format!("ferray svd: failed to build matrix: {e}"),
            }
        })?;
    let (_u, s, vt) = svd(&fa, false).map_err(|e| FerroError::NumericalInstability {
        message: format!("ferray svd failed: {e}"),
    })?;
    let s_nd = Array1::from_vec(s.iter().copied().collect());
    let vt_shape = vt.shape();
    let vt_nd = Array2::from_shape_vec((vt_shape[0], vt_shape[1]), vt.iter().copied().collect())
        .map_err(|e| FerroError::NumericalInstability {
            message: format!("ferray svd: Vt shape conversion failed: {e}"),
        })?;
    Ok((s_nd, vt_nd))
}

/// Maximum-likelihood empirical covariance of `xg` (rows = samples), centered on
/// the per-column mean and normalized by `n` (NOT `n-1`). Mirrors sklearn's
/// `empirical_covariance` / `np.cov(..., bias=1)` (`discriminant_analysis.py:77`,
/// `covariance/_empirical_covariance.py:109`).
///
/// Returns the `(p, p)` covariance and the per-column means (length `p`).
fn empirical_covariance<F: Float>(xg: &Array2<F>) -> Result<(Array2<F>, Array1<F>), FerroError> {
    let (n, p) = xg.dim();
    let nf = usize_to_f::<F>(n)?;
    let mut mean = Array1::<F>::zeros(p);
    for i in 0..n {
        for j in 0..p {
            mean[j] = mean[j] + xg[[i, j]];
        }
    }
    for j in 0..p {
        mean[j] = mean[j] / nf;
    }
    let mut cov = Array2::<F>::zeros((p, p));
    for a in 0..p {
        for b in 0..p {
            let mut acc = F::zero();
            for i in 0..n {
                acc = acc + (xg[[i, a]] - mean[a]) * (xg[[i, b]] - mean[b]);
            }
            cov[[a, b]] = acc / nf;
        }
    }
    Ok((cov, mean))
}

/// Ledoit-Wolf analytical shrinkage coefficient of `x` (rows = samples), the
/// transcription of sklearn's `ledoit_wolf_shrinkage`
/// (`covariance/_shrunk_covariance.py:299-402`) for the unblocked case
/// (`block_size=1000 >> n_features`, so `n_splits = 0` and the blocked loops
/// collapse to the single tail term `beta_ = Σ(X²ᵀ·X²)`, `delta_ =
/// Σ((Xᵀ·X)²)`). `x` is assumed already centered (`assume_centered=True`, the
/// caller centers in [`cov_shrunk`]).
///
/// Formula (`:365-401`): with `X² = x⊙x`, `emp_cov_trace = Σ_i X²[i,:]/n`,
/// `mu = Σ(emp_cov_trace)/p`, `beta_ = Σ(X²ᵀ·X²)`, `delta_ = Σ((Xᵀ·X)²)/n²`,
/// `beta = (1/(p·n))·(beta_/n − delta_)`,
/// `delta = (delta_ − 2·mu·Σ(emp_cov_trace) + p·mu²)/p`,
/// `beta = min(beta, delta)`, `shrinkage = 0 if beta==0 else beta/delta`.
fn ledoit_wolf_shrinkage<F: Float>(x: &Array2<F>) -> Result<F, FerroError> {
    let (n, p) = x.dim();
    // sklearn `:345-346`: for a single feature the result is shrinkage-invariant.
    if p == 1 {
        return Ok(F::zero());
    }
    let nf = usize_to_f::<F>(n)?;
    let pf = usize_to_f::<F>(p)?;

    // emp_cov_trace[j] = Σ_i x[i,j]² / n   (:365-366)
    let mut emp_cov_trace = Array1::<F>::zeros(p);
    for j in 0..p {
        let mut acc = F::zero();
        for i in 0..n {
            acc = acc + x[[i, j]] * x[[i, j]];
        }
        emp_cov_trace[j] = acc / nf;
    }
    // mu = Σ(emp_cov_trace) / p   (:367)
    let mut trace_sum = F::zero();
    for j in 0..p {
        trace_sum = trace_sum + emp_cov_trace[j];
    }
    let mu = trace_sum / pf;

    // beta_ = Σ over (a,b) of (X²ᵀ·X²)[a,b] = Σ_{a,b} Σ_i x[i,a]²·x[i,b]²  (:388-390)
    // delta_ = Σ over (a,b) of ((Xᵀ·X)[a,b])²  (:384-386)
    let mut beta_acc = F::zero();
    let mut delta_acc = F::zero();
    for a in 0..p {
        for b in 0..p {
            let mut g_ab = F::zero(); // (Xᵀ·X)[a,b] = Σ_i x[i,a]·x[i,b]
            let mut h_ab = F::zero(); // (X²ᵀ·X²)[a,b] = Σ_i x[i,a]²·x[i,b]²
            for i in 0..n {
                let xa = x[[i, a]];
                let xb = x[[i, b]];
                g_ab = g_ab + xa * xb;
                h_ab = h_ab + (xa * xa) * (xb * xb);
            }
            beta_acc = beta_acc + h_ab;
            delta_acc = delta_acc + g_ab * g_ab;
        }
    }
    // delta_ /= n²   (:387)
    let delta_ = delta_acc / (nf * nf);
    // beta = (1/(p·n)) · (beta_/n − delta_)   (:392)
    let beta = (F::one() / (pf * nf)) * (beta_acc / nf - delta_);
    // delta = (delta_ − 2·mu·Σ(emp_cov_trace) + p·mu²) / p   (:394-395)
    let two = F::one() + F::one();
    let mut delta = delta_ - two * mu * trace_sum + pf * mu * mu;
    delta = delta / pf;
    // beta = min(beta, delta)   (:399)
    let beta = if beta < delta { beta } else { delta };
    // shrinkage = 0 if beta==0 else beta/delta   (:401)
    if beta == F::zero() {
        Ok(F::zero())
    } else {
        Ok(beta / delta)
    }
}

/// Per-class covariance estimate with optional shrinkage — sklearn's `_cov`
/// (`discriminant_analysis.py:36-93`) for `covariance_estimator=None`:
/// - [`Shrinkage::None`] → empirical maximum-likelihood covariance (`:76-77`).
/// - [`Shrinkage::Fixed`]`(s)` → `shrunk_covariance(emp_cov, s)`
///   (`:78-79`, `_shrunk_covariance.py:153-156`):
///   `(1 − s)·emp_cov + s·(trace(emp_cov)/p)·I`.
/// - [`Shrinkage::Auto`] → Ledoit-Wolf on the StandardScaler-standardized data,
///   then rescaled (`:70-75`): standardize `Xs = (X − mean)/scale` (population
///   std, `ddof=0`, zeros → 1), `s = ledoit_wolf(Xs)`, then `cov[a,b] =
///   scale[a]·s[a,b]·scale[b]`.
fn cov_shrunk<F: Float>(xg: &Array2<F>, shrinkage: Shrinkage<F>) -> Result<Array2<F>, FerroError> {
    let (n, p) = xg.dim();
    match shrinkage {
        Shrinkage::None => {
            let (cov, _mean) = empirical_covariance(xg)?;
            Ok(cov)
        }
        Shrinkage::Fixed(s) => {
            let (emp, _mean) = empirical_covariance(xg)?;
            let pf = usize_to_f::<F>(p)?;
            let mut trace = F::zero();
            for j in 0..p {
                trace = trace + emp[[j, j]];
            }
            let mu = trace / pf;
            let mut out = Array2::<F>::zeros((p, p));
            for a in 0..p {
                for b in 0..p {
                    let diag = if a == b { mu } else { F::zero() };
                    out[[a, b]] = (F::one() - s) * emp[[a, b]] + s * diag;
                }
            }
            Ok(out)
        }
        Shrinkage::Auto => {
            // StandardScaler: center + divide by POPULATION std (ddof=0); zeros
            // replaced by 1.0 (sklearn StandardScaler `_handle_zeros_in_scale`).
            let nf = usize_to_f::<F>(n)?;
            let mut mean = Array1::<F>::zeros(p);
            for i in 0..n {
                for j in 0..p {
                    mean[j] = mean[j] + xg[[i, j]];
                }
            }
            for j in 0..p {
                mean[j] = mean[j] / nf;
            }
            let mut scale = Array1::<F>::zeros(p);
            for j in 0..p {
                let mut var = F::zero();
                for i in 0..n {
                    let d = xg[[i, j]] - mean[j];
                    var = var + d * d;
                }
                var = var / nf;
                let sd = var.sqrt();
                scale[j] = if sd == F::zero() { F::one() } else { sd };
            }
            let mut xs = Array2::<F>::zeros((n, p));
            for i in 0..n {
                for j in 0..p {
                    xs[[i, j]] = (xg[[i, j]] - mean[j]) / scale[j];
                }
            }
            // ledoit_wolf re-centers; Xs already has ~0 mean, but follow sklearn
            // exactly and re-center (assume_centered=False, `:357-358`).
            let mut xs_mean = Array1::<F>::zeros(p);
            for i in 0..n {
                for j in 0..p {
                    xs_mean[j] = xs_mean[j] + xs[[i, j]];
                }
            }
            for j in 0..p {
                xs_mean[j] = xs_mean[j] / nf;
            }
            let mut xc = Array2::<F>::zeros((n, p));
            for i in 0..n {
                for j in 0..p {
                    xc[[i, j]] = xs[[i, j]] - xs_mean[j];
                }
            }
            // shrinkage coefficient from the (centered) standardized data.
            let shr = ledoit_wolf_shrinkage(&xc)?;
            // emp_cov of the standardized data = Xcᵀ·Xc / n, then shrink:
            // (1 − shr)·emp + shr·(trace(emp)/p)·I  (`_shrunk_covariance.py`).
            let (emp, _m) = empirical_covariance(&xs)?;
            let pf = usize_to_f::<F>(p)?;
            let mut trace = F::zero();
            for j in 0..p {
                trace = trace + emp[[j, j]];
            }
            let mu = trace / pf;
            // rescale: cov[a,b] = scale[a] · shrunk[a,b] · scale[b]  (`:75`).
            let mut out = Array2::<F>::zeros((p, p));
            for a in 0..p {
                for b in 0..p {
                    let diag = if a == b { mu } else { F::zero() };
                    let shrunk = (F::one() - shr) * emp[[a, b]] + shr * diag;
                    out[[a, b]] = scale[a] * shrunk * scale[b];
                }
            }
            Ok(out)
        }
    }
}

/// Solve the multi-RHS least-squares problem `A @ x = b` (with `b` having
/// `nrhs` columns) through [`ferray::linalg::lstsq`]
/// (`ferray-linalg/src/solve.rs:208`, the LAPACK-`gelsd`-equivalent single-SVD
/// min-norm solver), bridging ndarray↔ferray at this boundary (R-SUBSTRATE-4),
/// mirroring the bridge in `linalg.rs::solve_lstsq`. Returns the `(n, nrhs)`
/// solution. Used by `_solve_lstsq` to compute `coef_ = lstsq(covariance_,
/// means_.T)[0].T` (`discriminant_analysis.py:416`).
///
/// `rcond` is `Some(F::epsilon())`, pinning the singular-value cutoff to scipy's
/// `cond=eps` default (matching `linalg.lstsq` `:416`), as in `solve_lstsq`.
fn lstsq_multi<F: LinalgFloat>(a: &Array2<F>, b: &Array2<F>) -> Result<Array2<F>, FerroError> {
    let (m, n) = a.dim();
    let (bm, nrhs) = b.dim();
    if bm != m {
        return Err(FerroError::ShapeMismatch {
            expected: vec![m, nrhs],
            actual: vec![bm, nrhs],
            context: "LDA lsqr: covariance/means row mismatch".into(),
        });
    }
    let a_flat: Vec<F> = a.iter().copied().collect();
    let fa =
        FerrayArray::<F, FerrayIx2>::from_vec(FerrayIx2::new([m, n]), a_flat).map_err(|e| {
            FerroError::NumericalInstability {
                message: format!("ferray lstsq: failed to build matrix A: {e}"),
            }
        })?;
    let b_flat: Vec<F> = b.iter().copied().collect();
    let fb = FerrayArray::<F, ferray::IxDyn>::from_vec(ferray::IxDyn::new(&[bm, nrhs]), b_flat)
        .map_err(|e| FerroError::NumericalInstability {
            message: format!("ferray lstsq: failed to build RHS B: {e}"),
        })?;
    let (sol, _residuals, _rank, _singular) = ferray::linalg::lstsq(&fa, &fb, Some(F::epsilon()))
        .map_err(|e| FerroError::NumericalInstability {
        message: format!("ferray lstsq solve failed: {e}"),
    })?;
    let sol_shape = sol.shape();
    let out = Array2::from_shape_vec((sol_shape[0], sol_shape[1]), sol.iter().copied().collect())
        .map_err(|e| FerroError::NumericalInstability {
        message: format!("ferray lstsq: solution shape conversion failed: {e}"),
    })?;
    Ok(out)
}

/// Lower-triangular Cholesky factor `L` of the symmetric-positive-definite `a`
/// (`a = L·Lᵀ`), on the ferray substrate ([`ferray::linalg::cholesky`],
/// `ferray-linalg/src/decomp/cholesky.rs:22`, the analog of
/// `scipy.linalg.cholesky(..., lower=True)`), bridging ndarray↔ferray at this
/// boundary (R-SUBSTRATE-4) exactly as [`svd_s_vt`]/[`lstsq_multi`] do. Used by
/// the generalized-eigen reduction in `_solve_eigen`
/// (`discriminant_analysis.py:475`, `linalg.eigh(Sb, Sw)`).
///
/// # Errors
///
/// Returns [`FerroError::NumericalInstability`] if the ferray build or the
/// factorization fails (e.g. `Sw` is not positive definite).
fn cholesky_lower<F: LinalgFloat>(a: &Array2<F>) -> Result<Array2<F>, FerroError> {
    let (m, n) = a.dim();
    let a_flat: Vec<F> = a.iter().copied().collect();
    let fa =
        FerrayArray::<F, FerrayIx2>::from_vec(FerrayIx2::new([m, n]), a_flat).map_err(|e| {
            FerroError::NumericalInstability {
                message: format!("ferray cholesky: failed to build matrix: {e}"),
            }
        })?;
    let l = cholesky(&fa).map_err(|e| FerroError::NumericalInstability {
        message: format!("ferray cholesky failed (Sw not positive definite?): {e}"),
    })?;
    let shape = l.shape();
    Array2::from_shape_vec((shape[0], shape[1]), l.iter().copied().collect()).map_err(|e| {
        FerroError::NumericalInstability {
            message: format!("ferray cholesky: shape conversion failed: {e}"),
        }
    })
}

/// Inverse of the square matrix `a`, on the ferray substrate
/// ([`ferray::linalg::inv`], `ferray-linalg/src/solve.rs:367`, the analog of
/// `numpy.linalg.inv`), bridging ndarray↔ferray at this boundary
/// (R-SUBSTRATE-4). Used to form `L⁻¹` for the generalized-eigen reduction in
/// `_solve_eigen`.
///
/// # Errors
///
/// Returns [`FerroError::NumericalInstability`] if the ferray build or the
/// inversion fails (singular matrix).
fn matrix_inverse<F: LinalgFloat>(a: &Array2<F>) -> Result<Array2<F>, FerroError> {
    let (m, n) = a.dim();
    let a_flat: Vec<F> = a.iter().copied().collect();
    let fa =
        FerrayArray::<F, FerrayIx2>::from_vec(FerrayIx2::new([m, n]), a_flat).map_err(|e| {
            FerroError::NumericalInstability {
                message: format!("ferray inv: failed to build matrix: {e}"),
            }
        })?;
    let ai = inv(&fa).map_err(|e| FerroError::NumericalInstability {
        message: format!("ferray inv failed (singular matrix?): {e}"),
    })?;
    let shape = ai.shape();
    Array2::from_shape_vec((shape[0], shape[1]), ai.iter().copied().collect()).map_err(|e| {
        FerroError::NumericalInstability {
            message: format!("ferray inv: shape conversion failed: {e}"),
        }
    })
}

/// Eigenvalues and eigenvectors of the symmetric matrix `a`, on the ferray
/// substrate ([`ferray::linalg::eigh`], `ferray-linalg/src/decomp/eigen.rs:105`,
/// the analog of `scipy.linalg.eigh` for the STANDARD symmetric problem),
/// bridging ndarray↔ferray at this boundary (R-SUBSTRATE-4). Returns
/// `(evals, evecs)` with eigenvalues in ASCENDING order (ferray/LAPACK
/// convention) and eigenvectors as columns of `evecs`. The generalized
/// `eigh(Sb, Sw)` of `_solve_eigen` (`discriminant_analysis.py:475`) is reduced
/// to this standard form via the Cholesky factor of `Sw` (see `solve_eigen`).
///
/// # Errors
///
/// Returns [`FerroError::NumericalInstability`] if the ferray build or the
/// eigendecomposition fails.
fn eigh_sym<F: LinalgFloat>(a: &Array2<F>) -> Result<(Array1<F>, Array2<F>), FerroError> {
    let (m, n) = a.dim();
    let a_flat: Vec<F> = a.iter().copied().collect();
    let fa =
        FerrayArray::<F, FerrayIx2>::from_vec(FerrayIx2::new([m, n]), a_flat).map_err(|e| {
            FerroError::NumericalInstability {
                message: format!("ferray eigh: failed to build matrix: {e}"),
            }
        })?;
    let (vals, vecs) = eigh(&fa).map_err(|e| FerroError::NumericalInstability {
        message: format!("ferray eigh failed: {e}"),
    })?;
    let vals_nd = Array1::from_vec(vals.iter().copied().collect());
    let vecs_shape = vecs.shape();
    let vecs_nd = Array2::from_shape_vec(
        (vecs_shape[0], vecs_shape[1]),
        vecs.iter().copied().collect(),
    )
    .map_err(|e| FerroError::NumericalInstability {
        message: format!("ferray eigh: eigenvector shape conversion failed: {e}"),
    })?;
    Ok((vals_nd, vecs_nd))
}

// ---------------------------------------------------------------------------
// Fit (sklearn _solve_svd)
// ---------------------------------------------------------------------------

impl<F: LinalgFloat + ScalarOperand> Fit<Array2<F>, Array1<usize>> for LDA<F> {
    type Fitted = FittedLDA<F>;
    type Error = FerroError;

    /// Fit the LDA model via sklearn's default `solver="svd"` path
    /// (`discriminant_analysis.py:487-559`): two SVDs whiten the within-class
    /// data and project onto the between-class subspace, yielding `scalings_`,
    /// `xbar_`, `coef_`, `intercept_` (embedding `log(priors_)`), and
    /// `explained_variance_ratio_`.
    ///
    /// # Errors
    ///
    /// - [`FerroError::InsufficientSamples`] if fewer than 2 samples / classes.
    /// - [`FerroError::InvalidParameter`] if `n_components` is zero or exceeds
    ///   `min(n_classes - 1, n_features)`.
    /// - [`FerroError::ShapeMismatch`] if `x` and `y` have different row counts.
    /// - [`FerroError::NumericalInstability`] if an SVD fails.
    #[allow(
        clippy::needless_range_loop,
        reason = "explicit index loops mirror sklearn's broadcasting per-column/per-class"
    )]
    fn fit(&self, x: &Array2<F>, y: &Array1<usize>) -> Result<FittedLDA<F>, FerroError> {
        let (n_samples, n_features) = x.dim();

        if n_samples != y.len() {
            return Err(FerroError::ShapeMismatch {
                expected: vec![n_samples],
                actual: vec![y.len()],
                context: "LDA: y length must match number of rows in X".into(),
            });
        }
        if n_samples < 2 {
            return Err(FerroError::InsufficientSamples {
                required: 2,
                actual: n_samples,
                context: "LDA requires at least 2 samples".into(),
            });
        }

        // Sorted unique classes (sklearn `classes_ = unique_labels(y)`, :592).
        let mut classes: Vec<usize> = y.to_vec();
        classes.sort_unstable();
        classes.dedup();
        let n_classes = classes.len();

        if n_classes < 2 {
            return Err(FerroError::InsufficientSamples {
                required: 2,
                actual: n_classes,
                context: "LDA requires at least 2 distinct classes".into(),
            });
        }
        // sklearn rejects n_samples == n_classes (:596-599).
        if n_samples == n_classes {
            return Err(FerroError::InsufficientSamples {
                required: n_classes + 1,
                actual: n_samples,
                context: "LDA: number of samples must exceed number of classes".into(),
            });
        }

        // Non-finite input validation (#2263). sklearn `LinearDiscriminantAnalysis.fit`
        // -> `self._validate_data(X, y, ensure_min_samples=2, ...)`
        // (`discriminant_analysis.py:589`) keeps the default
        // `force_all_finite=True`, so `check_array` rejects any NaN or +/-inf in
        // X with a `ValueError("Input X contains NaN.")` / `"... contains
        // infinity ..."` BEFORE the solver dispatch (svd/lsqr/eigen). `y` is
        // `Array1<usize>` here (integer class labels), finite by construction, so
        // only X needs the runtime check; LDA's `fit` takes no `sample_weight`.
        // `.iter().any(|v| !v.is_finite())` rejects both NaN and Inf (bounds-safe,
        // no panic, R-CODE-2). This is the shared fit entry — all three solvers
        // (svd default, lsqr, eigen) dispatch downstream of it, so the guard
        // covers every solver. The finite path is byte-identical.
        if x.iter().any(|v| !v.is_finite()) {
            return Err(FerroError::InvalidParameter {
                name: "X".into(),
                reason: "Input X contains NaN or infinity.".into(),
            });
        }

        // _max_components (sklearn :614-625).
        let max_components = (n_classes - 1).min(n_features);
        let user_max = match self.n_components {
            None => max_components,
            Some(0) => {
                return Err(FerroError::InvalidParameter {
                    name: "n_components".into(),
                    reason: "must be at least 1".into(),
                });
            }
            Some(k) if k > max_components => {
                return Err(FerroError::InvalidParameter {
                    name: "n_components".into(),
                    reason: format!(
                        "n_components ({k}) cannot be larger than min(n_features, n_classes - 1) = {max_components}"
                    ),
                });
            }
            Some(k) => k,
        };

        let n_f = usize_to_f::<F>(n_samples)?;

        // --- per-class means_ and class indices (sklearn `_class_means`) ------
        let mut means = Array2::<F>::zeros((n_classes, n_features));
        let mut class_indices: Vec<Vec<usize>> = vec![Vec::new(); n_classes];
        let mut class_pos = std::collections::HashMap::new();
        for (idx, &cls) in classes.iter().enumerate() {
            class_pos.insert(cls, idx);
        }
        for (i, &label) in y.iter().enumerate() {
            if let Some(&idx) = class_pos.get(&label) {
                class_indices[idx].push(i);
            }
        }
        for (idx, indices) in class_indices.iter().enumerate() {
            if indices.is_empty() {
                return Err(FerroError::InsufficientSamples {
                    required: 1,
                    actual: 0,
                    context: format!("LDA: class {} has no samples", classes[idx]),
                });
            }
            let cnt_f = usize_to_f::<F>(indices.len())?;
            for &i in indices {
                for j in 0..n_features {
                    means[[idx, j]] += x[[i, j]];
                }
            }
            for j in 0..n_features {
                means[[idx, j]] /= cnt_f;
            }
        }

        // --- priors_  (sklearn :601-605) --------------------------------------
        // `priors=None` (default) ⇒ empirical `n_k / n` inferred from the data
        // (`:601-603`). `Some(p)` ⇒ `p` used VERBATIM (`:605`,
        // `self.priors_ = xp.asarray(self.priors)`). sklearn would mis-index a
        // wrong-length array, so reject it up front (R-DEV-4 length check).
        let priors = match &self.priors {
            None => {
                let mut priors = Array1::<F>::zeros(n_classes);
                for idx in 0..n_classes {
                    priors[idx] = usize_to_f::<F>(class_indices[idx].len())? / n_f;
                }
                priors
            }
            Some(p) => {
                if p.len() != n_classes {
                    return Err(FerroError::ShapeMismatch {
                        expected: vec![n_classes],
                        actual: vec![p.len()],
                        context: "LDA: priors length must match number of classes".into(),
                    });
                }
                let mut p = p.clone();
                // sklearn rejects negative priors (:607-608,
                // `if xp.any(self.priors_ < 0): raise ValueError("priors must
                // be non-negative")`).
                if p.iter().any(|&v| v < <F as num_traits::Zero>::zero()) {
                    return Err(FerroError::InvalidParameter {
                        name: "priors".into(),
                        reason: "priors must be non-negative".into(),
                    });
                }
                // sklearn renormalizes (with a UserWarning) when the priors do
                // not sum to 1 (:610-612, `if xp.abs(xp.sum(self.priors_) - 1.0)
                // > 1e-5: warnings.warn(...); self.priors_ = self.priors_ /
                // self.priors_.sum()`). FerroError has no warning channel; the
                // crate emits warnings via `eprintln!` (cf. qda.rs collinearity
                // warning, `discriminant_analysis.py:947`). The observable
                // contract is the renormalized `priors_`.
                let s = p.sum();
                let tol_sum = F::from(1e-5).unwrap_or_else(F::epsilon);
                if (s - <F as num_traits::One>::one()).abs() > tol_sum {
                    eprintln!("The priors do not sum to 1. Renormalizing");
                    for v in p.iter_mut() {
                        *v /= s;
                    }
                }
                p
            }
        };

        // --- solver dispatch (sklearn :627-650) -------------------------------
        // `lsqr` and `eigen` need only means_/priors_/class_indices; resolve
        // them above this point, then branch. `svd` falls through to the
        // existing two-SVD path below (BYTE-IDENTICAL — `shrinkage` must be
        // `None` for svd, sklearn `NotImplementedError` `:628-629`).
        match self.solver {
            Solver::Lsqr => {
                return self.solve_lstsq(
                    x,
                    &classes,
                    &class_indices,
                    &means,
                    &priors,
                    user_max,
                    n_features,
                );
            }
            Solver::Eigen => {
                // The generalized-eigenvalue solver (sklearn `_solve_eigen`,
                // `discriminant_analysis.py:421-485`): generalized `eigh(Sb, Sw)`
                // reduced to a standard symmetric eigenproblem via the Cholesky
                // factor of `Sw`. Supports `shrinkage` (like lsqr). See
                // [`LDA::solve_eigen`].
                return self.solve_eigen(
                    x,
                    &classes,
                    &class_indices,
                    &means,
                    &priors,
                    user_max,
                    n_features,
                );
            }
            Solver::Svd => {
                // sklearn: svd + shrinkage != None → NotImplementedError
                // ("shrinkage not supported with 'svd' solver.", `:628-629`).
                if !matches!(self.shrinkage, Shrinkage::None) {
                    return Err(FerroError::InvalidParameter {
                        name: "shrinkage".into(),
                        reason: "shrinkage not supported with svd solver".into(),
                    });
                }
            }
        }

        // --- xbar_ = priors_ @ means_  (sklearn :517) -------------------------
        let mut xbar = Array1::<F>::zeros(n_features);
        for j in 0..n_features {
            let mut acc = <F as num_traits::Zero>::zero();
            for idx in 0..n_classes {
                acc += priors[idx] * means[[idx, j]];
            }
            xbar[j] = acc;
        }

        // --- covariance_  (sklearn :509-510, `_class_cov` :128-172) -----------
        // When `store_covariance` is set, compute the shared within-class
        // covariance `Σ_k priors_[k] · cov(X_k)`, where `cov(X_k)` is the
        // MAXIMUM-LIKELIHOOD empirical covariance of class k's samples —
        // `empirical_covariance` calls `np.cov(Xg.T, bias=1)`
        // (`covariance/_empirical_covariance.py:109`), i.e. centered on the
        // class mean and normalized by `n_k` (NOT `n_k - 1`). Verified against
        // the live oracle: class-0 of the dispatch fixture yields the documented
        // `[[0.4296875, …], …]` only under the `bias=1` (÷n_k) normalization.
        let covariance = if self.store_covariance {
            let mut cov = Array2::<F>::zeros((n_features, n_features));
            for (idx, indices) in class_indices.iter().enumerate() {
                let nk = usize_to_f::<F>(indices.len())?;
                let prior_k = priors[idx];
                // cov(X_k)[a, b] = (1/n_k) Σ_i (x_ia - μ_ka)(x_ib - μ_kb).
                for a in 0..n_features {
                    for b in 0..n_features {
                        let mut acc = <F as num_traits::Zero>::zero();
                        for &i in indices {
                            acc += (x[[i, a]] - means[[idx, a]]) * (x[[i, b]] - means[[idx, b]]);
                        }
                        cov[[a, b]] += prior_k * (acc / nk);
                    }
                }
            }
            Some(cov)
        } else {
            None
        };

        // --- Xc = each sample minus its class mean (stacked; sklearn :512-519) -
        let mut xc = Array2::<F>::zeros((n_samples, n_features));
        for (idx, indices) in class_indices.iter().enumerate() {
            for &i in indices {
                for j in 0..n_features {
                    xc[[i, j]] = x[[i, j]] - means[[idx, j]];
                }
            }
        }

        // --- std = population std of Xc per column (ddof=0; sklearn :522-524) --
        // numpy std: sqrt(mean((Xc - mean(Xc))^2)). Xc columns already have ~0
        // mean by construction, but follow numpy exactly (subtract the column
        // mean) for ULP fidelity.
        let mut std = Array1::<F>::zeros(n_features);
        for j in 0..n_features {
            let mut col_mean = <F as num_traits::Zero>::zero();
            for i in 0..n_samples {
                col_mean += xc[[i, j]];
            }
            col_mean /= n_f;
            let mut var = <F as num_traits::Zero>::zero();
            for i in 0..n_samples {
                let d = xc[[i, j]] - col_mean;
                var += d * d;
            }
            var /= n_f;
            let s = var.sqrt();
            std[j] = if s == <F as num_traits::Zero>::zero() {
                <F as num_traits::One>::one()
            } else {
                s
            };
        }

        // --- Xw = sqrt(1/(n-c)) * (Xc / std)  (sklearn :525-528) --------------
        let denom = usize_to_f::<F>(n_samples - n_classes)?;
        let fac_sqrt = (<F as num_traits::One>::one() / denom).sqrt();
        let mut xw = Array2::<F>::zeros((n_samples, n_features));
        for i in 0..n_samples {
            for j in 0..n_features {
                xw[[i, j]] = fac_sqrt * (xc[[i, j]] / std[j]);
            }
        }

        // --- first SVD: within whitening (sklearn :530-534) -------------------
        // sklearn's svd-solver rank threshold `tol` (constructor default `1e-4`,
        // `discriminant_analysis.py:354,362`), now configurable via
        // `LDA::with_tol`. Default `1e-4` ⇒ byte-identical to the prior hardcode.
        let tol = self.tol;
        let (s1, vt1) = svd_s_vt::<F>(&xw)?;
        let rank1 = s1.iter().filter(|&&v| v > tol).count();
        if rank1 == 0 {
            return Err(FerroError::NumericalInstability {
                message: "LDA: within-class scatter has rank 0 (all features constant)".into(),
            });
        }
        // scalings = (Vt[:rank]/std).T / S[:rank]   -> (n_features, rank1)
        let mut scalings1 = Array2::<F>::zeros((n_features, rank1));
        for k in 0..rank1 {
            let sk = s1[k];
            for j in 0..n_features {
                scalings1[[j, k]] = (vt1[[k, j]] / std[j]) / sk;
            }
        }

        // --- between-class scaled centers (sklearn :535-541) ------------------
        // Xb[i] = sqrt(n * priors_[i] * fac2) * (means_[i] - xbar_)  then @ scalings.
        let fac2 = if n_classes == 1 {
            <F as num_traits::One>::one()
        } else {
            <F as num_traits::One>::one() / usize_to_f::<F>(n_classes - 1)?
        };
        let mut xb_centers = Array2::<F>::zeros((n_classes, n_features));
        for idx in 0..n_classes {
            let w = (n_f * priors[idx] * fac2).sqrt();
            for j in 0..n_features {
                xb_centers[[idx, j]] = w * (means[[idx, j]] - xbar[j]);
            }
        }
        let xb = xb_centers.dot(&scalings1); // (n_classes, rank1)

        // --- second SVD: between-class projection (sklearn :545-555) ----------
        let (s2, vt2) = svd_s_vt::<F>(&xb)?;

        // explained_variance_ratio_ = (S2^2 / sum(S2^2))[:max_components] (:550-552)
        let mut sum_sq = <F as num_traits::Zero>::zero();
        for &v in s2.iter() {
            sum_sq += v * v;
        }
        let evr_len = user_max.min(s2.len());
        let mut explained_variance_ratio = Array1::<F>::zeros(evr_len);
        for k in 0..evr_len {
            explained_variance_ratio[k] = if sum_sq > <F as num_traits::Zero>::zero() {
                (s2[k] * s2[k]) / sum_sq
            } else {
                <F as num_traits::Zero>::zero()
            };
        }

        // rank2 = sum(S2 > tol * S2[0])  (sklearn :554)
        let s2_0 = if s2.is_empty() {
            <F as num_traits::Zero>::zero()
        } else {
            s2[0]
        };
        let rank2 = s2.iter().filter(|&&v| v > tol * s2_0).count();
        if rank2 == 0 {
            return Err(FerroError::NumericalInstability {
                message: "LDA: between-class scatter has rank 0 (classes coincide)".into(),
            });
        }

        // scalings_ = scalings @ Vt2.T[:, :rank2]   -> (n_features, rank2)
        // Vt2 is (k2, rank1); Vt2.T is (rank1, k2); take first rank2 columns.
        let mut scalings = Array2::<F>::zeros((n_features, rank2));
        for j in 0..n_features {
            for c in 0..rank2 {
                let mut acc = <F as num_traits::Zero>::zero();
                for k in 0..rank1 {
                    // Vt2.T[k, c] = Vt2[c, k]
                    acc += scalings1[[j, k]] * vt2[[c, k]];
                }
                scalings[[j, c]] = acc;
            }
        }

        // --- coef_ / intercept_  (sklearn :556-559) ---------------------------
        // coef = (means_ - xbar_) @ scalings_     (n_classes, rank2)
        let mut centered_means = Array2::<F>::zeros((n_classes, n_features));
        for idx in 0..n_classes {
            for j in 0..n_features {
                centered_means[[idx, j]] = means[[idx, j]] - xbar[j];
            }
        }
        let coef_lowrank = centered_means.dot(&scalings); // (n_classes, rank2)

        // intercept_ = -0.5 * sum(coef^2, axis=1) + log(priors_)
        let neg_half = -half::<F>();
        let mut intercept = Array1::<F>::zeros(n_classes);
        for idx in 0..n_classes {
            let mut sq = <F as num_traits::Zero>::zero();
            for c in 0..rank2 {
                sq += coef_lowrank[[idx, c]] * coef_lowrank[[idx, c]];
            }
            intercept[idx] = neg_half * sq + priors[idx].ln();
        }

        // coef_ = coef @ scalings_.T              (n_classes, n_features)
        let coef = coef_lowrank.dot(&scalings.t());

        // intercept_ -= xbar_ @ coef_.T           (subtract per class)
        for idx in 0..n_classes {
            let mut dot = <F as num_traits::Zero>::zero();
            for j in 0..n_features {
                dot += xbar[j] * coef[[idx, j]];
            }
            intercept[idx] -= dot;
        }

        Ok(FittedLDA {
            scalings,
            means,
            xbar,
            priors,
            coef,
            intercept,
            explained_variance_ratio,
            covariance,
            classes,
            max_components: user_max,
            n_features,
        })
    }
}

// ---------------------------------------------------------------------------
// Fit (sklearn _solve_lstsq) — the lsqr solver
// ---------------------------------------------------------------------------

impl<F: LinalgFloat + ScalarOperand> LDA<F> {
    /// The least-squares solver (sklearn's `_solve_lstsq`,
    /// `discriminant_analysis.py:365-419`), dispatched from [`Fit::fit`] when
    /// [`Solver::Lsqr`] is selected.
    ///
    /// Computes (`:412-418`):
    /// - `covariance_ = Σ_k priors_[k] · cov(X_k)` where `cov(X_k)` applies the
    ///   configured [`Shrinkage`] to class `k`'s empirical covariance
    ///   (`_class_cov` `:128-172`, `_cov` `:36-93`).
    /// - `coef_ = lstsq(covariance_, means_.T)[0].T` (`:416`), via
    ///   [`ferray::linalg::lstsq`] (multi-RHS).
    /// - `intercept_ = -½·diag(means_ @ coef_.T) + log(priors_)` (`:417-418`).
    ///
    /// Unlike sklearn's `svd` solver, the lsqr `coef_` is the FULL-space
    /// discriminant `(n_classes, n_features)` — NO `scalings_`/`xbar_`/
    /// `explained_variance_ratio_` and NO `transform` (sklearn raises
    /// `NotImplementedError` for `transform` under lsqr, `:676-679`); here
    /// [`Transform`] returns an error because `scalings_` is the zero matrix /
    /// `xbar_` is zero, and `transform` slices to `max_components` of a meaningless
    /// projection — we instead document that `transform` is unsupported by
    /// recording `max_components = 0` so the projection is empty (mirroring
    /// sklearn's "dimensionality reduction is not supported" for lsqr).
    ///
    /// `covariance_` is ALWAYS populated for lsqr (sklearn `:413`, the attribute
    /// is set regardless of `store_covariance`), exposed via
    /// [`FittedLDA::covariance`].
    ///
    /// `decision_function`/`predict`/`predict_proba` work identically to the svd
    /// path because they consume `coef_`/`intercept_` only.
    ///
    /// NOTE: the binary-collapse of `coef_`/`intercept_` to a single row
    /// (sklearn `:651-657`) is NOT applied here, matching the existing svd path
    /// (open prereq blocker #600); `coef_` stays `(n_classes, n_features)`.
    ///
    /// # Errors
    ///
    /// - [`FerroError::InvalidParameter`] if [`Shrinkage::Fixed`]`(s)` has
    ///   `s ∉ [0, 1]` (sklearn `Interval(Real, 0, 1, closed="both")`, `:339`).
    /// - [`FerroError::NumericalInstability`] if the least-squares solve fails.
    #[allow(
        clippy::too_many_arguments,
        reason = "the lsqr solver consumes the same pre-resolved fit state (classes/indices/means/priors/dims) the svd path computes; threading them avoids recomputation"
    )]
    fn solve_lstsq(
        &self,
        x: &Array2<F>,
        classes: &[usize],
        class_indices: &[Vec<usize>],
        means: &Array2<F>,
        priors: &Array1<F>,
        _user_max: usize,
        n_features: usize,
    ) -> Result<FittedLDA<F>, FerroError> {
        let n_classes = classes.len();

        // Validate Fixed(s): sklearn Interval(Real, 0, 1, closed="both") (:339).
        if let Shrinkage::Fixed(s) = self.shrinkage
            && (s < <F as num_traits::Zero>::zero() || s > <F as num_traits::One>::one())
        {
            return Err(FerroError::InvalidParameter {
                name: "shrinkage".into(),
                reason: "shrinkage float must be in [0, 1]".into(),
            });
        }

        // covariance_ = Σ_k priors_[k] · cov(X_k)  (sklearn _class_cov :167-172).
        let mut covariance = Array2::<F>::zeros((n_features, n_features));
        for (idx, indices) in class_indices.iter().enumerate() {
            // Gather class-k rows.
            let nk = indices.len();
            let mut xg = Array2::<F>::zeros((nk, n_features));
            for (r, &i) in indices.iter().enumerate() {
                for j in 0..n_features {
                    xg[[r, j]] = x[[i, j]];
                }
            }
            let cov_k = cov_shrunk(&xg, self.shrinkage)?;
            let prior_k = priors[idx];
            for a in 0..n_features {
                for b in 0..n_features {
                    covariance[[a, b]] += prior_k * cov_k[[a, b]];
                }
            }
        }

        // coef_ = lstsq(covariance_, means_.T)[0].T  (sklearn :416).
        // Solve `covariance_ @ X = means_.T` for X of shape (n_features,
        // n_classes); X = lstsq(...)[0]; coef_ = X.T = (n_classes, n_features).
        let means_t = means.t().to_owned(); // (n_features, n_classes)
        let sol = lstsq_multi(&covariance, &means_t)?; // (n_features, n_classes)
        let coef = sol.t().to_owned(); // (n_classes, n_features)

        // intercept_ = -0.5 * diag(means_ @ coef_.T) + log(priors_)  (:417-418).
        // diag(means_ @ coef_.T)[k] = Σ_j means_[k,j] · coef_[k,j].
        let neg_half = -half::<F>();
        let mut intercept = Array1::<F>::zeros(n_classes);
        for k in 0..n_classes {
            let mut dot = <F as num_traits::Zero>::zero();
            for j in 0..n_features {
                dot += means[[k, j]] * coef[[k, j]];
            }
            intercept[k] = neg_half * dot + priors[k].ln();
        }

        // lsqr does NOT support dimensionality reduction (sklearn :372-373,
        // :676-679): no scalings_/xbar_/explained_variance_ratio_. Set them to
        // empty/zero and `max_components = 0` so `transform` yields a `(n, 0)`
        // projection (the lsqr "no transform" contract).
        Ok(FittedLDA {
            scalings: Array2::<F>::zeros((n_features, 0)),
            means: means.to_owned(),
            xbar: Array1::<F>::zeros(n_features),
            priors: priors.to_owned(),
            coef,
            intercept,
            explained_variance_ratio: Array1::<F>::zeros(0),
            covariance: Some(covariance),
            classes: classes.to_vec(),
            max_components: 0,
            n_features,
        })
    }
}

// ---------------------------------------------------------------------------
// Fit (sklearn _solve_eigen) — the eigen solver
// ---------------------------------------------------------------------------

impl<F: LinalgFloat + ScalarOperand> LDA<F> {
    /// The generalized-eigenvalue solver (sklearn's `_solve_eigen`,
    /// `discriminant_analysis.py:421-485`), dispatched from [`Fit::fit`] when
    /// [`Solver::Eigen`] is selected.
    ///
    /// Computes (`:466-485`):
    /// - `Sw = Σ_k priors_[k] · cov(X_k)` — the within-class scatter (the
    ///   shrinkage-aware `_class_cov`, `:467-471`), stored as `covariance_`.
    /// - `St = cov(WHOLE X)` — the total scatter of all of `X` (the SAME `_cov`
    ///   with the configured shrinkage applied to the full centered `X`, `:472`).
    /// - `Sb = St - Sw` — the between-class scatter (`:473`).
    /// - the GENERALIZED symmetric-definite eigenproblem `eigh(Sb, Sw)` (`:475`),
    ///   sorted by DESCENDING eigenvalue (`:479`,
    ///   `evecs = evecs[:, argsort(evals)[::-1]]`).
    /// - `explained_variance_ratio_ = sort(evals / Σ evals)[::-1][:max_components]`
    ///   (`:476-478`).
    /// - `scalings_ = evecs` (`:481`), `coef_ = (means_ @ evecs) @ evecs.T`
    ///   (`:482`), `intercept_ = -½·diag(means_ @ coef_.T) + log(priors_)`
    ///   (`:483-485`).
    ///
    /// # The Cholesky reduction
    ///
    /// ferray exposes the STANDARD symmetric eigensolver
    /// ([`ferray::linalg::eigh`]) and [`ferray::linalg::cholesky`], not a
    /// generalized solver, so the generalized problem `Sb·v = λ·Sw·v` (with `Sw`
    /// SPD) is reduced to standard form: let `Sw = L·Lᵀ` (Cholesky); then
    /// `M = L⁻¹·Sb·L⁻ᵀ` is symmetric and `eigh(M)` gives the same eigenvalues
    /// `λ`, with generalized eigenvectors `v = L⁻ᵀ·w` (`w` the standard
    /// eigenvectors of `M`). `M` is symmetrized (`M = (M + Mᵀ)/2`) to kill
    /// rounding asymmetry before [`eigh`]. Because `coef_ = (means_@evecs)@evecsᵀ`
    /// is invariant to the per-column SIGN and ORDER of `evecs`, `coef_`/
    /// `intercept_` (hence `predict`/`predict_proba`/`decision_function`) match
    /// sklearn exactly regardless of the eigenvector sign/order ambiguity. The
    /// explained-variance ratio is sorted by eigenvalue, so it is order-stable.
    ///
    /// `scalings_` is `(n_features, n_features)` (sklearn keeps ALL columns for
    /// the eigen solver, `:481`); `transform` slices to `[:, :max_components]`.
    /// Eigen has NO `xbar_` (sklearn `_solve_eigen` does not set it, `:466-485`),
    /// so `transform` is the un-centered `X @ scalings_` (`:687`).
    ///
    /// # Errors
    ///
    /// - [`FerroError::InvalidParameter`] if [`Shrinkage::Fixed`]`(s)` has
    ///   `s ∉ [0, 1]` (sklearn `Interval(Real, 0, 1, closed="both")`, `:339`).
    /// - [`FerroError::NumericalInstability`] if the Cholesky factorization,
    ///   inversion, or eigendecomposition fails.
    #[allow(
        clippy::too_many_arguments,
        reason = "the eigen solver consumes the same pre-resolved fit state (X/classes/indices/means/priors/dims) the svd path computes; threading them avoids recomputation"
    )]
    #[allow(
        clippy::needless_range_loop,
        reason = "explicit index loops mirror sklearn's matrix arithmetic per-row/per-column"
    )]
    fn solve_eigen(
        &self,
        x: &Array2<F>,
        classes: &[usize],
        class_indices: &[Vec<usize>],
        means: &Array2<F>,
        priors: &Array1<F>,
        max_components: usize,
        n_features: usize,
    ) -> Result<FittedLDA<F>, FerroError> {
        let n_classes = classes.len();

        // Validate Fixed(s): sklearn Interval(Real, 0, 1, closed="both") (:339).
        if let Shrinkage::Fixed(s) = self.shrinkage
            && (s < <F as num_traits::Zero>::zero() || s > <F as num_traits::One>::one())
        {
            return Err(FerroError::InvalidParameter {
                name: "shrinkage".into(),
                reason: "shrinkage float must be in [0, 1]".into(),
            });
        }

        // Sw = Σ_k priors_[k] · cov(X_k)  — within-class scatter (covariance_,
        // sklearn :467-471, `_class_cov` :128-172, `_cov` :36-93 shrinkage-aware).
        let mut sw = Array2::<F>::zeros((n_features, n_features));
        for (idx, indices) in class_indices.iter().enumerate() {
            let nk = indices.len();
            let mut xg = Array2::<F>::zeros((nk, n_features));
            for (r, &i) in indices.iter().enumerate() {
                for j in 0..n_features {
                    xg[[r, j]] = x[[i, j]];
                }
            }
            let cov_k = cov_shrunk(&xg, self.shrinkage)?;
            let prior_k = priors[idx];
            for a in 0..n_features {
                for b in 0..n_features {
                    sw[[a, b]] += prior_k * cov_k[[a, b]];
                }
            }
        }

        // St = _cov(WHOLE X, shrinkage)  — total scatter of all of X (sklearn
        // :472). Same `_cov` (shrinkage applied to the FULL centered X).
        let st = cov_shrunk(&x.to_owned(), self.shrinkage)?;

        // Sb = St - Sw  — between-class scatter (sklearn :473).
        let mut sb = Array2::<F>::zeros((n_features, n_features));
        for a in 0..n_features {
            for b in 0..n_features {
                sb[[a, b]] = st[[a, b]] - sw[[a, b]];
            }
        }

        // --- generalized eigh(Sb, Sw) via Cholesky reduction (sklearn :475) ---
        // Sw = L·Lᵀ; M = L⁻¹·Sb·L⁻ᵀ (symmetric); (evals, W) = eigh(M);
        // generalized eigenvectors evecs = L⁻ᵀ·W.
        let l = cholesky_lower(&sw)?; // lower-triangular, Sw = L·Lᵀ
        let l_inv = matrix_inverse(&l)?; // L⁻¹  (n_features, n_features)
        // M = L⁻¹ · Sb · L⁻ᵀ
        let l_inv_t = l_inv.t().to_owned();
        let m = l_inv.dot(&sb).dot(&l_inv_t);
        // Symmetrize M = (M + Mᵀ)/2 to kill rounding asymmetry before eigh.
        let mut m_sym = Array2::<F>::zeros((n_features, n_features));
        let half_f = half::<F>();
        for a in 0..n_features {
            for b in 0..n_features {
                m_sym[[a, b]] = (m[[a, b]] + m[[b, a]]) * half_f;
            }
        }
        // eigh(M) — ascending eigenvalues, eigenvectors as columns.
        let (evals_asc, w) = eigh_sym(&m_sym)?;
        // Generalized eigenvectors: evecs = L⁻ᵀ · W  (n_features, n_features).
        let evecs_asc = l_inv_t.dot(&w);

        // --- sort by DESCENDING eigenvalue (sklearn :479) ---------------------
        // argsort(evals)[::-1]: indices ordered by descending eigenvalue.
        let n = evals_asc.len();
        let mut order: Vec<usize> = (0..n).collect();
        // evals_asc is ascending; reverse it for descending order.
        order.reverse();
        let mut evals_desc = Array1::<F>::zeros(n);
        let mut evecs = Array2::<F>::zeros((n_features, n));
        for (new_c, &old_c) in order.iter().enumerate() {
            evals_desc[new_c] = evals_asc[old_c];
            for j in 0..n_features {
                evecs[[j, new_c]] = evecs_asc[[j, old_c]];
            }
        }

        // --- explained_variance_ratio_  (sklearn :476-478) --------------------
        // sort(evals / sum(evals))[::-1][:max_components]. The ratio is over ALL
        // eigenvalues, sorted descending, then truncated to max_components.
        let mut sum_evals = <F as num_traits::Zero>::zero();
        for &v in evals_asc.iter() {
            sum_evals += v;
        }
        let evr_len = max_components.min(n);
        let mut explained_variance_ratio = Array1::<F>::zeros(evr_len);
        for k in 0..evr_len {
            // evals_desc is already the descending sort of evals; dividing by the
            // (sign-stable) sum preserves the sort order sklearn applies.
            explained_variance_ratio[k] = if sum_evals != <F as num_traits::Zero>::zero() {
                evals_desc[k] / sum_evals
            } else {
                <F as num_traits::Zero>::zero()
            };
        }

        // --- scalings_ / coef_ / intercept_  (sklearn :481-485) ---------------
        // scalings_ = evecs (ALL columns; sklearn keeps the full (n_features,
        // n_features) for the eigen solver, :481).
        let scalings = evecs.clone();
        // coef_ = (means_ @ evecs) @ evecs.T   (n_classes, n_features).
        // This is invariant to the per-column sign/order of evecs, so it matches
        // sklearn exactly despite the eigenvector sign/order ambiguity.
        let coef = means.dot(&evecs).dot(&evecs.t());
        // intercept_ = -0.5 * diag(means_ @ coef_.T) + log(priors_)  (:483-485).
        // diag(means_ @ coef_.T)[k] = Σ_j means_[k,j] · coef_[k,j].
        let neg_half = -half::<F>();
        let mut intercept = Array1::<F>::zeros(n_classes);
        for k in 0..n_classes {
            let mut dot = <F as num_traits::Zero>::zero();
            for j in 0..n_features {
                dot += means[[k, j]] * coef[[k, j]];
            }
            intercept[k] = neg_half * dot + priors[k].ln();
        }

        Ok(FittedLDA {
            scalings,
            means: means.to_owned(),
            // Eigen has NO xbar_ (sklearn `_solve_eigen` never sets it,
            // :466-485); `transform` is the un-centered `X @ scalings_` (:687).
            xbar: Array1::<F>::zeros(n_features),
            priors: priors.to_owned(),
            coef,
            intercept,
            explained_variance_ratio,
            // covariance_ = Sw, ALWAYS populated for eigen (sklearn :467-469,
            // the attribute is set regardless of store_covariance).
            covariance: Some(sw),
            classes: classes.to_vec(),
            max_components,
            n_features,
        })
    }
}

// ---------------------------------------------------------------------------
// Transform (sklearn svd transform)
// ---------------------------------------------------------------------------

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

    /// Project `x` onto the discriminant axes: `((X - xbar_) @ scalings_)[:, :n]`
    /// where `n = max_components`. Mirrors sklearn's svd-solver `transform`
    /// (`discriminant_analysis.py:684-685,689`).
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if `x.ncols()` does not match the
    /// number of features seen during fitting.
    fn transform(&self, x: &Array2<F>) -> Result<Array2<F>, FerroError> {
        if x.ncols() != self.n_features {
            return Err(FerroError::ShapeMismatch {
                expected: vec![x.nrows(), self.n_features],
                actual: vec![x.nrows(), x.ncols()],
                context: "FittedLDA::transform".into(),
            });
        }
        // (X - xbar_) @ scalings_
        let mut xc = x.to_owned();
        for mut row in xc.rows_mut() {
            for j in 0..self.n_features {
                row[j] = row[j] - self.xbar[j];
            }
        }
        let projected = xc.dot(&self.scalings);
        // Slice to [:, :max_components] (sklearn :689).
        let keep = self.max_components.min(projected.ncols());
        Ok(projected.slice(ndarray::s![.., ..keep]).to_owned())
    }
}

// ---------------------------------------------------------------------------
// Predict (argmax of the affine decision_function)
// ---------------------------------------------------------------------------

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

    /// Classify samples by argmax of the affine `decision_function`
    /// (`classes_[argmax(X @ coef_.T + intercept_)]`), mirroring sklearn's
    /// `predict` (the `LinearClassifierMixin`, `discriminant_analysis.py:739`).
    /// The argmax follows numpy's first-max-wins tie-breaking.
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if the number of features does not
    /// match the model.
    fn predict(&self, x: &Array2<F>) -> Result<Array1<usize>, FerroError> {
        let decision = self.decision_function(x)?;
        let n_samples = decision.nrows();
        let n_classes = decision.ncols();
        let mut predictions = Array1::<usize>::zeros(n_samples);
        for i in 0..n_samples {
            let mut best_idx = 0usize;
            let mut best = decision[[i, 0]];
            for c in 1..n_classes {
                let v = decision[[i, c]];
                // numpy argmax: strictly-greater wins; ties keep first index.
                if v > best {
                    best = v;
                    best_idx = c;
                }
            }
            predictions[i] = self.classes[best_idx];
        }
        Ok(predictions)
    }
}

// ---------------------------------------------------------------------------
// Introspection
// ---------------------------------------------------------------------------

impl<F: Float + Send + Sync + 'static> HasClasses for FittedLDA<F> {
    fn classes(&self) -> &[usize] {
        &self.classes
    }

    fn n_classes(&self) -> usize {
        self.classes.len()
    }
}

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

impl<F: LinalgFloat + ScalarOperand> PipelineEstimator<F> for LDA<F> {
    /// Fit LDA using the pipeline interface.
    ///
    /// # Errors
    ///
    /// Propagates errors from [`Fit::fit`].
    fn fit_pipeline(
        &self,
        x: &Array2<F>,
        y: &Array1<F>,
    ) -> Result<Box<dyn FittedPipelineEstimator<F>>, FerroError> {
        let y_usize: Array1<usize> = y.mapv(|v| v.to_usize().unwrap_or(0));
        let fitted = self.fit(x, &y_usize)?;
        Ok(Box::new(FittedLDAPipeline(fitted)))
    }
}

/// Wrapper for pipeline integration that converts predictions to float.
struct FittedLDAPipeline<F>(FittedLDA<F>);

impl<F: Float + Send + Sync + 'static> FittedPipelineEstimator<F> for FittedLDAPipeline<F> {
    /// Predict via the pipeline interface, returning float class labels.
    fn predict_pipeline(&self, x: &Array2<F>) -> Result<Array1<F>, FerroError> {
        let preds = self.0.predict(x)?;
        Ok(preds.mapv(|v| NumCast::from(v).unwrap_or_else(F::nan)))
    }
}

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

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

    // ------------------------------------------------------------------
    // Helpers
    // ------------------------------------------------------------------

    fn linearly_separable_2d() -> (Array2<f64>, Array1<usize>) {
        // Two well-separated Gaussian clusters.
        let x = Array2::from_shape_vec(
            (8, 2),
            vec![
                1.0, 1.0, 1.5, 1.2, 0.8, 0.9, 1.1, 1.3, // class 0
                6.0, 6.0, 6.2, 5.8, 5.9, 6.1, 6.3, 5.7, // class 1
            ],
        )
        .unwrap();
        let y = array![0, 0, 0, 0, 1, 1, 1, 1];
        (x, y)
    }

    fn three_class_data() -> (Array2<f64>, Array1<usize>) {
        let x = Array2::from_shape_vec(
            (9, 2),
            vec![
                0.0, 0.0, 0.5, 0.1, 0.1, 0.5, // class 0
                5.0, 0.0, 5.2, 0.3, 4.8, 0.1, // class 1
                0.0, 5.0, 0.1, 5.2, 0.3, 4.8, // class 2
            ],
        )
        .unwrap();
        let y = array![0, 0, 0, 1, 1, 1, 2, 2, 2];
        (x, y)
    }

    // ------------------------------------------------------------------

    #[test]
    fn test_lda_fit_returns_fitted() {
        let (x, y) = linearly_separable_2d();
        let lda = LDA::<f64>::new(Some(1));
        let fitted = lda.fit(&x, &y).unwrap();
        // scalings_ is (n_features, rank2); for a binary 2-feature set rank2 = 1.
        assert_eq!(fitted.scalings().ncols(), 1);
        assert_eq!(fitted.scalings().nrows(), 2);
    }

    #[test]
    fn test_lda_default_n_components() {
        // With 2 classes the default n_components = min(1, n_features) = 1.
        let (x, y) = linearly_separable_2d();
        let lda = LDA::<f64>::default();
        let fitted = lda.fit(&x, &y).unwrap();
        // transform output is truncated to max_components = 1.
        assert_eq!(fitted.transform(&x).unwrap().ncols(), 1);
    }

    #[test]
    fn test_lda_transform_shape() {
        let (x, y) = linearly_separable_2d();
        let lda = LDA::<f64>::new(Some(1));
        let fitted = lda.fit(&x, &y).unwrap();
        let proj = fitted.transform(&x).unwrap();
        assert_eq!(proj.dim(), (8, 1));
    }

    #[test]
    fn test_lda_predict_accuracy_binary() {
        let (x, y) = linearly_separable_2d();
        let lda = LDA::<f64>::new(Some(1));
        let fitted = lda.fit(&x, &y).unwrap();
        let preds = fitted.predict(&x).unwrap();
        let correct = preds.iter().zip(y.iter()).filter(|(p, a)| *p == *a).count();
        assert_eq!(correct, 8, "All 8 samples should be classified correctly");
    }

    #[test]
    fn test_lda_predict_three_classes() {
        let (x, y) = three_class_data();
        let lda = LDA::<f64>::new(Some(2));
        let fitted = lda.fit(&x, &y).unwrap();
        let preds = fitted.predict(&x).unwrap();
        let correct = preds.iter().zip(y.iter()).filter(|(p, a)| *p == *a).count();
        assert!(correct >= 7, "Expected at least 7/9 correct, got {correct}");
    }

    #[test]
    fn test_lda_explained_variance_ratio_positive() {
        let (x, y) = linearly_separable_2d();
        let lda = LDA::<f64>::new(Some(1));
        let fitted = lda.fit(&x, &y).unwrap();
        for &v in fitted.explained_variance_ratio() {
            assert!(v >= 0.0);
        }
    }

    #[test]
    fn test_lda_explained_variance_ratio_le_1() {
        let (x, y) = three_class_data();
        let lda = LDA::<f64>::new(Some(2));
        let fitted = lda.fit(&x, &y).unwrap();
        let total: f64 = fitted.explained_variance_ratio().iter().sum();
        assert!(total <= 1.0 + 1e-9, "total={total}");
    }

    /// R-CHAR-3 oracle pin for `explained_variance_ratio_` (REQ-13). Expected
    /// values are the live sklearn 1.5.2
    /// `LinearDiscriminantAnalysis().fit(X,y).explained_variance_ratio_` on the
    /// 3-class / 2-feature balanced set (same data as `divergence_lda_fit.rs`):
    /// ```text
    /// python3 -c "import numpy as np; \
    ///   from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as L; \
    ///   X=np.array([[0.,0.],[1.,.5],[.5,1.],[1.,1.],[4.,4.],[5.,4.5],[4.5,5.],[5.,5.],\
    ///               [0.,5.],[1.,6.],[.5,5.5],[1.,5.]]); \
    ///   y=np.array([0,0,0,0,1,1,1,1,2,2,2,2]); \
    ///   print(repr(L().fit(X,y).explained_variance_ratio_.tolist()))"
    /// # [0.6428683117561941, 0.3571316882438059]
    /// ```
    #[test]
    fn test_lda_explained_variance_ratio_oracle() {
        const SK_EVR: [f64; 2] = [0.6428683117561941, 0.3571316882438059];
        let x = Array2::from_shape_vec(
            (12, 2),
            vec![
                0.0, 0.0, 1.0, 0.5, 0.5, 1.0, 1.0, 1.0, 4.0, 4.0, 5.0, 4.5, 4.5, 5.0, 5.0, 5.0,
                0.0, 5.0, 1.0, 6.0, 0.5, 5.5, 1.0, 5.0,
            ],
        )
        .unwrap();
        let y = array![0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2];
        let fitted = LDA::<f64>::new(Some(2)).fit(&x, &y).unwrap();
        let evr = fitted.explained_variance_ratio();
        assert_eq!(evr.len(), 2);
        for k in 0..2 {
            assert_abs_diff_eq!(evr[k], SK_EVR[k], epsilon = 1e-9);
        }
    }

    /// R-CHAR-3 oracle pin for `coef_`/`intercept_`/`xbar_` (REQ-8). Live
    /// sklearn 1.5.2 attributes on the same 3-class / 2-feature set:
    /// ```text
    /// python3 -c "... ; m=L().fit(X,y); \
    ///   print(repr(m.coef_.tolist())); print(repr(m.intercept_.tolist())); \
    ///   print(repr(m.xbar_.tolist()))"
    /// # coef_ [[2.2582417582417564, -14.02747252747253],
    /// #        [13.335164835164827, -2.950549450549442],
    /// #        [-15.593406593406584, 16.978021978021978]]
    /// # intercept_ [25.208393205837393, -32.94545294800878, -56.65081009086592]
    /// # xbar_ [1.958333333333333, 3.541666666666666]
    /// ```
    #[test]
    fn test_lda_coef_intercept_xbar_oracle() {
        const SK_COEF: [[f64; 2]; 3] = [
            [2.2582417582417564, -14.02747252747253],
            [13.335164835164827, -2.950549450549442],
            [-15.593406593406584, 16.978021978021978],
        ];
        const SK_INTERCEPT: [f64; 3] = [25.208393205837393, -32.94545294800878, -56.65081009086592];
        const SK_XBAR: [f64; 2] = [1.958333333333333, 3.541666666666666];
        let x = Array2::from_shape_vec(
            (12, 2),
            vec![
                0.0, 0.0, 1.0, 0.5, 0.5, 1.0, 1.0, 1.0, 4.0, 4.0, 5.0, 4.5, 4.5, 5.0, 5.0, 5.0,
                0.0, 5.0, 1.0, 6.0, 0.5, 5.5, 1.0, 5.0,
            ],
        )
        .unwrap();
        let y = array![0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2];
        let fitted = LDA::<f64>::new(Some(2)).fit(&x, &y).unwrap();
        for i in 0..3 {
            for (j, &expected) in SK_COEF[i].iter().enumerate() {
                assert_abs_diff_eq!(fitted.coef()[[i, j]], expected, epsilon = 1e-9);
            }
            assert_abs_diff_eq!(fitted.intercept()[i], SK_INTERCEPT[i], epsilon = 1e-9);
        }
        for (j, &expected) in SK_XBAR.iter().enumerate() {
            assert_abs_diff_eq!(fitted.xbar()[j], expected, epsilon = 1e-12);
        }
    }

    #[test]
    fn test_lda_classes_accessor() {
        let (x, y) = linearly_separable_2d();
        let lda = LDA::<f64>::new(Some(1));
        let fitted = lda.fit(&x, &y).unwrap();
        assert_eq!(fitted.classes(), &[0usize, 1]);
    }

    #[test]
    fn test_lda_means_shape() {
        // means_ is now in the ORIGINAL feature space (n_classes, n_features).
        let (x, y) = three_class_data();
        let lda = LDA::<f64>::new(Some(2));
        let fitted = lda.fit(&x, &y).unwrap();
        assert_eq!(fitted.means().dim(), (3, 2));
    }

    #[test]
    fn test_lda_transform_shape_mismatch() {
        let (x, y) = linearly_separable_2d();
        let lda = LDA::<f64>::new(Some(1));
        let fitted = lda.fit(&x, &y).unwrap();
        let x_bad = Array2::<f64>::zeros((3, 5));
        assert!(fitted.transform(&x_bad).is_err());
    }

    #[test]
    fn test_lda_predict_shape_mismatch() {
        let (x, y) = linearly_separable_2d();
        let lda = LDA::<f64>::new(Some(1));
        let fitted = lda.fit(&x, &y).unwrap();
        let x_bad = Array2::<f64>::zeros((3, 5));
        assert!(fitted.predict(&x_bad).is_err());
    }

    #[test]
    fn test_lda_error_zero_n_components() {
        let (x, y) = linearly_separable_2d();
        let lda = LDA::<f64>::new(Some(0));
        assert!(lda.fit(&x, &y).is_err());
    }

    #[test]
    fn test_lda_error_n_components_too_large() {
        let (x, y) = linearly_separable_2d(); // 2 classes → max 1 component
        let lda = LDA::<f64>::new(Some(5));
        assert!(lda.fit(&x, &y).is_err());
    }

    #[test]
    fn test_lda_error_single_class() {
        let x =
            Array2::from_shape_vec((4, 2), vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]).unwrap();
        let y = array![0usize, 0, 0, 0];
        let lda = LDA::<f64>::new(None);
        assert!(lda.fit(&x, &y).is_err());
    }

    #[test]
    fn test_lda_error_shape_mismatch_fit() {
        let x = Array2::<f64>::zeros((4, 2));
        let y = array![0usize, 1]; // wrong length
        let lda = LDA::<f64>::new(None);
        assert!(lda.fit(&x, &y).is_err());
    }

    #[test]
    fn test_lda_error_insufficient_samples() {
        let x = Array2::<f64>::zeros((1, 2));
        let y = array![0usize];
        let lda = LDA::<f64>::new(None);
        assert!(lda.fit(&x, &y).is_err());
    }

    #[test]
    fn test_lda_scalings_accessor() {
        let (x, y) = linearly_separable_2d();
        let lda = LDA::<f64>::new(Some(1));
        let fitted = lda.fit(&x, &y).unwrap();
        assert_eq!(fitted.scalings().nrows(), 2);
    }

    #[test]
    fn test_lda_pipeline_estimator() {
        use ferrolearn_core::pipeline::PipelineEstimator;

        let (x, y_usize) = linearly_separable_2d();
        let y_f64 = y_usize.mapv(|v| v as f64);
        let lda = LDA::<f64>::new(Some(1));
        let fitted = lda.fit_pipeline(&x, &y_f64).unwrap();
        let preds = fitted.predict_pipeline(&x).unwrap();
        assert_eq!(preds.len(), 8);
    }

    #[test]
    fn test_lda_n_components_getter() {
        let lda = LDA::<f64>::new(Some(2));
        assert_eq!(lda.n_components(), Some(2));
        let lda_none = LDA::<f64>::new(None);
        assert_eq!(lda_none.n_components(), None);
    }

    #[test]
    fn test_lda_priors_builder_default_none() {
        // Default (sklearn `priors=None`, discriminant_analysis.py:359).
        let lda = LDA::<f64>::new(None);
        assert!(lda.priors().is_none());
        // with_priors stores the vector verbatim.
        let lda = lda.with_priors(array![0.7, 0.3]);
        let p = lda.priors().cloned().unwrap_or_default();
        assert_eq!(p.len(), 2);
        assert_abs_diff_eq!(p[0], 0.7, epsilon = 1e-12);
        assert_abs_diff_eq!(p[1], 0.3, epsilon = 1e-12);
    }

    /// Re-oracled (was `test_lda_transform_then_predict_consistent`, which
    /// asserted the OLD nearest-centroid algorithm: `predict ==
    /// argmin ‖transform(x) - projected_mean‖`). The SVD solver's `predict` is
    /// the argmax of the affine `decision_function = X @ coef_.T + intercept_`
    /// (`discriminant_analysis.py:739`), NOT nearest-centroid in projected
    /// space, so this now checks the new contract.
    #[test]
    fn test_lda_predict_matches_decision_argmax() {
        let (x, y) = linearly_separable_2d();
        let lda = LDA::<f64>::new(Some(1));
        let fitted = lda.fit(&x, &y).unwrap();
        let dec = fitted.decision_function(&x).unwrap();
        let preds = fitted.predict(&x).unwrap();
        let n_samples = dec.nrows();
        let n_classes = dec.ncols();
        for i in 0..n_samples {
            let mut best = 0usize;
            let mut best_v = dec[[i, 0]];
            for c in 1..n_classes {
                if dec[[i, c]] > best_v {
                    best_v = dec[[i, c]];
                    best = c;
                }
            }
            assert_eq!(preds[i], fitted.classes()[best]);
        }
    }

    #[test]
    fn test_lda_projected_class_separation() {
        let (x, y) = linearly_separable_2d();
        let lda = LDA::<f64>::new(Some(1));
        let fitted = lda.fit(&x, &y).unwrap();
        let projected = fitted.transform(&x).unwrap();

        // Means of class 0 and class 1 in projected space should be far apart.
        let mean0: f64 = projected
            .rows()
            .into_iter()
            .zip(y.iter())
            .filter(|&(_, label)| *label == 0)
            .map(|(row, _)| row[0])
            .sum::<f64>()
            / 4.0;
        let mean1: f64 = projected
            .rows()
            .into_iter()
            .zip(y.iter())
            .filter(|&(_, label)| *label == 1)
            .map(|(row, _)| row[0])
            .sum::<f64>()
            / 4.0;

        assert!(
            (mean0 - mean1).abs() > 0.5,
            "Projected means should differ, got {mean0} vs {mean1}"
        );
    }

    #[test]
    fn test_lda_transform_known_data() {
        // With perfectly separated 1-D data the centered/whitened transform
        // should still place the two classes on opposite sides.
        let x = Array2::from_shape_vec((4, 1), vec![-2.0, -1.0, 1.0, 2.0]).unwrap();
        let y = array![0usize, 0, 1, 1];
        let lda = LDA::<f64>::new(Some(1));
        let fitted = lda.fit(&x, &y).unwrap();
        let proj = fitted.transform(&x).unwrap();
        let sign0 = proj[[0, 0]].signum();
        let sign1 = proj[[2, 0]].signum();
        assert_ne!(
            sign0 as i32, sign1 as i32,
            "Classes should be on opposite sides"
        );
    }

    #[test]
    fn test_lda_predict_proba_rows_sum_to_one() {
        let (x, y) = three_class_data();
        let lda = LDA::<f64>::new(Some(2));
        let fitted = lda.fit(&x, &y).unwrap();
        let proba = fitted.predict_proba(&x).unwrap();
        assert_eq!(proba.dim(), (9, 3));
        for row in proba.rows() {
            let s: f64 = row.iter().sum();
            assert_abs_diff_eq!(s, 1.0, epsilon = 1e-12);
        }
    }
}