ferrolearn_linear/

linear_svc.rs

//! Linear Support Vector Classifier.
//!
//! This module provides [`LinearSVC`], a liblinear-faithful linear support
//! vector classifier that operates directly in the primal space without the
//! overhead of a kernel function. The fit minimizes the L2-regularized
//! hinge (or squared-hinge) classification objective
//!
//! ```text
//!   min_w   0.5 * ||w||^2  +  C * sum_i  L(y_i, w . x_i)
//! ```
//!
//! with `y_i ∈ {-1, +1}` per one-vs-rest sub-problem (NO `1/n` averaging — the
//! summed loss is scaled by `C`, matching `sklearn/svm/_base.py`
//! `_fit_liblinear`). The solver is liblinear's dual coordinate descent for
//! classification (`solve_l2r_l1l2_svc` in
//! `sklearn/svm/src/liblinear/linear.cpp:819`), which converges to the unique
//! minimizer of the strongly convex objective.
//!
//! Unlike [`SVC`](crate::svm::SVC) with a [`LinearKernel`](crate::svm::LinearKernel),
//! `LinearSVC` avoids computing and caching the full kernel matrix, making it
//! significantly faster for high-dimensional data.
//!
//! # Examples
//!
//! ```
//! use ferrolearn_linear::linear_svc::LinearSVC;
//! use ferrolearn_core::{Fit, Predict};
//! use ndarray::{array, Array2};
//!
//! let x = Array2::from_shape_vec((6, 2), vec![
//!     1.0, 1.0, 1.0, 2.0, 2.0, 1.0,
//!     5.0, 5.0, 5.0, 6.0, 6.0, 5.0,
//! ]).unwrap();
//! let y = array![0usize, 0, 0, 1, 1, 1];
//!
//! let model = LinearSVC::<f64>::new();
//! let fitted = model.fit(&x, &y).unwrap();
//! let preds = fitted.predict(&x).unwrap();
//! assert_eq!(preds.len(), 6);
//! ```
//!
//! ## REQ status
//!
//! Binary (R-DEFER-2): SHIPPED = impl + non-test production consumer + tests +
//! green oracle verification; NOT-STARTED = open blocker `#`. `LinearSVC`/
//! `FittedLinearSVC`/`LinearSVCLoss` are boundary estimator types re-exported at
//! the crate root (`pub use linear_svc::{…}` in `lib.rs`) and registered as the
//! PyO3 `RsLinearSVC` estimator (`ferrolearn-python/src/extras.rs`); under
//! S5/R-DEFER-1 those ARE the non-test production-consumer surface. See
//! `.design/linear/linear_svc.md`.
//!
//! | REQ | Status | Evidence |
//! |---|---|---|
//! | REQ-1 (fit parity — coef_/intercept_ vs liblinear oracle) | SHIPPED | `fn solve_binary_dual` minimizes liblinear's `0.5·‖w‖² + C·Σ L` via the dual CD (`solve_l2r_l1l2_svc`, `linear.cpp:819`); `fn fit` maps `classes_[1]→+1` and extracts `coef_ = w[:n_features]`, `intercept_ = intercept_scaling·w_last` (`_base.py:1240-1245`). Pinned by `tests/divergence_linear_svc_fit.rs::linear_svc_coef_parity` (live oracle `coef_ [[0.12835213611984458, 0.12835213611984475]]`, `intercept_ [-1.1943776585907158]`, C=1.0, squared_hinge, fit_intercept=True). Consumer: `pub use linear_svc::{…}` (`lib.rs`) + `RsLinearSVC` (PyO3). |
//! | REQ-2 (decision_function shape `(n,)` + values) | SHIPPED | `fn decision_function` returns [`DecisionScores::Binary`] = 1-D `X·w + b` for the binary case (sklearn ravels the single-column score to `(n,)`, `linear_model/_base.py:365`) and [`DecisionScores::Multiclass`] `(n, n_classes)` otherwise. Pinned by `tests/divergence_linear_svc_fit.rs::linear_svc_decision_function` (live oracle 1-D `(8,)` values). Consumer: `fn predict` reads the binary scores' sign. |
//! | REQ-3 (predict + classes_) | SHIPPED | `fn predict` uses the sign of the binary decision (`>= 0 → classes_[1]`) / argmax of the OvR scores; `HasClasses::classes` = sorted unique `y` (`classes_ = np.unique(y)`, `_classes.py:311`). The labels are downstream of the liblinear-parity fit and pinned against the live oracle by `linear_svc_predict_parity in tests/divergence_linear_svc_fit.rs` (#620; 8×2 set: `predict [0,0,0,0,1,1,1,1]`, `classes_ [0,1]`). |
//! | REQ-4 (loss {hinge, squared_hinge}) | SHIPPED | The dual CD solves BOTH the true `hinge` (`U=C`, `diag=0`) and `squared_hinge` (`U=∞`, `diag=0.5/C`) optima (`solve_l2r_l1l2_svc`, `linear.cpp:849-858`). The `hinge` optimum is pinned against the live oracle by `linear_svc_hinge_coef_parity in tests/divergence_linear_svc_fit.rs` (#621; 8×2 set, `loss='hinge'`, `C=1.0`: `coef_ [[0.15384615383852776, 0.15384615383915584]]`, `intercept_ [-1.4615384615168394]`). |
//! | REQ-5 (penalty {l1, l2}) | SHIPPED | `LinearSVC<F>` exposes `pub penalty: LinearSVCPenalty` (default `L2`) + `#[must_use] with_penalty`. `penalty=l1` routes to `fn solve_binary_l1r_l2` — liblinear's feature-major coordinate descent (`solve_l1r_l2_svc`, `linear.cpp:1467`, solver type 5, `_base.py:1014`) minimizing `‖w‖₁ + C·Σ max(0,1−yf)²` (sparse `coef_`); `penalty=l2` keeps `fn solve_binary_dual`. `liblinear_solver_type` rejects `('l1','hinge')` (`_base.py:1013`). Pinned by `test_l1_penalty_smoke` (live oracle 8×2 `l1,squared_hinge,dual=False,C=1`: `coef_ [[0.1283185834966579, 0.12831858464059265]]`, `intercept_ [-1.2079646017762715]`; ferrolearn lands within ~1.2e-9) + `test_unsupported_combinations_rejected`. Consumer: `pub use linear_svc::{…}` (`lib.rs`) + `RsLinearSVC` (PyO3). |
//! | REQ-6 (multi_class {ovr, crammer_singer}) | SHIPPED | `LinearSVC<F>` exposes `pub multi_class: MultiClass` (`Ovr`/`CrammerSinger`, default `Ovr`, `_classes.py:239`) + `#[must_use] with_multi_class`. `multi_class=CrammerSinger` selects liblinear solver type 4 REGARDLESS of penalty/loss/dual (`_get_liblinear_solver_type`, `_base.py:1017,1020-1021`), so `fn fit` short-circuits to `fn fit_crammer_singer` (ignoring penalty/loss/dual) which runs ONE joint solve `fn solve_crammer_singer` — a faithful transcription of `Solver_MCSVM_CS` (`linear.cpp:510`, the class at `:493-787`): flattened `w[feature*nr_class + m]`, `alpha[i*nr_class + m]` with `Σ_m alpha=0`, `C[i] = weighted_C[y_i]`, per-sample shrinking (`active_size_i`/`alpha_index`/`be_shrunk`), the simplex inner solve `fn cs_solve_sub_problem` (sort-descending breakpoint, `linear.cpp:541-564`), and the two-level `eps_shrink = max(10·tol, 1)` stopping (`linear.cpp:738-753`). Extraction: `coef_[m][feature] = w[feature*nr_class + m]`, `intercept_[m] = intercept_scaling·w[n_features*nr_class + m]` (`_base.py:1240-1245`). BINARY: collapse to a single weight vector `coef_ = row_1 − row_0`, `intercept_ = int_1 − int_0` (`_classes.py:340-344`), `is_binary=true`. ferrolearn sweeps natural order (no `bounded_rand_int` shuffle, `linear.cpp:629`); the CS optimum is unique so the limit is identical (documented RNG/shrink-path boundary). Pinned by `test_crammer_singer_smoke in linear_svc.rs` (live oracle 3-class set `coef [[-0.06762,-0.24341],[0.30048,0.02171],[-0.23286,0.22171]]`, `int [0.91078,-0.62206,-0.28873]`, predict all-correct; binary 8×2 collapse `coef [[0.15504,0.15504]]`, `int [-1.48062]`; within 1e-2). The rigorous oracle pin in `tests/divergence_linear_svc_fit.rs` is the critic's next step. Consumer: `pub use linear_svc::{…}` (`lib.rs`) + `RsLinearSVC` (PyO3). |
//! | REQ-7 (fit_intercept + intercept_scaling) | SHIPPED | `LinearSVC<F>` exposes `pub fit_intercept: bool` (default true) + `pub intercept_scaling: F` (default 1.0) + `#[must_use]` builders. When fitting an intercept the design matrix is augmented with a penalized constant column = `intercept_scaling`, and `intercept_ = intercept_scaling·w_last` (`_base.py:1188-1198,:1240-1245`); `intercept_scaling > 0` is validated. Pinned by `linear_svc_coef_parity` + module `test_fit_intercept_false_zero_intercept`/`test_invalid_intercept_scaling`. |
//! | REQ-8 (dual param) | SHIPPED | `LinearSVC<F>` exposes `pub dual: DualMode` (default `Auto`) + `#[must_use] with_dual`. `fn resolve_dual` resolves `Auto→bool` (`_validate_dual_parameter`, `_classes.py:13-29`: `n<f`→prefer dual, else→prefer primal, with fallback) against `fn liblinear_solver_type` (the `_get_liblinear_solver_type` matrix, `_base.py:995-1018`), and `fn fit` validates the resolved combo (`hinge+dual=false`, `l1+dual=true`, `l1+hinge` all rejected → `FerroError::InvalidParameter`). R-DEV-7: the resolved `dual` is **observably immaterial for `penalty=l2`** — the l2 dual CD and l2 primal minimize the same `0.5·‖w‖² + C·Σ L` and reach the same `coef_`/`intercept_`, so `penalty=l2` keeps `fn solve_binary_dual` regardless of `dual`. Pinned by `test_unsupported_combinations_rejected` + `test_dual_auto_resolution`. Consumer: `pub use linear_svc::{…}` (`lib.rs`) + `RsLinearSVC` (PyO3). |
//! | REQ-9 (class_weight) | SHIPPED | `LinearSVC<F>` exposes `pub class_weight: ClassWeight<F>` (`None`/`Balanced`/`Explicit`, default `None`) + `#[must_use] with_class_weight`. `fn compute_class_weight` (mirroring `sklearn.utils.compute_class_weight`, `class_weight.py:63-81`, as called at `_base.py:1179`) expands per-class weights; `fn fit` scales `C` per class: binary `cp = C·weights[idx(classes[1])]`, `cn = C·weights[idx(classes[0])]` (`train_one(Cp=weighted_C[1], Cn=weighted_C[0])`, `linear.cpp:2543-2551`), OvR class `k` `cp = C·weights[k]`, `cn = C` base (the negative rest is UNWEIGHTED, `linear.cpp:2559-2571`). `SolverConfig` now carries `(cp, cn)`; `solve_binary_dual`/`solve_binary_l1r_l2` apply the per-sample `C_[i] = (y_i>0?cp:cn)` (`diag[i]`/`upper_bound[i]`/`C[i]`, `linear.cpp:843-858`,`:1504-1509`). When `cp == cn` (no class_weight) the math is identical to before (the 9 divergence pins stay green). Pinned by `test_class_weight_smoke in linear_svc.rs` (live oracle 8×2 imbalanced set, `squared_hinge,dual=True,C=1.0`: `None coef [[0.10056,0.15957]] int [-1.26346]`; `balanced coef [[0.09937,0.16666]] int [-1.21320]` weights `[0.6667,2.0]`; `{0:1,1:5} coef [[0.11059,0.17164]] int [-1.29547]`; ferrolearn within 1e-2). The rigorous oracle pin in `tests/divergence_linear_svc_fit.rs` is the critic's next step. Consumer: `pub use linear_svc::{…}` (`lib.rs`) + `RsLinearSVC` (PyO3). |
//! | REQ-10 (C-scaling convention) | SHIPPED | the `c / n_f` division is removed; the dual CD uses `upper_bound = C` (hinge) / `diag = 0.5/C` (squared_hinge), so `coef_` tracks `C` like liblinear. Pinned by `linear_svc_coef_c_dependence` (C=0.1 → `0.0784651864625997`, C=1.0 → `0.12835213611984458`). |
//! | REQ-11 (n_iter_/n_features_in_ + param validation) | SHIPPED | `fn n_features_in` (returns the stored `n_features`, set by `_validate_data`, `_classes.py:302`) and `fn n_iter` (the max dual-CD outer-iteration count across the binary/OvR fits, `n_iter_ = n_iter_.max().item()`, `_classes.py:338`) on `FittedLinearSVC`; `fn fit` validates `tol > 0` (`Interval(Real, 0.0, None, closed="neither")`, `_classes.py:237`). Pinned by `linear_svc_attrs_and_tol_validation in tests/divergence_linear_svc_fit.rs` (#627). `n_features_in_` (oracle `2`) and the `tol <= 0` reject are exact; `n_iter_` is the documented shuffle-path RNG boundary (ferrolearn sweeps natural order, sklearn's liblinear shuffles `index` each sweep, cf. SGD), so the pin bounds `n_iter` in `[1, max_iter]` rather than exact-matching. |
//! | REQ-12 (ferray substrate) | NOT-STARTED | open prereq blocker #628. Imports `ndarray`, not `ferray-core`/`ferray::linalg` (R-SUBSTRATE). |
//! | REQ-13 (non-finite input rejected) | SHIPPED | `fn fit` rejects any NaN/+/-inf in X BEFORE the dual coordinate-descent solve with `FerroError::InvalidParameter`, mirroring sklearn's `_validate_data(force_all_finite=True)` (`svm/_classes.py:302`, the liblinear `_fit_liblinear` path; cf. `svm/_base.py:190`) → `ValueError("Input X contains NaN.")` / `"... contains infinity ..."`. `y` is `Array1<usize>` (integer labels), finite by type, and ferrolearn's `LinearSVC::fit` takes no `sample_weight` argument, so X is the only runtime check. The guard sits ahead of both the OvR and Crammer-Singer dispatch, 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): `LinearSVC().fit` raises `ValueError` for NaN/+inf/-inf in X (`tests/divergence_linear_nonfinite_batch4.rs::linsvc_*`). Non-test consumer: `pub use linear_svc::{…}` (`lib.rs`) + `RsLinearSVC` (PyO3). (#2263) |

use ferrolearn_core::error::FerroError;
use ferrolearn_core::introspection::{HasClasses, HasCoefficients};
use ferrolearn_core::traits::{Fit, Predict};
use ndarray::{Array1, Array2, ScalarOperand};
use num_traits::Float;

/// Penalty (regularizer) norm for [`LinearSVC`].
///
/// Mirrors `sklearn.svm.LinearSVC`'s `penalty` parameter
/// (`sklearn/svm/_classes.py:51-54`): `'l2'` (default) is the standard SVC
/// `0.5·‖w‖²` regularizer; `'l1'` is the `‖w‖₁` regularizer, which yields a
/// sparse `coef_`. The penalty interacts with [`LinearSVCLoss`] / [`DualMode`]
/// via liblinear's solver-selection matrix (`_get_liblinear_solver_type`,
/// `_base.py:1011-1018`): `'l1'` is only supported with `squared_hinge` +
/// `dual=False` (solver type 5).
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum LinearSVCPenalty {
    /// `‖w‖₁` regularizer — sparse `coef_`. Only valid with `squared_hinge` +
    /// `dual=false` (liblinear solver type 5, `_base.py:1014`).
    L1,
    /// `0.5·‖w‖²` regularizer (default), the standard SVC penalty.
    #[default]
    L2,
}

/// Dual / primal optimization-problem selector for [`LinearSVC`].
///
/// Mirrors `sklearn.svm.LinearSVC`'s `dual` parameter
/// (`sklearn/svm/_classes.py:62-71`, default `"auto"`). `Auto` resolves to a
/// concrete `bool` via `_validate_dual_parameter` (`_classes.py:13-29`): when
/// `n_samples < n_features` it prefers `dual=true` (falling back to `false` if
/// that penalty×loss combination has no dual solver); otherwise it prefers
/// `dual=false` (falling back to `true` if there is no primal solver, e.g.
/// `hinge`). The resolved `dual` selects the liblinear solver type
/// (`_get_liblinear_solver_type`, `_base.py:1011-1018`).
///
/// Under R-DEV-7 the resolved `dual` is **observably immaterial for
/// `penalty=l2`**: the l2 dual coordinate descent and the l2 primal both
/// minimize the same strongly convex `0.5·‖w‖² + C·Σ L` and reach the same
/// `coef_`/`intercept_`. It is load-bearing only for the unsupported-combination
/// rejects and for selecting the genuinely different `l1` primal solver.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum DualMode {
    /// Resolve to `true`/`false` automatically via `_validate_dual_parameter`
    /// (`_classes.py:13-29`); the default.
    #[default]
    Auto,
    /// Solve the dual optimization problem.
    True,
    /// Solve the primal optimization problem.
    False,
}

/// Multiclass strategy for [`LinearSVC`].
///
/// Mirrors `sklearn.svm.LinearSVC`'s `multi_class` parameter
/// (`sklearn/svm/_classes.py:239`, constraint `{"ovr", "crammer_singer"}`,
/// default `"ovr"`). `Ovr` trains one binary classifier per class (the default,
/// using the penalty/loss/dual solver matrix); `CrammerSinger` runs the joint
/// Crammer-Singer multiclass SVM solver (`Solver_MCSVM_CS`,
/// `liblinear/linear.cpp:493-787`, solver type 4).
///
/// When `CrammerSinger` is selected, `_get_liblinear_solver_type` returns 4
/// **regardless of penalty/loss/dual** (`_base.py:1017,1020-1021`), so the
/// penalty/loss/dual parameters are ignored and the joint solver runs.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum MultiClass {
    /// One-vs-rest: one binary sub-problem per class (the default).
    #[default]
    Ovr,
    /// Joint Crammer-Singer multiclass SVM (`Solver_MCSVM_CS`, solver type 4,
    /// `linear.cpp:493-787`). Ignores penalty/loss/dual (`_base.py:1017`).
    CrammerSinger,
}

/// Loss function for [`LinearSVC`].
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum LinearSVCLoss {
    /// Standard hinge loss: `max(0, 1 - y * f(x))`. liblinear solver
    /// `L2R_L1LOSS_SVC_DUAL` (type 3): box `0 <= alpha <= C`, `diag = 0`.
    Hinge,
    /// Squared hinge loss: `max(0, 1 - y * f(x))^2` (default). liblinear solver
    /// `L2R_L2LOSS_SVC_DUAL` (type 1): box `0 <= alpha <= +inf`,
    /// `diag = 0.5 / C`.
    SquaredHinge,
}

/// Per-class weighting strategy for [`LinearSVC`].
///
/// Mirrors `sklearn.svm.LinearSVC`'s `class_weight` parameter
/// (`sklearn/svm/_classes.py:118-124`, constraint `{None, dict, 'balanced'}`):
/// it scales the inverse-regularization `C` per class so the effective penalty
/// for class `i` is `class_weight[i]·C` (`compute_class_weight`,
/// `sklearn/svm/_base.py:1179`; `weighted_C[i] = C·class_weight[i]`,
/// `liblinear/linear.cpp:2496-2507`). The expanded per-class weights are
/// computed by [`compute_class_weight`] following
/// `sklearn.utils.compute_class_weight` semantics
/// (`sklearn/utils/class_weight.py:63-81`).
///
/// This mirrors `ferrolearn_linear::sgd::ClassWeight` for cross-estimator
/// consistency, but is defined locally (no cross-import of `sgd` internals).
#[derive(Debug, Clone, Default)]
pub enum ClassWeight<F> {
    /// Uniform weights (all classes weighted `1.0`). The default
    /// (`class_weight=None`, `class_weight.py:63-65`).
    #[default]
    None,
    /// Balanced weights `n_samples / (n_classes · count_c)` per class `c`,
    /// matching `sklearn.utils.compute_class_weight("balanced", ...)`
    /// (`class_weight.py:66-74`).
    Balanced,
    /// Explicit class-label -> weight map. Classes absent from the map default
    /// to `1.0`, matching the dict branch of `compute_class_weight`
    /// (`class_weight.py:75-81`).
    Explicit(Vec<(usize, F)>),
}

/// Compute the expanded per-class weight vector aligned to `classes`
/// (sorted ascending, matching sklearn's `classes_ = np.unique(y)`).
///
/// Faithful to `sklearn.utils.compute_class_weight`
/// (`sklearn/utils/class_weight.py:63-81`), as called by `_fit_liblinear`
/// (`compute_class_weight(class_weight, classes=classes_, y=y)`,
/// `sklearn/svm/_base.py:1179`):
/// - `None` -> all `1.0` (`:63-65`).
/// - `Balanced` -> `n_samples / (n_classes · count_c)` per class `c`,
///   where `count_c` is the number of samples with label `c` (`:66-74`).
/// - `Explicit(map)` -> `1.0` default, overridden by the map entries matched by
///   class label (`:75-81`).
///
/// `classes` is the sorted unique label set; `y` is the per-sample label array.
/// Mirrors `ferrolearn_linear::sgd::compute_class_weight` exactly.
fn compute_class_weight<F: Float>(cw: &ClassWeight<F>, classes: &[usize], y: &[usize]) -> Vec<F> {
    match cw {
        ClassWeight::None => vec![F::one(); classes.len()],
        ClassWeight::Balanced => {
            // `recip_freq = len(y) / (n_classes * bincount(y_ind))`
            // (`class_weight.py:73`), indexed per class.
            let n_samples = F::from(y.len()).unwrap_or_else(F::zero);
            let n_classes = F::from(classes.len()).unwrap_or_else(F::one);
            classes
                .iter()
                .map(|&c| {
                    let count = y.iter().filter(|&&label| label == c).count();
                    let count_f = F::from(count).unwrap_or_else(F::one);
                    if count_f > F::zero() {
                        n_samples / (n_classes * count_f)
                    } else {
                        F::one()
                    }
                })
                .collect()
        }
        ClassWeight::Explicit(map) => classes
            .iter()
            .map(|&c| {
                map.iter()
                    .find(|(label, _)| *label == c)
                    .map_or_else(F::one, |(_, w)| *w)
            })
            .collect(),
    }
}

/// Confidence scores returned by [`FittedLinearSVC::decision_function`].
///
/// Mirrors `sklearn.svm.LinearSVC.decision_function`
/// (`LinearClassifierMixin.decision_function`,
/// `sklearn/linear_model/_base.py:341-365`): the binary case collapses the
/// single-column score matrix to a 1-D `(n_samples,)` array
/// (`return xp.reshape(scores, (-1,)) if scores.shape[1] == 1 else scores`,
/// `_base.py:365`), the multiclass case returns `(n_samples, n_classes)`.
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum DecisionScores<F> {
    /// Binary scores `X · w + b` for the positive class (`classes_[1]`), shape
    /// `(n_samples,)`. `> 0` predicts `classes_[1]`.
    Binary(Array1<F>),
    /// One-vs-rest scores, shape `(n_samples, n_classes)`. The argmax of each
    /// row agrees with [`Predict`].
    Multiclass(Array2<F>),
}

impl<F: Clone> DecisionScores<F> {
    /// Number of samples scored (the leading axis length in both variants).
    #[must_use]
    pub fn n_samples(&self) -> usize {
        match self {
            DecisionScores::Binary(v) => v.len(),
            DecisionScores::Multiclass(m) => m.nrows(),
        }
    }

    /// Borrow the binary 1-D scores, if this is the binary case.
    #[must_use]
    pub fn as_binary(&self) -> Option<&Array1<F>> {
        match self {
            DecisionScores::Binary(v) => Some(v),
            DecisionScores::Multiclass(_) => None,
        }
    }

    /// Borrow the multiclass `(n_samples, n_classes)` scores, if this is the
    /// multiclass case.
    #[must_use]
    pub fn as_multiclass(&self) -> Option<&Array2<F>> {
        match self {
            DecisionScores::Multiclass(m) => Some(m),
            DecisionScores::Binary(_) => None,
        }
    }
}

/// Linear Support Vector Classifier (liblinear dual CD).
///
/// Solves the L2-regularized hinge or squared-hinge objective
/// `0.5*||w||^2 + C * sum_i L(y_i, w.x_i)` via liblinear's dual coordinate
/// descent. Supports binary and multiclass (one-vs-rest) classification.
/// Mirrors `sklearn.svm.LinearSVC`.
///
/// # Type Parameters
///
/// - `F`: The floating-point type (`f32` or `f64`).
#[derive(Debug, Clone)]
pub struct LinearSVC<F> {
    /// Inverse regularization strength. Larger values allow more
    /// misclassification. Must be strictly positive.
    pub c: F,
    /// Maximum number of dual coordinate descent iterations.
    pub max_iter: usize,
    /// Convergence tolerance on the projected-gradient span.
    pub tol: F,
    /// Loss function to use.
    pub loss: LinearSVCLoss,
    /// Regularizer norm (`l1` / `l2`). Default `l2`. `l1` is only supported with
    /// `squared_hinge` + `dual=false` (liblinear solver type 5,
    /// `_get_liblinear_solver_type`, `_base.py:1014`).
    pub penalty: LinearSVCPenalty,
    /// Dual / primal selector. Default `Auto`, resolved via
    /// `_validate_dual_parameter` (`_classes.py:13-29`). Observably immaterial
    /// for `penalty=l2` (R-DEV-7 dual-invariance), load-bearing for the
    /// unsupported-combination rejects and the `l1` primal solver.
    pub dual: DualMode,
    /// Whether to fit an intercept. When `true`, the design matrix is augmented
    /// with a synthetic constant column equal to [`intercept_scaling`]; the
    /// augmented weight is penalized like any feature (liblinear convention).
    ///
    /// [`intercept_scaling`]: Self::intercept_scaling
    pub fit_intercept: bool,
    /// Value of the synthetic intercept feature when [`fit_intercept`] is
    /// `true`. Must be strictly positive. `intercept_ = intercept_scaling *
    /// w_last`.
    ///
    /// [`fit_intercept`]: Self::fit_intercept
    pub intercept_scaling: F,
    /// Per-class scaling of `C`. Default [`ClassWeight::None`] (all classes
    /// weighted `1.0`). The effective penalty for class `i` is
    /// `class_weight[i]·C` (`compute_class_weight`, `_base.py:1179`;
    /// `weighted_C[i] = C·class_weight[i]`, `linear.cpp:2496-2507`).
    pub class_weight: ClassWeight<F>,
    /// Multiclass strategy. Default [`MultiClass::Ovr`] (one-vs-rest).
    /// [`MultiClass::CrammerSinger`] runs the joint Crammer-Singer solver,
    /// ignoring penalty/loss/dual (`_base.py:1017`, `_classes.py:239`).
    pub multi_class: MultiClass,
}

impl<F: Float> LinearSVC<F> {
    /// Create a new `LinearSVC` with scikit-learn's default settings.
    ///
    /// Defaults (matching `sklearn.svm.LinearSVC`, `_classes.py`):
    /// `C = 1.0`, `max_iter = 1000`, `tol = 1e-4`, `loss = SquaredHinge`,
    /// `penalty = L2`, `dual = Auto`, `fit_intercept = true`,
    /// `intercept_scaling = 1.0`, `class_weight = None`, `multi_class = Ovr`.
    #[must_use]
    pub fn new() -> Self {
        // 1e-4/1.0 are exactly representable in f32/f64; the defensive fallback
        // for `from(1e-4)` is never taken (no `.unwrap()` in lib code).
        let one = F::one();
        Self {
            c: one,
            max_iter: 1000,
            tol: F::from(1e-4).unwrap_or_else(|| {
                let ten = F::from(10).unwrap_or(one);
                one / (ten * ten * ten * ten)
            }),
            loss: LinearSVCLoss::SquaredHinge,
            penalty: LinearSVCPenalty::L2,
            dual: DualMode::Auto,
            fit_intercept: true,
            intercept_scaling: one,
            class_weight: ClassWeight::None,
            multi_class: MultiClass::Ovr,
        }
    }

    /// Set the regularization parameter C.
    #[must_use]
    pub fn with_c(mut self, c: F) -> Self {
        self.c = c;
        self
    }

    /// Set the maximum number of iterations.
    #[must_use]
    pub fn with_max_iter(mut self, max_iter: usize) -> Self {
        self.max_iter = max_iter;
        self
    }

    /// Set the convergence tolerance.
    #[must_use]
    pub fn with_tol(mut self, tol: F) -> Self {
        self.tol = tol;
        self
    }

    /// Set the loss function.
    #[must_use]
    pub fn with_loss(mut self, loss: LinearSVCLoss) -> Self {
        self.loss = loss;
        self
    }

    /// Set the penalty (regularizer) norm (sklearn `penalty`). `l1` requires
    /// `squared_hinge` + `dual=false` (`_base.py:1014`).
    #[must_use]
    pub fn with_penalty(mut self, penalty: LinearSVCPenalty) -> Self {
        self.penalty = penalty;
        self
    }

    /// Set the dual / primal selector (sklearn `dual`). `Auto` (default)
    /// resolves via `_validate_dual_parameter` (`_classes.py:13-29`).
    #[must_use]
    pub fn with_dual(mut self, dual: DualMode) -> Self {
        self.dual = dual;
        self
    }

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

    /// Set the intercept scaling (sklearn `intercept_scaling`). Must be
    /// strictly positive when `fit_intercept` is `true`.
    #[must_use]
    pub fn with_intercept_scaling(mut self, intercept_scaling: F) -> Self {
        self.intercept_scaling = intercept_scaling;
        self
    }

    /// Set the per-class `C` scaling (sklearn `class_weight`,
    /// `_classes.py:118-124`). [`ClassWeight::None`] (default) leaves every
    /// class at `1.0`; [`ClassWeight::Balanced`] uses
    /// `n_samples / (n_classes · count_c)`; [`ClassWeight::Explicit`] takes a
    /// `(label, weight)` map (unlisted classes default to `1.0`).
    #[must_use]
    pub fn with_class_weight(mut self, class_weight: ClassWeight<F>) -> Self {
        self.class_weight = class_weight;
        self
    }

    /// Set the multiclass strategy (sklearn `multi_class`, `_classes.py:239`).
    /// [`MultiClass::Ovr`] (default) trains one-vs-rest binary classifiers;
    /// [`MultiClass::CrammerSinger`] runs the joint Crammer-Singer solver
    /// (ignoring penalty/loss/dual, `_base.py:1017`).
    #[must_use]
    pub fn with_multi_class(mut self, multi_class: MultiClass) -> Self {
        self.multi_class = multi_class;
        self
    }
}

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

/// Fitted Linear Support Vector Classifier.
///
/// Stores the learned weight vectors, intercepts, and class labels.
/// For binary classification a single weight vector is stored; for
/// multiclass, one per class (one-vs-rest).
#[derive(Debug, Clone)]
pub struct FittedLinearSVC<F> {
    /// Weight vectors: one per binary sub-problem.
    /// Binary: `[w]`, Multiclass: `[w_0, w_1, ..., w_{k-1}]`.
    weight_vectors: Vec<Array1<F>>,
    /// Intercept for each sub-problem.
    intercepts: Vec<F>,
    /// Sorted unique class labels.
    classes: Vec<usize>,
    /// Whether this is a binary problem.
    is_binary: bool,
    /// Number of features.
    n_features: usize,
    /// Maximum dual-CD outer-iteration count across sub-problem fits
    /// (`n_iter_ = n_iter_.max().item()`, `_classes.py:338`).
    n_iter: usize,
}

impl<F: Float> FittedLinearSVC<F> {
    /// Returns the weight vectors (one per binary sub-problem).
    #[must_use]
    pub fn weight_vectors(&self) -> &[Array1<F>] {
        &self.weight_vectors
    }

    /// Returns the intercepts (one per binary sub-problem).
    #[must_use]
    pub fn intercepts(&self) -> &[F] {
        &self.intercepts
    }

    /// Number of features seen during fit (`n_features_in_`).
    ///
    /// Mirrors sklearn's `n_features_in_`, set by `_validate_data`
    /// (`sklearn/svm/_classes.py:302`); equals `X.ncols()`.
    #[must_use]
    pub fn n_features_in(&self) -> usize {
        self.n_features
    }

    /// Maximum number of dual coordinate-descent outer iterations across the
    /// (binary or one-vs-rest) sub-problem fits.
    ///
    /// Mirrors sklearn's `n_iter_ = n_iter_.max().item()`
    /// (`sklearn/svm/_classes.py:338`). The exact value is shuffle-path
    /// dependent (sklearn's liblinear shuffles the active index each sweep;
    /// ferrolearn sweeps natural order), so it is bounded in `[1, max_iter]`
    /// rather than exact-matching the oracle.
    #[must_use]
    pub fn n_iter(&self) -> usize {
        self.n_iter
    }
}

impl<F: Float + ScalarOperand + Send + Sync + 'static> FittedLinearSVC<F> {
    /// Raw signed distance from the decision boundary. Mirrors sklearn
    /// `LinearSVC.decision_function`.
    ///
    /// Binary: [`DecisionScores::Binary`] of shape `(n_samples,)` containing
    /// `X @ w + b` for the positive class (`classes_[1]`); sklearn ravels the
    /// single-column score to 1-D (`linear_model/_base.py:365`).
    /// Multiclass: [`DecisionScores::Multiclass`] of shape
    /// `(n_samples, n_classes)` of one-vs-rest scores; the 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<DecisionScores<F>, FerroError> {
        let n_features = x.ncols();
        if n_features != self.n_features {
            return Err(FerroError::ShapeMismatch {
                expected: vec![self.n_features],
                actual: vec![n_features],
                context: "number of features must match fitted model".into(),
            });
        }
        let n_samples = x.nrows();
        if self.is_binary {
            // sklearn collapses the single-column binary score matrix to a 1-D
            // (n_samples,) array (`linear_model/_base.py:365`).
            let scores = x.dot(&self.weight_vectors[0]) + self.intercepts[0];
            Ok(DecisionScores::Binary(scores))
        } else {
            let n_classes = self.classes.len();
            let mut out = Array2::<F>::zeros((n_samples, n_classes));
            for c in 0..n_classes {
                for i in 0..n_samples {
                    out[[i, c]] = x.row(i).dot(&self.weight_vectors[c]) + self.intercepts[c];
                }
            }
            Ok(DecisionScores::Multiclass(out))
        }
    }
}

/// Resolved per-fit solver configuration passed to [`solve_binary_dual`]
/// (groups the dual-CD knobs to keep the solver signature small).
#[derive(Debug, Clone, Copy)]
struct SolverConfig<F> {
    /// Per-sample penalty for the positive (`y_i > 0`) group: `Cp = C·w[+]`
    /// (`train_one(Cp, Cn)`, `linear.cpp:2543-2571`; `C_[i] = (y_i>0 ? Cp :
    /// Cn)`, `linear.cpp:843-858`, `:1504-1509`).
    cp: F,
    /// Per-sample penalty for the negative (`y_i <= 0`) group: `Cn = C·w[-]`
    /// (binary) or the base `C` (multiclass OvR; the negative group is the
    /// unweighted rest, `linear.cpp:2559-2571`).
    cn: F,
    /// Maximum dual-CD outer iterations.
    max_iter: usize,
    /// Projected-gradient span stopping tolerance.
    tol: F,
    /// Hinge / squared-hinge loss selector.
    loss: LinearSVCLoss,
    /// Whether the design matrix is augmented with a penalized bias column.
    fit_intercept: bool,
    /// Value of the synthetic bias column when `fit_intercept` is set.
    intercept_scaling: F,
}

/// Solve a single binary L2-SVM via liblinear's dual coordinate descent.
///
/// Minimizes `0.5 * ||w||^2 + C * sum_i L(y_i, w.x_i)` (NO `1/n` averaging,
/// matching `sklearn/svm/_base.py` `_fit_liblinear`) with `y_i ∈ {-1, +1}`.
/// The augmented weight vector has length `w_size`; when `fit_intercept` is
/// set, `w_size = n_features + 1` and the trailing weight multiplies the
/// synthetic constant column `intercept_scaling` (penalized like any feature).
///
/// This is liblinear's `solve_l2r_l1l2_svc` (`linear.cpp:819`): the dual is
///
/// ```text
///   min_alpha  0.5 * alpha^T (Q + diag) alpha  -  e^T alpha,
///     s.t.     0 <= alpha_i <= U,
/// ```
///
/// where `Q_ij = y_i y_j x_i.x_j`, `w = sum_i alpha_i y_i x_i`,
/// `QD[i] = diag + ||x_i||^2`. **hinge** (`L2R_L1LOSS_SVC_DUAL`): `diag = 0`,
/// `U = C` (`linear.cpp:852-858`). **squared_hinge** (`L2R_L2LOSS_SVC_DUAL`):
/// `diag = 0.5/C`, `U = +inf` (`linear.cpp:849-850`).
///
/// Returns `(w_augmented, n_iter, converged)`.
fn solve_binary_dual<F: Float + 'static>(
    x: &Array2<F>,
    y_signed: &Array1<F>,
    cfg: &SolverConfig<F>,
) -> (Vec<F>, usize, bool) {
    let SolverConfig {
        cp,
        cn,
        max_iter,
        tol,
        loss,
        fit_intercept,
        intercept_scaling,
    } = *cfg;

    let (n_samples, n_features) = x.dim();
    let w_size = if fit_intercept {
        n_features + 1
    } else {
        n_features
    };

    let inf = F::infinity();
    let two = F::one() + F::one();
    let half = F::one() / two;
    let tiny = F::from(1.0e-12).unwrap_or_else(F::epsilon);

    // Per-sample penalty `C_[i] = (y_i > 0 ? Cp : Cn)` (`linear.cpp:843-858`,
    // `GETI(i) ≡ i`); `class_weight` makes Cp/Cn differ. Per-sample diag /
    // upper_bound follow the solver type: squared_hinge → `diag[i] = 0.5/C_[i]`,
    // `U[i] = +inf`; hinge → `diag[i] = 0`, `U[i] = C_[i]`.
    let mut diag = vec![F::zero(); n_samples];
    let mut upper_bound = vec![inf; n_samples];
    for i in 0..n_samples {
        let c_i = if y_signed[i] > F::zero() { cp } else { cn };
        match loss {
            LinearSVCLoss::Hinge => {
                diag[i] = F::zero();
                upper_bound[i] = c_i;
            }
            LinearSVCLoss::SquaredHinge => {
                diag[i] = half / c_i;
                upper_bound[i] = inf;
            }
        }
    }

    // QD[i] = diag[i] + ||x_i||^2 (including the augmented bias column).
    let mut qd = vec![F::zero(); n_samples];
    for (i, qd_i) in qd.iter_mut().enumerate() {
        let mut acc = diag[i];
        let row = x.row(i);
        for &v in row.iter() {
            acc = acc + v * v;
        }
        if fit_intercept {
            acc = acc + intercept_scaling * intercept_scaling;
        }
        *qd_i = acc;
    }

    let mut alpha = vec![F::zero(); n_samples];
    let mut w = vec![F::zero(); w_size];
    // alpha starts at 0 so w starts at 0; nothing to accumulate.

    let mut index: Vec<usize> = (0..n_samples).collect();
    let mut active_size = n_samples;
    let mut pgmax_old = inf;
    let mut pgmin_old = -inf;

    // w . x_i over the (augmented) design matrix.
    let dot_w_xi = |w: &[F], i: usize| -> F {
        let mut acc = F::zero();
        let row = x.row(i);
        for (j, &v) in row.iter().enumerate() {
            acc = acc + v * w[j];
        }
        if fit_intercept {
            acc = acc + intercept_scaling * w[n_features];
        }
        acc
    };

    let mut n_iter: usize = 0;
    let mut converged = false;

    // liblinear shuffles `index` each sweep; the minimizer is unique so order
    // only affects the path, not the limit. We sweep in natural order for
    // determinism (no RNG), reaching the same converged optimum.
    for iter in 0..max_iter {
        n_iter = iter + 1;
        let mut pgmax_new = -inf;
        let mut pgmin_new = inf;

        let mut s = 0;
        while s < active_size {
            let i = index[s];
            let yi = y_signed[i];

            // G = y_i*(w.x_i) - 1 + diag_i*alpha_i (`linear.cpp:909-921`).
            let g = yi * dot_w_xi(&w, i) - F::one() + diag[i] * alpha[i];

            // Projected gradient + shrinking (`linear.cpp:923-949`).
            let mut pg = F::zero();
            if alpha[i] == F::zero() {
                if g > pgmax_old {
                    active_size -= 1;
                    index.swap(s, active_size);
                    continue; // re-process the swapped-in element at `s`
                } else if g < F::zero() {
                    pg = g;
                }
            } else if alpha[i] == upper_bound[i] {
                if g < pgmin_old {
                    active_size -= 1;
                    index.swap(s, active_size);
                    continue;
                } else if g > F::zero() {
                    pg = g;
                }
            } else {
                pg = g;
            }

            if pg > pgmax_new {
                pgmax_new = pg;
            }
            if pg < pgmin_new {
                pgmin_new = pg;
            }

            if pg.abs() > tiny {
                let alpha_old = alpha[i];
                // alpha_i <- clamp(alpha_i - G/QD[i], 0, U) (`linear.cpp:957`).
                let mut new_alpha = alpha[i] - g / qd[i];
                if new_alpha < F::zero() {
                    new_alpha = F::zero();
                } else if new_alpha > upper_bound[i] {
                    new_alpha = upper_bound[i];
                }
                alpha[i] = new_alpha;
                let d = (alpha[i] - alpha_old) * yi;
                if d != F::zero() {
                    let row = x.row(i);
                    for (j, &v) in row.iter().enumerate() {
                        w[j] = w[j] + d * v;
                    }
                    if fit_intercept {
                        w[n_features] = w[n_features] + d * intercept_scaling;
                    }
                }
            }

            s += 1;
        }

        // Stopping: PGmax_new - PGmin_new <= tol on the full set
        // (`linear.cpp:972-990`). Absolute tol (the dual SVC solver receives
        // `param->eps` directly, `linear.cpp:2364`).
        if pgmax_new - pgmin_new <= tol {
            if active_size == n_samples {
                converged = true;
                break;
            }
            active_size = n_samples;
            pgmax_old = inf;
            pgmin_old = -inf;
            continue;
        }

        pgmax_old = pgmax_new;
        pgmin_old = pgmin_new;
        if pgmax_old <= F::zero() {
            pgmax_old = inf;
        }
        if pgmin_old >= F::zero() {
            pgmin_old = -inf;
        }
    }

    (w, n_iter, converged)
}

/// Resolve `(penalty, loss, dual)` to a liblinear solver "magic number" for the
/// one-vs-rest `multi_class='ovr'` case, mirroring `_get_liblinear_solver_type`
/// (`sklearn/svm/_base.py:995-1049`). Returns `Err` for the unsupported
/// combinations sklearn raises `ValueError` on.
///
/// The supported (`multi_class='ovr'`) entries of `_solver_type_dict`
/// (`_base.py:1013-1014`):
///
/// ```text
///   hinge:         { l2: { dual=true: 3 } }
///   squared_hinge: { l1: { dual=false: 5 }, l2: { dual=false: 2, dual=true: 1 } }
/// ```
///
/// (`crammer_singer` = 4 is REQ-6/#623, out of scope; this function assumes the
/// existing `'ovr'` multi-class.) The `Err` strings mirror sklearn's
/// `error_string` (`_base.py:1033-1043`).
fn liblinear_solver_type(
    penalty: LinearSVCPenalty,
    loss: LinearSVCLoss,
    dual: bool,
) -> Result<u8, FerroError> {
    match (loss, penalty, dual) {
        // hinge: { l2: { dual=true: 3 } }
        (LinearSVCLoss::Hinge, LinearSVCPenalty::L2, true) => Ok(3),
        (LinearSVCLoss::Hinge, LinearSVCPenalty::L2, false) => Err(FerroError::InvalidParameter {
            name: "dual".into(),
            reason: "The combination of penalty='l2' and loss='hinge' are not \
                         supported when dual=false"
                .into(),
        }),
        // hinge + l1: penalty has no entry under `hinge` → combination unsupported.
        (LinearSVCLoss::Hinge, LinearSVCPenalty::L1, _) => Err(FerroError::InvalidParameter {
            name: "penalty".into(),
            reason: "The combination of penalty='l1' and loss='hinge' is not supported".into(),
        }),
        // squared_hinge: { l1: { dual=false: 5 } }
        (LinearSVCLoss::SquaredHinge, LinearSVCPenalty::L1, false) => Ok(5),
        (LinearSVCLoss::SquaredHinge, LinearSVCPenalty::L1, true) => {
            Err(FerroError::InvalidParameter {
                name: "dual".into(),
                reason: "The combination of penalty='l1' and loss='squared_hinge' are not \
                         supported when dual=true"
                    .into(),
            })
        }
        // squared_hinge: { l2: { dual=false: 2, dual=true: 1 } }
        (LinearSVCLoss::SquaredHinge, LinearSVCPenalty::L2, false) => Ok(2),
        (LinearSVCLoss::SquaredHinge, LinearSVCPenalty::L2, true) => Ok(1),
    }
}

/// Resolve the [`DualMode`] to a concrete `bool`, mirroring
/// `_validate_dual_parameter` (`sklearn/svm/_classes.py:13-29`).
///
/// For [`DualMode::Auto`]: when `n_samples < n_features` try `dual=true` (fall
/// back to `false` if that penalty×loss combination has no dual solver); else
/// (`n_samples >= n_features`) try `dual=false` (fall back to `true` if there is
/// no primal solver, e.g. `hinge`). Resolution is checked against
/// [`liblinear_solver_type`] so it is automatically consistent with the
/// solver-selection matrix.
fn resolve_dual(
    dual: DualMode,
    penalty: LinearSVCPenalty,
    loss: LinearSVCLoss,
    n_samples: usize,
    n_features: usize,
) -> bool {
    match dual {
        DualMode::True => true,
        DualMode::False => false,
        DualMode::Auto => {
            if n_samples < n_features {
                // Prefer dual=true; fall back to false if unsupported.
                liblinear_solver_type(penalty, loss, true).is_ok()
            } else {
                // Prefer dual=false; fall back to true if no primal solver.
                liblinear_solver_type(penalty, loss, false).is_err()
            }
        }
    }
}

/// Solve a single binary L1-regularized L2-loss (squared-hinge) SVM via
/// liblinear's feature-major coordinate descent, `solve_l1r_l2_svc`
/// (`sklearn/svm/src/liblinear/linear.cpp:1467`). Minimizes
///
/// ```text
///   ‖w‖₁  +  C · Σ_i  max(0, 1 − y_i·(w·x_i))²
/// ```
///
/// (the `l1`-penalty objective — a genuinely different, sparse optimum from the
/// l2 dual). The augmented intercept column (value `intercept_scaling`) is
/// appended when `fit_intercept`, penalized in `‖w‖₁` like any feature
/// (`coef_ = w[:n_features]`, `intercept_ = intercept_scaling·w[n_features]`).
///
/// State (`linear.cpp:1488-1526`): `b[i] = 1 − y_i·(w·x_i)` (running residual),
/// `xj_sq[j] = Σ_i C·(y_i·x_ij)²`. Per feature `j`:
/// `G_loss = −2·Σ_{i: b[i]>0} C·(y_i·x_ij)·b[i]`,
/// `H = max(2·Σ_{i: b[i]>0} C·(y_i·x_ij)², 1e-12)` (`linear.cpp:1542-1562`).
/// `Gp = G_loss+1`, `Gn = G_loss−1`. Newton direction with soft-threshold
/// (`linear.cpp:1589-1595`): `d = −Gp/H` if `Gp < H·w[j]`, `d = −Gn/H` if
/// `Gn > H·w[j]`, else `d = −w[j]`. Then a backtracking line search
/// (`sigma=0.01`, ≤20 steps, halving `d`, `linear.cpp:1600-1661`) updating
/// `b[]`, then `w[j] += d`. Shrinking via `active_size`/`Gmax_old`
/// (`linear.cpp:1567-1579, 1691-1705`); stop when
/// `Gnorm1_new ≤ eps·Gnorm1_init` on the full active set (`linear.cpp:1691`).
///
/// `C[i] = (y_i > 0 ? Cp : Cn)` per-sample (`class_weight` scales `C` per class,
/// `linear.cpp:1504-1509`; `Cp = C·w[+]`, `Cn = C·w[-]` binary / base `C`
/// multiclass, `:2543-2571`). liblinear shuffles
/// `index` each sweep (`bounded_rand_int`, `linear.cpp:1535`); ferrolearn sweeps
/// NATURAL ORDER for determinism (no RNG) — the l1 optimum is unique so the
/// limit is identical (the documented RNG-path boundary, as `solve_binary_dual`).
/// We use `eps = tol` directly: liblinear scales `primal_solver_tol`
/// (`linear.cpp:2321,2374`) but the unique optimum is `tol`-invariant at the
/// limit (the test drives `tol=1e-10` + huge `max_iter`).
///
/// Returns `(w_augmented, n_iter, converged)`.
#[allow(
    clippy::too_many_lines,
    reason = "faithful transcription of liblinear solve_l1r_l2_svc (linear.cpp:1467)"
)]
fn solve_binary_l1r_l2<F: Float + 'static>(
    x: &Array2<F>,
    y_signed: &Array1<F>,
    cfg: &SolverConfig<F>,
) -> (Vec<F>, usize, bool) {
    let SolverConfig {
        cp,
        cn,
        max_iter,
        tol,
        fit_intercept,
        intercept_scaling,
        ..
    } = *cfg;

    let (n_samples, n_features) = x.dim();

    // Per-sample penalty `C[i] = (y_i > 0 ? Cp : Cn)` (`solve_l1r_l2_svc`,
    // `linear.cpp:1504-1509`); `class_weight` makes Cp/Cn differ.
    let c_of = |i: usize| -> F { if y_signed[i] > F::zero() { cp } else { cn } };
    let w_size = if fit_intercept {
        n_features + 1
    } else {
        n_features
    };

    let inf = F::infinity();
    let two = F::one() + F::one();
    let sigma = F::from(0.01).unwrap_or_else(|| F::one() / (two * two * two * two * two * two));
    let tiny = F::from(1.0e-12).unwrap_or_else(F::epsilon);
    let max_num_linesearch = 20usize;
    let nl = F::from(n_samples).unwrap_or_else(F::one);

    // `yx(i, j)` = y_i·x_ij over the augmented design matrix (the j-th feature
    // column entry for sample i). liblinear stores `x->value *= y[ind]`
    // (`linear.cpp:1520`); we recompute it lazily for determinism / clarity.
    let yx = |i: usize, j: usize| -> F {
        let yi = y_signed[i];
        if j < n_features {
            yi * x[[i, j]]
        } else {
            yi * intercept_scaling
        }
    };

    // b[i] = 1 − y_i·(w·x_i). w starts at 0 so b starts at 1 (`linear.cpp:1500`).
    let mut b = vec![F::one(); n_samples];
    let mut w = vec![F::zero(); w_size];

    // xj_sq[j] = Σ_i C[i]·(y_i·x_ij)² (`linear.cpp:1523`, per-sample C[i]).
    let mut xj_sq = vec![F::zero(); w_size];
    for (j, xj_sq_j) in xj_sq.iter_mut().enumerate() {
        let mut acc = F::zero();
        for i in 0..n_samples {
            let val = yx(i, j);
            acc = acc + c_of(i) * val * val;
        }
        *xj_sq_j = acc;
    }

    let mut index: Vec<usize> = (0..w_size).collect();
    let mut active_size = w_size;
    let mut gmax_old = inf;
    let mut gnorm1_init = -F::one();

    let mut n_iter: usize = 0;
    let mut converged = false;

    while n_iter < max_iter {
        let mut gmax_new = F::zero();
        let mut gnorm1_new = F::zero();

        // liblinear shuffles `index[0..active_size]` here (`linear.cpp:1533-1537`);
        // ferrolearn sweeps natural order (no RNG); the unique optimum is
        // unchanged at the limit.

        let mut s = 0;
        while s < active_size {
            let j = index[s];

            // G_loss = −2·Σ_{i: b[i]>0} C·(y_i·x_ij)·b[i];
            // H = 2·Σ_{i: b[i]>0} C·(y_i·x_ij)² (`linear.cpp:1542-1561`).
            let mut g_loss = F::zero();
            let mut h = F::zero();
            for (i, &bi) in b.iter().enumerate() {
                if bi > F::zero() {
                    let val = yx(i, j);
                    let tmp = c_of(i) * val;
                    g_loss = g_loss - tmp * bi;
                    h = h + tmp * val;
                }
            }
            g_loss = g_loss * two;
            let g = g_loss;
            h = h * two;
            if h < tiny {
                h = tiny;
            }

            let gp = g + F::one();
            let gn = g - F::one();
            let wj = w[j];

            // Violation + shrinking (`linear.cpp:1564-1587`).
            let mut violation = F::zero();
            if wj == F::zero() {
                if gp < F::zero() {
                    violation = -gp;
                } else if gn > F::zero() {
                    violation = gn;
                } else if gp > gmax_old / nl && gn < -(gmax_old / nl) {
                    active_size -= 1;
                    index.swap(s, active_size);
                    continue; // re-process the swapped-in element at `s`
                }
            } else if wj > F::zero() {
                violation = gp.abs();
            } else {
                violation = gn.abs();
            }

            if violation > gmax_new {
                gmax_new = violation;
            }
            gnorm1_new = gnorm1_new + violation;

            // Newton direction with soft-threshold (`linear.cpp:1589-1595`).
            let mut d = if gp < h * wj {
                -gp / h
            } else if gn > h * wj {
                -gn / h
            } else {
                -wj
            };

            if d.abs() < tiny {
                s += 1;
                continue;
            }

            // Backtracking line search (`linear.cpp:1600-1661`).
            let mut delta = (wj + d).abs() - wj.abs() + g * d;
            let mut d_old = F::zero();
            let mut num_linesearch = 0usize;
            while num_linesearch < max_num_linesearch {
                let d_diff = d_old - d;
                let mut cond = (wj + d).abs() - wj.abs() - sigma * delta;

                let appxcond = xj_sq[j] * d * d + g_loss * d + cond;
                if appxcond <= F::zero() {
                    for (i, bi) in b.iter_mut().enumerate() {
                        *bi = *bi + d_diff * yx(i, j);
                    }
                    break;
                }

                let mut loss_old = F::zero();
                let mut loss_new = F::zero();
                if num_linesearch == 0 {
                    for (i, bi) in b.iter_mut().enumerate() {
                        if *bi > F::zero() {
                            loss_old = loss_old + c_of(i) * *bi * *bi;
                        }
                        let b_new = *bi + d_diff * yx(i, j);
                        *bi = b_new;
                        if b_new > F::zero() {
                            loss_new = loss_new + c_of(i) * b_new * b_new;
                        }
                    }
                } else {
                    for (i, bi) in b.iter_mut().enumerate() {
                        let b_new = *bi + d_diff * yx(i, j);
                        *bi = b_new;
                        if b_new > F::zero() {
                            loss_new = loss_new + c_of(i) * b_new * b_new;
                        }
                    }
                }

                cond = cond + loss_new - loss_old;
                if cond <= F::zero() {
                    break;
                }
                d_old = d;
                d = d / two;
                delta = delta / two;
                num_linesearch += 1;
            }

            w[j] = w[j] + d;

            // Recompute b[] if the line search took the maximum steps
            // (`linear.cpp:1665-1682`).
            if num_linesearch >= max_num_linesearch {
                for bi in b.iter_mut() {
                    *bi = F::one();
                }
                for (jj, &wjj) in w.iter().enumerate() {
                    if wjj == F::zero() {
                        continue;
                    }
                    for (i, bi) in b.iter_mut().enumerate() {
                        *bi = *bi - wjj * yx(i, jj);
                    }
                }
            }

            s += 1;
        }

        if n_iter == 0 {
            gnorm1_init = gnorm1_new;
        }
        n_iter += 1;

        // Stop when Gnorm1_new ≤ eps·Gnorm1_init on the full active set
        // (`linear.cpp:1691-1702`).
        if gnorm1_new <= tol * gnorm1_init {
            if active_size == w_size {
                converged = true;
                break;
            }
            active_size = w_size;
            gmax_old = inf;
            continue;
        }

        gmax_old = gmax_new;
    }

    (w, n_iter, converged)
}

/// Solve the joint Crammer-Singer multiclass SVM, transcribing liblinear's
/// `Solver_MCSVM_CS` (`sklearn/svm/src/liblinear/linear.cpp:493-787`,
/// solver type 4, `_base.py:1017`).
///
/// Minimizes the Crammer-Singer objective over a single joint weight matrix
/// `w` flattened as `(w_size × nr_class)` with `w[feature*nr_class + m]`. The
/// dual variables `alpha[i*nr_class + m]` satisfy `Σ_m alpha[i,m] = 0`,
/// `alpha[i,m] <= C[i]` if `y_i == m` else `alpha[i,m] <= 0`. Per sample,
/// `C[i] = W[i] · weighted_C[y_i]` (`linear.cpp:521-522`); here `W[i] = 1` and
/// `weighted_C[c] = C · class_weight[c]` (`= C` when `class_weight=None`).
///
/// `y_class[i]` is the class index (`0..nr_class`, the sorted-`classes`
/// position) of sample `i`. When `fit_intercept`, the design matrix is augmented
/// with the constant column `intercept_scaling` at feature index `n_features`
/// (it IS part of `QD` and `w`, `linear.cpp:512`/`608-618` over the augmented
/// row).
///
/// State / sweep / shrinking (`linear.cpp:576-754`) is transcribed faithfully,
/// including per-sample shrinking (`active_size_i`, `alpha_index`) and the
/// two-level stopping (`eps_shrink = max(10·eps, 1)`). liblinear shuffles the
/// `index` set each sweep (`bounded_rand_int`, `linear.cpp:629`); ferrolearn
/// sweeps NATURAL ORDER (no RNG) — the Crammer-Singer optimum is unique, so the
/// limit is identical (the documented RNG/shrink-path boundary, as the dual/l1
/// solvers). `eps = tol` (the solver receives `param->eps` directly,
/// `linear.cpp:2535`).
///
/// Returns `(w_flat, n_iter, converged)` where `w_flat[feature*nr_class + m]`.
#[allow(
    clippy::too_many_lines,
    clippy::too_many_arguments,
    reason = "faithful transcription of liblinear Solver_MCSVM_CS (linear.cpp:493-787); \
              the args mirror the solver's ctor + Solve() inputs (prob/nr_class/weighted_C/\
              eps/max_iter/bias) and grouping them would obscure the 1:1 transcription"
)]
fn solve_crammer_singer<F: Float + 'static>(
    x: &Array2<F>,
    y_class: &[usize],
    nr_class: usize,
    weighted_c: &[F],
    max_iter: usize,
    tol: F,
    fit_intercept: bool,
    intercept_scaling: F,
) -> (Vec<F>, usize, bool) {
    let (l, n_features) = x.dim();
    let w_size = if fit_intercept {
        n_features + 1
    } else {
        n_features
    };

    let inf = F::infinity();
    let ten = F::from(10).unwrap_or_else(|| {
        let two = F::one() + F::one();
        two + two + two + two + two
    });
    let tiny = F::from(1.0e-12).unwrap_or_else(F::epsilon);

    // `C[i] = W[i] · weighted_C[y_i]`; `W[i] = 1` (`linear.cpp:521-522`).
    let c_per_sample: Vec<F> = y_class.iter().map(|&yi| weighted_c[yi]).collect();

    // `x_val(i, feat)` over the augmented row (feature index `n_features` is the
    // constant intercept_scaling column when fit_intercept).
    let x_val = |i: usize, feat: usize| -> F {
        if feat < n_features {
            x[[i, feat]]
        } else {
            intercept_scaling
        }
    };

    // alpha[i*nr_class + m], w[feature*nr_class + m] (`linear.cpp:580-602`).
    let mut alpha = vec![F::zero(); l * nr_class];
    let mut w = vec![F::zero(); w_size * nr_class];

    // alpha_index[i*nr_class + m] = m; QD[i] = ||x_i||^2 over the augmented row;
    // active_size_i[i] = nr_class; y_index[i] = y_class[i] (`linear.cpp:603-622`).
    let mut alpha_index = vec![0usize; l * nr_class];
    let mut qd = vec![F::zero(); l];
    let mut active_size_i = vec![nr_class; l];
    let mut y_index = y_class.to_vec();
    let mut index: Vec<usize> = (0..l).collect();
    for i in 0..l {
        for m in 0..nr_class {
            alpha_index[i * nr_class + m] = m;
        }
        let mut acc = F::zero();
        for feat in 0..w_size {
            let v = x_val(i, feat);
            acc = acc + v * v;
        }
        qd[i] = acc;
    }

    // Scratch buffers reused per sample (`B`, `G`, `alpha_new`,
    // `linear.cpp:518-519,581`).
    let mut b_buf = vec![F::zero(); nr_class];
    let mut g_buf = vec![F::zero(); nr_class];
    let mut alpha_new = vec![F::zero(); nr_class];

    let mut active_size = l;
    // eps_shrink = max(10·eps, 1) (`linear.cpp:590`).
    let mut eps_shrink = (ten * tol).max(F::one());
    let mut start_from_all = true;

    let mut iter = 0usize;
    let mut converged = false;

    while iter < max_iter {
        let mut stopping = -inf;

        // liblinear shuffles index[0..active_size] here (`linear.cpp:627-631`);
        // ferrolearn sweeps natural order (no RNG); the unique optimum is
        // unchanged at the limit.

        let mut s = 0;
        while s < active_size {
            let i = index[s];
            let ai = qd[i];

            if ai > F::zero() {
                let asi = active_size_i[i];
                // G[m] = (m==y_i ? 0 : 1) + w_{alpha_index[m]} · x_i
                // (`linear.cpp:641-653`).
                for g in g_buf.iter_mut().take(asi) {
                    *g = F::one();
                }
                if y_index[i] < asi {
                    g_buf[y_index[i]] = F::zero();
                }
                for feat in 0..w_size {
                    let xv = x_val(i, feat);
                    if xv == F::zero() {
                        continue;
                    }
                    let base = feat * nr_class;
                    for m in 0..asi {
                        let idx = alpha_index[i * nr_class + m];
                        g_buf[m] = g_buf[m] + w[base + idx] * xv;
                    }
                }

                // minG over {alpha_i[idx]<0} ∪ {y_i if alpha_i[y_i]<C[i]};
                // maxG over all active m (`linear.cpp:655-666`).
                let mut min_g = inf;
                let mut max_g = -inf;
                for m in 0..asi {
                    let idx = alpha_index[i * nr_class + m];
                    if alpha[i * nr_class + idx] < F::zero() && g_buf[m] < min_g {
                        min_g = g_buf[m];
                    }
                    if g_buf[m] > max_g {
                        max_g = g_buf[m];
                    }
                }
                if y_index[i] < asi
                    && alpha[i * nr_class + y_class[i]] < c_per_sample[i]
                    && g_buf[y_index[i]] < min_g
                {
                    min_g = g_buf[y_index[i]];
                }

                // Per-sample shrinking via be_shrunk (`linear.cpp:668-697`).
                let mut m = 0;
                while m < active_size_i[i] {
                    let idx_m = alpha_index[i * nr_class + m];
                    if cs_be_shrunk(
                        c_per_sample[i],
                        m,
                        y_index[i],
                        alpha[i * nr_class + idx_m],
                        g_buf[m],
                        min_g,
                    ) {
                        active_size_i[i] -= 1;
                        while active_size_i[i] > m {
                            let asi_top = active_size_i[i];
                            let idx_top = alpha_index[i * nr_class + asi_top];
                            if !cs_be_shrunk(
                                c_per_sample[i],
                                asi_top,
                                y_index[i],
                                alpha[i * nr_class + idx_top],
                                g_buf[asi_top],
                                min_g,
                            ) {
                                alpha_index.swap(i * nr_class + m, i * nr_class + asi_top);
                                g_buf.swap(m, asi_top);
                                if y_index[i] == asi_top {
                                    y_index[i] = m;
                                } else if y_index[i] == m {
                                    y_index[i] = asi_top;
                                }
                                break;
                            }
                            active_size_i[i] -= 1;
                        }
                    }
                    m += 1;
                }

                if active_size_i[i] <= 1 {
                    active_size -= 1;
                    index.swap(s, active_size);
                    continue; // re-process the swapped-in element at `s`
                }

                if max_g - min_g <= tiny {
                    s += 1;
                    continue;
                }
                stopping = stopping.max(max_g - min_g);

                // B[m] = G[m] - Ai·alpha_i[idx] (`linear.cpp:704-705`).
                let asi = active_size_i[i];
                for m in 0..asi {
                    let idx = alpha_index[i * nr_class + m];
                    b_buf[m] = g_buf[m] - ai * alpha[i * nr_class + idx];
                }

                cs_solve_sub_problem(ai, y_index[i], c_per_sample[i], asi, &b_buf, &mut alpha_new);

                // Apply d = alpha_new[m] - alpha_i[idx] and update w
                // (`linear.cpp:708-728`).
                let mut d_nz: Vec<(usize, F)> = Vec::new();
                for m in 0..asi {
                    let idx = alpha_index[i * nr_class + m];
                    let d = alpha_new[m] - alpha[i * nr_class + idx];
                    alpha[i * nr_class + idx] = alpha_new[m];
                    if d.abs() >= tiny {
                        d_nz.push((idx, d));
                    }
                }
                if !d_nz.is_empty() {
                    for feat in 0..w_size {
                        let xv = x_val(i, feat);
                        if xv == F::zero() {
                            continue;
                        }
                        let base = feat * nr_class;
                        for &(idx, d) in &d_nz {
                            w[base + idx] = w[base + idx] + d * xv;
                        }
                    }
                }
            }

            s += 1;
        }

        iter += 1;

        // Two-level stopping (`linear.cpp:738-753`).
        if stopping < eps_shrink {
            if stopping < tol && start_from_all {
                converged = true;
                break;
            }
            active_size = l;
            for asi in active_size_i.iter_mut() {
                *asi = nr_class;
            }
            eps_shrink = (eps_shrink / (F::one() + F::one())).max(tol);
            start_from_all = true;
        } else {
            start_from_all = false;
        }
    }

    (w, iter, converged)
}

/// `Solver_MCSVM_CS::be_shrunk` (`linear.cpp:566-574`): shrink class `m` of
/// sample `i` when its dual variable is at its bound and the gradient is below
/// `minG`. `bound = C[i]` if `m == yi` (the class index), else `0`.
fn cs_be_shrunk<F: Float>(c_yi: F, m: usize, yi: usize, alpha_i: F, g_m: F, min_g: F) -> bool {
    let bound = if m == yi { c_yi } else { F::zero() };
    alpha_i == bound && g_m < min_g
}

/// `Solver_MCSVM_CS::solve_sub_problem` (`linear.cpp:541-564`): the per-sample
/// simplex-projection inner solve. `B[..active_i]` is the gradient offset
/// buffer; the result is written to `alpha_new[..active_i]`.
///
/// Clones `B → D`, adds `Ai·C_yi` to `D[yi]` (when `yi < active_i`), sorts `D`
/// DESCENDING, then finds the breakpoint `beta` and projects each coordinate:
/// `alpha_new[r] = min(C_yi, (beta-B[r])/Ai)` for `r == yi`, else
/// `min(0, (beta-B[r])/Ai)`.
fn cs_solve_sub_problem<F: Float + 'static>(
    ai: F,
    yi: usize,
    c_yi: F,
    active_i: usize,
    b_buf: &[F],
    alpha_new: &mut [F],
) {
    // clone(D, B, active_i); D[yi] += Ai·C_yi (`linear.cpp:546-548`).
    let mut d: Vec<F> = b_buf[..active_i].to_vec();
    if yi < active_i {
        d[yi] = d[yi] + ai * c_yi;
    }
    // qsort DESCENDING (compare_double returns -1 when a > b, `linear.cpp:532-538`).
    d.sort_by(|a, b| match b.partial_cmp(a) {
        Some(ord) => ord,
        None => core::cmp::Ordering::Equal,
    });

    // beta = D[0] - Ai·C_yi; for r=1; r<active_i && beta<r·D[r]; r++ { beta+=D[r]; }
    // beta /= r (`linear.cpp:551-554`).
    let mut beta = d[0] - ai * c_yi;
    let mut r = 1usize;
    while r < active_i {
        let r_f = F::from(r).unwrap_or_else(F::one);
        if beta < r_f * d[r] {
            beta = beta + d[r];
            r += 1;
        } else {
            break;
        }
    }
    let r_f = F::from(r).unwrap_or_else(F::one);
    beta = beta / r_f;

    // Project each coordinate (`linear.cpp:556-562`).
    for (rr, an) in alpha_new.iter_mut().enumerate().take(active_i) {
        let cand = (beta - b_buf[rr]) / ai;
        *an = if rr == yi {
            c_yi.min(cand)
        } else {
            F::zero().min(cand)
        };
    }
}

impl<F: Float + Send + Sync + ScalarOperand + 'static> Fit<Array2<F>, Array1<usize>>
    for LinearSVC<F>
{
    type Fitted = FittedLinearSVC<F>;
    type Error = FerroError;

    /// Fit the linear SVC model using liblinear's dual coordinate descent.
    ///
    /// Minimizes `0.5*||w||^2 + C * sum_i L(y_i, w.x_i)` (no `1/n`). When
    /// `fit_intercept` is set, the design matrix is augmented with a constant
    /// column equal to `intercept_scaling`, the augmented weight is penalized
    /// like any feature, and `intercept_ = intercept_scaling * w_last`
    /// (`_base.py:1188-1198,:1240-1245`).
    ///
    /// # Errors
    ///
    /// - [`FerroError::ShapeMismatch`] — sample count mismatch.
    /// - [`FerroError::InvalidParameter`] — `C` not positive, or
    ///   `intercept_scaling` not positive when fitting an intercept.
    /// - [`FerroError::InsufficientSamples`] — fewer than 2 distinct classes.
    fn fit(&self, x: &Array2<F>, y: &Array1<usize>) -> Result<FittedLinearSVC<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: "y length must match number of samples in X".into(),
            });
        }

        if self.c <= F::zero() {
            return Err(FerroError::InvalidParameter {
                name: "C".into(),
                reason: "must be positive".into(),
            });
        }

        // `_parameter_constraints["tol"] = Interval(Real, 0.0, None,
        // closed="neither")` (`_classes.py:237`) → `tol <= 0` raises.
        if self.tol <= F::zero() {
            return Err(FerroError::InvalidParameter {
                name: "tol".into(),
                reason: "must be positive".into(),
            });
        }

        // liblinear raises when intercept_scaling <= 0 with fit_intercept
        // (`_base.py:1190-1196`).
        if self.fit_intercept && self.intercept_scaling <= F::zero() {
            return Err(FerroError::InvalidParameter {
                name: "intercept_scaling".into(),
                reason: "must be greater than 0 when fit_intercept is true".into(),
            });
        }

        // Non-finite input validation (#2263). sklearn `LinearSVC.fit`
        // -> `self._validate_data(X, y, ...)` (`svm/_classes.py:302`, the shared
        // liblinear `_fit_liblinear` path; cf. `svm/_base.py:190`) 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 dual coordinate-descent solve. `y`
        // is `Array1<usize>` (integer class labels), finite by type, and
        // ferrolearn's `LinearSVC::fit` takes no `sample_weight` argument, so X is
        // the only runtime check. `.iter().any(|v| !v.is_finite())` rejects both
        // NaN and Inf (bounds-safe, no panic, R-CODE-2); 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(),
            });
        }

        let y_vec: Vec<usize> = y.to_vec();
        let mut classes: Vec<usize> = y_vec.clone();
        classes.sort_unstable();
        classes.dedup();

        if classes.len() < 2 {
            return Err(FerroError::InsufficientSamples {
                required: 2,
                actual: classes.len(),
                context: "LinearSVC requires at least 2 distinct classes".into(),
            });
        }

        // `class_weight` scales `C` per class: `weighted_C[i] = C·w[class i]`
        // (`compute_class_weight(class_weight, classes=classes_, y=y)`,
        // `_base.py:1179`; `linear.cpp:2496-2507`). `weights` is aligned to
        // `classes` (sorted unique = LabelEncoder order). Used by BOTH the OvR
        // path (below) and the Crammer-Singer joint solver.
        let weights = compute_class_weight(&self.class_weight, &classes, &y_vec);

        // `multi_class='crammer_singer'` selects liblinear solver type 4
        // REGARDLESS of penalty/loss/dual (`_get_liblinear_solver_type` returns
        // 4 for crammer_singer, `_base.py:1017,1020-1021`), so penalty/loss/dual
        // are IGNORED and the joint Crammer-Singer solver runs.
        if self.multi_class == MultiClass::CrammerSinger {
            return self.fit_crammer_singer(x, &classes, &y_vec, &weights);
        }

        // Resolve `dual` (auto → bool via `_validate_dual_parameter`,
        // `_classes.py:13-29`) then validate the penalty×loss×dual combination
        // against the liblinear solver matrix (`_get_liblinear_solver_type`,
        // `_base.py:995-1049`). Unsupported combinations raise (ValueError →
        // `InvalidParameter`). The resolved solver type selects the per-
        // sub-problem solver below.
        let dual = resolve_dual(self.dual, self.penalty, self.loss, n_samples, n_features);
        let _solver_type = liblinear_solver_type(self.penalty, self.loss, dual)?;

        // Solve one binary sub-problem (positive group penalty `cp`, negative
        // group penalty `cn`) and split the augmented weight vector into
        // (coef_, intercept_) per `_base.py:1240-1245`. For `penalty=l1`
        // use liblinear's L1 coordinate descent (`solve_l1r_l2_svc`, solver
        // type 5, `linear.cpp:1467`); for `penalty=l2` keep the dual CD — the
        // l2 optimum is dual-invariant (R-DEV-7), so dual=True and dual=False
        // reach the same `coef_`/`intercept_`.
        let penalty = self.penalty;
        let solve_one = |y_signed: &Array1<F>, cp: F, cn: F| -> (Array1<F>, F, usize, bool) {
            let cfg = SolverConfig {
                cp,
                cn,
                max_iter: self.max_iter,
                tol: self.tol,
                loss: self.loss,
                fit_intercept: self.fit_intercept,
                intercept_scaling: self.intercept_scaling,
            };
            let (w, n_iter, converged) = match penalty {
                LinearSVCPenalty::L1 => solve_binary_l1r_l2(x, y_signed, &cfg),
                LinearSVCPenalty::L2 => solve_binary_dual(x, y_signed, &cfg),
            };
            let coef = Array1::from_iter(w.iter().take(n_features).copied());
            let intercept = if self.fit_intercept {
                self.intercept_scaling * w[n_features]
            } else {
                F::zero()
            };
            (coef, intercept, n_iter, converged)
        };

        let mut any_unconverged = false;
        // n_iter_ = n_iter_.max() across the sub-problem fits (`_classes.py:338`).
        let mut max_n_iter: usize = 0;

        let fitted = if classes.len() == 2 {
            // Binary classification. liblinear's `prob.y` are the LabelEncoder
            // indices, mapped `y[i] = +1 if y_ind > 0 else -1`
            // (`linear.cpp:861-871`), so the positive class is `classes_[1]`
            // and `coef_` is for `classes_[1]` (`_base.py` / `_classes.py:311`).
            let y_signed: Array1<F> = y.mapv(|label| {
                if label == classes[1] {
                    F::one()
                } else {
                    -F::one()
                }
            });

            // Binary (`train_one(Cp=weighted_C[1], Cn=weighted_C[0])`,
            // `linear.cpp:2543-2551`): Cp = C·w[classes[1]] (the +1 class),
            // Cn = C·w[classes[0]] (the −1 class).
            let cp = self.c * weights[1];
            let cn = self.c * weights[0];
            let (coef, intercept, n_iter, converged) = solve_one(&y_signed, cp, cn);
            if !converged {
                any_unconverged = true;
            }
            max_n_iter = max_n_iter.max(n_iter);

            FittedLinearSVC {
                weight_vectors: vec![coef],
                intercepts: vec![intercept],
                classes,
                is_binary: true,
                n_features,
                n_iter: max_n_iter,
            }
        } else {
            // Multiclass: one-vs-rest. Each class is the positive (+1) class of
            // its own binary sub-problem.
            let mut weight_vectors = Vec::with_capacity(classes.len());
            let mut intercepts = Vec::with_capacity(classes.len());

            for (k, &cls) in classes.iter().enumerate() {
                let y_signed: Array1<F> =
                    y.mapv(|label| if label == cls { F::one() } else { -F::one() });
                // Multiclass OvR (`train_one(Cp=weighted_C[k], Cn=param->C)`,
                // `linear.cpp:2559-2571`): Cp = C·w[class k] (the +1 class),
                // Cn = C (the BASE C — the negative rest is UNWEIGHTED).
                let cp = self.c * weights[k];
                let cn = self.c;
                let (coef, intercept, n_iter, converged) = solve_one(&y_signed, cp, cn);
                if !converged {
                    any_unconverged = true;
                }
                max_n_iter = max_n_iter.max(n_iter);
                weight_vectors.push(coef);
                intercepts.push(intercept);
            }

            FittedLinearSVC {
                weight_vectors,
                intercepts,
                classes,
                is_binary: false,
                n_features,
                n_iter: max_n_iter,
            }
        };

        // liblinear warns when any sub-problem reaches `max_iter` without
        // satisfying the stopping criterion (`_base.py:1234-1238`). The crate's
        // warning channel is `eprintln!` (cf. qda.rs / lda.rs warnings).
        if any_unconverged {
            eprintln!("Liblinear failed to converge, increase the number of iterations.");
        }

        Ok(fitted)
    }
}

impl<F: Float + Send + Sync + ScalarOperand + 'static> LinearSVC<F> {
    /// Fit via the joint Crammer-Singer multiclass solver (solver type 4,
    /// `_base.py:1017`). Runs ONE joint solve over all classes
    /// ([`solve_crammer_singer`], `Solver_MCSVM_CS`, `linear.cpp:493-787`),
    /// extracts per-class `coef_`/`intercept_` from the flattened
    /// `w[feature*nr_class + m]`, and applies the binary collapse
    /// (`_classes.py:340-344`).
    ///
    /// `classes` is the sorted unique label set; `y_vec` the per-sample labels;
    /// `weights` the per-class `class_weight` multipliers aligned to `classes`.
    fn fit_crammer_singer(
        &self,
        x: &Array2<F>,
        classes: &[usize],
        y_vec: &[usize],
        weights: &[F],
    ) -> Result<FittedLinearSVC<F>, FerroError> {
        let (_n_samples, n_features) = x.dim();
        let nr_class = classes.len();

        // Map each sample's label to its class index (`0..nr_class`, the sorted
        // `classes` position — liblinear's LabelEncoder order). `y_class[i]`
        // indexes `w[..][m]` and `weighted_C[m]`.
        let y_class: Vec<usize> = y_vec
            .iter()
            .map(|label| classes.iter().position(|c| c == label).unwrap_or(0))
            .collect();

        // weighted_C[c] = C · class_weight[c] (`linear.cpp:2496-2507`,
        // `weighted_C[i] = param->C` then `*= weight`). `weights` already holds
        // the per-class multipliers (1.0 when class_weight=None).
        let weighted_c: Vec<F> = weights.iter().map(|&wc| self.c * wc).collect();

        let (w, n_iter, converged) = solve_crammer_singer(
            x,
            &y_class,
            nr_class,
            &weighted_c,
            self.max_iter,
            self.tol,
            self.fit_intercept,
            self.intercept_scaling,
        );

        if !converged {
            eprintln!("Liblinear failed to converge, increase the number of iterations.");
        }

        // Extract coef_[m][feature] = w[feature*nr_class + m] and
        // intercept_[m] = intercept_scaling · w[n_features*nr_class + m]
        // (the augmented-column weight), per `_base.py:1240-1245`.
        let mut weight_vectors: Vec<Array1<F>> = Vec::with_capacity(nr_class);
        let mut intercepts: Vec<F> = Vec::with_capacity(nr_class);
        for m in 0..nr_class {
            let mut row = Array1::<F>::zeros(n_features);
            for (feat, r) in row.iter_mut().enumerate() {
                *r = w[feat * nr_class + m];
            }
            weight_vectors.push(row);
            let intercept = if self.fit_intercept {
                self.intercept_scaling * w[n_features * nr_class + m]
            } else {
                F::zero()
            };
            intercepts.push(intercept);
        }

        // Binary collapse (`_classes.py:340-344`): with 2 classes, the joint
        // solve yields two rows; sklearn collapses them to a single weight
        // vector `coef_ = row_1 - row_0` and `intercept_ = int_1 - int_0`, so
        // the binary sign decision path is used.
        if nr_class == 2 {
            let collapsed_coef = &weight_vectors[1] - &weight_vectors[0];
            let collapsed_intercept = if self.fit_intercept {
                intercepts[1] - intercepts[0]
            } else {
                F::zero()
            };
            return Ok(FittedLinearSVC {
                weight_vectors: vec![collapsed_coef],
                intercepts: vec![collapsed_intercept],
                classes: classes.to_vec(),
                is_binary: true,
                n_features,
                n_iter,
            });
        }

        Ok(FittedLinearSVC {
            weight_vectors,
            intercepts,
            classes: classes.to_vec(),
            is_binary: false,
            n_features,
            n_iter,
        })
    }
}

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

    /// Predict class labels for the given feature matrix.
    ///
    /// Binary: `sign(X @ w + b)` mapped to class labels (`>= 0 → classes_[1]`).
    /// Multiclass: argmax of decision values across one-vs-rest classifiers.
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if the number of features
    /// does not match the fitted model.
    fn predict(&self, x: &Array2<F>) -> Result<Array1<usize>, FerroError> {
        let n_features = x.ncols();
        if n_features != self.n_features {
            return Err(FerroError::ShapeMismatch {
                expected: vec![self.n_features],
                actual: vec![n_features],
                context: "number of features must match fitted model".into(),
            });
        }

        let n_samples = x.nrows();
        let mut predictions = Array1::<usize>::zeros(n_samples);

        if self.is_binary {
            let scores = x.dot(&self.weight_vectors[0]) + self.intercepts[0];
            for i in 0..n_samples {
                predictions[i] = if scores[i] >= F::zero() {
                    self.classes[1]
                } else {
                    self.classes[0]
                };
            }
        } else {
            // Multiclass: pick class with highest decision value.
            for i in 0..n_samples {
                let mut best_class = 0;
                let mut best_score = F::neg_infinity();
                for (c, w) in self.weight_vectors.iter().enumerate() {
                    let score = x.row(i).dot(w) + self.intercepts[c];
                    if score > best_score {
                        best_score = score;
                        best_class = c;
                    }
                }
                predictions[i] = self.classes[best_class];
            }
        }

        Ok(predictions)
    }
}

impl<F: Float + Send + Sync + ScalarOperand + 'static> HasCoefficients<F> for FittedLinearSVC<F> {
    /// Returns the coefficient vector of the first (or only) binary sub-problem.
    fn coefficients(&self) -> &Array1<F> {
        &self.weight_vectors[0]
    }

    /// Returns the intercept of the first (or only) binary sub-problem.
    fn intercept(&self) -> F {
        self.intercepts[0]
    }
}

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

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

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

    #[test]
    fn test_default_constructor() {
        let m = LinearSVC::<f64>::new();
        assert_eq!(m.max_iter, 1000);
        assert!(m.c == 1.0);
        assert_eq!(m.loss, LinearSVCLoss::SquaredHinge);
        assert!(m.fit_intercept);
        assert!(m.intercept_scaling == 1.0);
    }

    #[test]
    fn test_builder_setters() {
        let m = LinearSVC::<f64>::new()
            .with_c(10.0)
            .with_max_iter(500)
            .with_tol(1e-6)
            .with_loss(LinearSVCLoss::Hinge)
            .with_fit_intercept(false)
            .with_intercept_scaling(2.0);
        assert!(m.c == 10.0);
        assert_eq!(m.max_iter, 500);
        assert_eq!(m.loss, LinearSVCLoss::Hinge);
        assert!(!m.fit_intercept);
        assert!(m.intercept_scaling == 2.0);
    }

    #[test]
    fn test_binary_classification() {
        let x = Array2::from_shape_vec(
            (8, 2),
            vec![
                1.0, 1.0, 1.0, 2.0, 2.0, 1.0, 2.0, 2.0, 8.0, 8.0, 8.0, 9.0, 9.0, 8.0, 9.0, 9.0,
            ],
        )
        .unwrap();
        let y = array![0, 0, 0, 0, 1, 1, 1, 1];

        let model = LinearSVC::<f64>::new().with_c(1.0).with_max_iter(5000);
        let fitted = model.fit(&x, &y).unwrap();
        let preds = fitted.predict(&x).unwrap();

        let correct: usize = preds.iter().zip(y.iter()).filter(|(p, a)| p == a).count();
        assert!(correct >= 6, "expected at least 6 correct, got {correct}");
    }

    #[test]
    fn test_binary_hinge_loss() {
        let x = Array2::from_shape_vec(
            (6, 2),
            vec![1.0, 1.0, 1.0, 2.0, 2.0, 1.0, 8.0, 8.0, 8.0, 9.0, 9.0, 8.0],
        )
        .unwrap();
        let y = array![0, 0, 0, 1, 1, 1];

        let model = LinearSVC::<f64>::new()
            .with_loss(LinearSVCLoss::Hinge)
            .with_max_iter(5000);
        let fitted = model.fit(&x, &y).unwrap();
        let preds = fitted.predict(&x).unwrap();

        let correct: usize = preds.iter().zip(y.iter()).filter(|(p, a)| p == a).count();
        assert!(correct >= 4, "expected at least 4 correct, got {correct}");
    }

    #[test]
    fn test_multiclass_classification() {
        let x = Array2::from_shape_vec(
            (9, 2),
            vec![
                0.0, 0.0, 0.5, 0.0, 0.0, 0.5, 10.0, 0.0, 10.5, 0.0, 10.0, 0.5, 0.0, 10.0, 0.5,
                10.0, 0.0, 10.5,
            ],
        )
        .unwrap();
        let y = array![0, 0, 0, 1, 1, 1, 2, 2, 2];

        let model = LinearSVC::<f64>::new().with_c(10.0).with_max_iter(5000);
        let fitted = model.fit(&x, &y).unwrap();

        assert_eq!(fitted.n_classes(), 3);
        assert_eq!(fitted.classes(), &[0, 1, 2]);

        let preds = fitted.predict(&x).unwrap();
        let correct: usize = preds.iter().zip(y.iter()).filter(|(p, a)| p == a).count();
        assert!(correct >= 7, "expected at least 7 correct, got {correct}");
    }

    #[test]
    fn test_binary_decision_function_is_1d() {
        let x = Array2::from_shape_vec(
            (6, 2),
            vec![1.0, 1.0, 1.0, 2.0, 2.0, 1.0, 8.0, 8.0, 8.0, 9.0, 9.0, 8.0],
        )
        .unwrap();
        let y = array![0, 0, 0, 1, 1, 1];

        let fitted = LinearSVC::<f64>::new()
            .with_max_iter(5000)
            .fit(&x, &y)
            .unwrap();
        let df = fitted.decision_function(&x).unwrap();
        // sklearn ravels the binary decision_function to (n,) (`_base.py:365`).
        let binary = df.as_binary().expect("binary decision is 1-D");
        assert_eq!(binary.len(), 6);
        assert!(df.as_multiclass().is_none());
        assert_eq!(df.n_samples(), 6);
    }

    #[test]
    fn test_multiclass_decision_function_is_2d() {
        let x = Array2::from_shape_vec(
            (9, 2),
            vec![
                0.0, 0.0, 0.5, 0.0, 0.0, 0.5, 10.0, 0.0, 10.5, 0.0, 10.0, 0.5, 0.0, 10.0, 0.5,
                10.0, 0.0, 10.5,
            ],
        )
        .unwrap();
        let y = array![0, 0, 0, 1, 1, 1, 2, 2, 2];

        let fitted = LinearSVC::<f64>::new()
            .with_c(10.0)
            .with_max_iter(5000)
            .fit(&x, &y)
            .unwrap();
        let df = fitted.decision_function(&x).unwrap();
        let scores = df.as_multiclass().expect("multiclass decision is 2-D");
        assert_eq!(scores.dim(), (9, 3));
        assert!(df.as_binary().is_none());
    }

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

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

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

        let model = LinearSVC::<f64>::new().with_c(0.0);
        assert!(model.fit(&x, &y).is_err());

        let model_neg = LinearSVC::<f64>::new().with_c(-1.0);
        assert!(model_neg.fit(&x, &y).is_err());
    }

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

        // fit_intercept=true + intercept_scaling<=0 is rejected
        // (`_base.py:1190-1196`).
        let model = LinearSVC::<f64>::new().with_intercept_scaling(0.0);
        assert!(model.fit(&x, &y).is_err());

        // But with fit_intercept=false, intercept_scaling is ignored.
        let model = LinearSVC::<f64>::new()
            .with_fit_intercept(false)
            .with_intercept_scaling(0.0);
        assert!(model.fit(&x, &y).is_ok());
    }

    #[test]
    fn test_fit_intercept_false_zero_intercept() {
        let x = Array2::from_shape_vec(
            (6, 2),
            vec![1.0, 1.0, 1.0, 2.0, 2.0, 1.0, 8.0, 8.0, 8.0, 9.0, 9.0, 8.0],
        )
        .unwrap();
        let y = array![0, 0, 0, 1, 1, 1];

        let fitted = LinearSVC::<f64>::new()
            .with_fit_intercept(false)
            .with_max_iter(5000)
            .fit(&x, &y)
            .unwrap();
        assert!(fitted.intercept() == 0.0);
    }

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

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

    #[test]
    fn test_has_coefficients() {
        let x = Array2::from_shape_vec(
            (6, 2),
            vec![1.0, 1.0, 1.0, 2.0, 2.0, 1.0, 8.0, 8.0, 8.0, 9.0, 9.0, 8.0],
        )
        .unwrap();
        let y = array![0, 0, 0, 1, 1, 1];

        let model = LinearSVC::<f64>::new().with_max_iter(5000);
        let fitted = model.fit(&x, &y).unwrap();
        assert_eq!(fitted.coefficients().len(), 2);
    }

    #[test]
    #[allow(
        clippy::assertions_on_constants,
        reason = "assert!(false) reports the unexpected-Err fit path without a gated panic!/expect"
    )]
    fn test_l1_penalty_smoke() {
        // Smoke check that the `penalty=l1` path (solve_binary_l1r_l2,
        // `linear.cpp:1467`, solver type 5, `_base.py:1014`) lands near the live
        // sklearn 1.5.2 oracle. The rigorous oracle pin is the critic's to add.
        //
        // Oracle (live sklearn 1.5.2; values per R-CHAR-3 — NEVER copied from
        // ferrolearn):
        //   python3 -c "import numpy as np; from sklearn.svm import LinearSVC; \
        //     X=np.array([[1.,1.],[1.,2.],[2.,1.],[2.,2.],[8.,8.],[8.,9.],[9.,8.],[9.,9.]]); \
        //     y=np.array([0,0,0,0,1,1,1,1]); \
        //     m=LinearSVC(penalty='l1',loss='squared_hinge',dual=False,C=1.0, \
        //       fit_intercept=True,max_iter=200000,tol=1e-10).fit(X,y); \
        //     print(m.coef_.tolist(), m.intercept_.tolist())"
        //   # coef [[0.1283185834966579, 0.12831858464059265]] int [-1.2079646017762715]
        const SK_COEF_0: f64 = 0.1283185834966579;
        const SK_COEF_1: f64 = 0.12831858464059265;
        const SK_INTERCEPT: f64 = -1.2079646017762715;

        let x = array![
            [1.0, 1.0],
            [1.0, 2.0],
            [2.0, 1.0],
            [2.0, 2.0],
            [8.0, 8.0],
            [8.0, 9.0],
            [9.0, 8.0],
            [9.0, 9.0],
        ];
        let y = array![0usize, 0, 0, 0, 1, 1, 1, 1];

        let result = LinearSVC::<f64>::new()
            .with_penalty(LinearSVCPenalty::L1)
            .with_loss(LinearSVCLoss::SquaredHinge)
            .with_dual(DualMode::False)
            .with_c(1.0)
            .with_max_iter(200_000)
            .with_tol(1e-10)
            .fit(&x, &y);

        let Ok(fitted) = result else {
            assert!(false, "l1 fit must succeed");
            return;
        };

        let coef = fitted.coefficients();
        assert!(
            (coef[0] - SK_COEF_0).abs() < 1e-2,
            "l1 coef[0] {} vs oracle {SK_COEF_0}",
            coef[0]
        );
        assert!(
            (coef[1] - SK_COEF_1).abs() < 1e-2,
            "l1 coef[1] {} vs oracle {SK_COEF_1}",
            coef[1]
        );
        assert!(
            (fitted.intercept() - SK_INTERCEPT).abs() < 1e-2,
            "l1 intercept {} vs oracle {SK_INTERCEPT}",
            fitted.intercept()
        );
    }

    #[test]
    #[allow(
        clippy::assertions_on_constants,
        reason = "assert!(false) reports the unexpected-Err fit path without a gated panic!/expect"
    )]
    fn test_class_weight_smoke() {
        // Smoke check that `class_weight` (scaling `C` per class —
        // `compute_class_weight`, `_base.py:1179`; `weighted_C[i] = C·w[i]`,
        // `linear.cpp:2496-2507`) lands near the live sklearn 1.5.2 oracle. The
        // rigorous oracle pin is the critic's to add.
        //
        // Oracle (live sklearn 1.5.2; values per R-CHAR-3 — NEVER copied from
        // ferrolearn):
        //   python3 -c "import numpy as np; from sklearn.svm import LinearSVC; \
        //     X=np.array([[1.,1.],[1.,2.],[2.,1.],[2.,2.],[1.5,1.5],[2.,1.5],[8.,8.],[9.,9.]]); \
        //     y=np.array([0,0,0,0,0,0,1,1]); \
        //     for cw in (None,'balanced',{0:1,1:5}): \
        //       m=LinearSVC(C=1.0,loss='squared_hinge',dual=True,fit_intercept=True, \
        //         max_iter=200000,tol=1e-10,class_weight=cw).fit(X,y); \
        //       print(cw, m.coef_.tolist(), m.intercept_.tolist())"
        //   # None       coef [[0.10056447415875154, 0.15957404219329038]] int [-1.263461307484969]
        //   # balanced   coef [[0.09936888940946959, 0.16666283617002833]] int [-1.2132032194327564]
        //   #            (compute_class_weight('balanced',...) = [0.6667, 2.0])
        //   # {0:1,1:5}  coef [[0.11058720549912869, 0.17164468739390437]] int [-1.2954689964246575]
        let x = array![
            [1.0, 1.0],
            [1.0, 2.0],
            [2.0, 1.0],
            [2.0, 2.0],
            [1.5, 1.5],
            [2.0, 1.5],
            [8.0, 8.0],
            [9.0, 9.0],
        ];
        let y = array![0usize, 0, 0, 0, 0, 0, 1, 1];

        // (class_weight, expected coef_, expected intercept_) from the oracle.
        let cases: [(ClassWeight<f64>, [f64; 2], f64); 3] = [
            (
                ClassWeight::None,
                [0.100_564_474_158_751_54, 0.159_574_042_193_290_38],
                -1.263_461_307_484_969,
            ),
            (
                ClassWeight::Balanced,
                [0.099_368_889_409_469_59, 0.166_662_836_170_028_33],
                -1.213_203_219_432_756_4,
            ),
            (
                ClassWeight::Explicit(vec![(0, 1.0), (1, 5.0)]),
                [0.110_587_205_499_128_69, 0.171_644_687_393_904_37],
                -1.295_468_996_424_657_5,
            ),
        ];

        for (cw, exp_coef, exp_int) in cases {
            let result = LinearSVC::<f64>::new()
                .with_loss(LinearSVCLoss::SquaredHinge)
                .with_dual(DualMode::True)
                .with_c(1.0)
                .with_class_weight(cw)
                .with_max_iter(200_000)
                .with_tol(1e-10)
                .fit(&x, &y);

            let Ok(fitted) = result else {
                assert!(false, "class_weight fit must succeed");
                return;
            };

            let coef = fitted.coefficients();
            assert!(
                (coef[0] - exp_coef[0]).abs() < 1e-2,
                "coef[0] {} vs oracle {}",
                coef[0],
                exp_coef[0]
            );
            assert!(
                (coef[1] - exp_coef[1]).abs() < 1e-2,
                "coef[1] {} vs oracle {}",
                coef[1],
                exp_coef[1]
            );
            assert!(
                (fitted.intercept() - exp_int).abs() < 1e-2,
                "intercept {} vs oracle {exp_int}",
                fitted.intercept()
            );
        }
    }

    #[test]
    fn test_unsupported_combinations_rejected() {
        let x = array![[1.0, 1.0], [2.0, 2.0], [8.0, 8.0], [9.0, 9.0]];
        let y = array![0usize, 0, 1, 1];

        // l1 + hinge is unsupported for any dual (`_base.py:1013` has no `l1`
        // under `hinge`).
        assert!(
            LinearSVC::<f64>::new()
                .with_penalty(LinearSVCPenalty::L1)
                .with_loss(LinearSVCLoss::Hinge)
                .fit(&x, &y)
                .is_err(),
            "l1 + hinge must be rejected"
        );

        // l2 + hinge + dual=false is unsupported (`hinge: {l2: {dual=true: 3}}`,
        // no dual=false entry).
        assert!(
            LinearSVC::<f64>::new()
                .with_penalty(LinearSVCPenalty::L2)
                .with_loss(LinearSVCLoss::Hinge)
                .with_dual(DualMode::False)
                .fit(&x, &y)
                .is_err(),
            "l2 + hinge + dual=false must be rejected"
        );

        // l1 + squared_hinge + dual=true is unsupported (`squared_hinge: {l1:
        // {dual=false: 5}}`, no dual=true entry).
        assert!(
            LinearSVC::<f64>::new()
                .with_penalty(LinearSVCPenalty::L1)
                .with_loss(LinearSVCLoss::SquaredHinge)
                .with_dual(DualMode::True)
                .fit(&x, &y)
                .is_err(),
            "l1 + squared_hinge + dual=true must be rejected"
        );
    }

    #[test]
    fn test_dual_auto_resolution() {
        // n_samples (4) >= n_features (2): auto prefers dual=false. hinge+l2:
        // dual=false has no solver, so auto must fall back to dual=true
        // (type 3, `_classes.py:22-27`) — the fit succeeds rather than rejecting.
        let x = array![[1.0, 1.0], [2.0, 2.0], [8.0, 8.0], [9.0, 9.0]];
        let y = array![0usize, 0, 1, 1];

        assert!(
            LinearSVC::<f64>::new()
                .with_loss(LinearSVCLoss::Hinge)
                .with_dual(DualMode::Auto)
                .fit(&x, &y)
                .is_ok(),
            "auto must fall back to dual=true for hinge+l2"
        );
    }

    #[test]
    #[allow(
        clippy::assertions_on_constants,
        reason = "assert!(false) reports the unexpected-Err fit path without a gated panic!/expect"
    )]
    fn test_crammer_singer_smoke() {
        // Smoke check that the joint Crammer-Singer solver (solve_crammer_singer,
        // Solver_MCSVM_CS, `linear.cpp:493-787`, solver type 4, `_base.py:1017`)
        // lands near the live sklearn 1.5.2 oracle for BOTH multiclass and the
        // binary collapse (`_classes.py:340-344`). The rigorous oracle pin is the
        // critic's to add.
        //
        // Oracle (live sklearn 1.5.2; values per R-CHAR-3 — NEVER copied from
        // ferrolearn):
        //   python3 -c "import numpy as np; from sklearn.svm import LinearSVC; \
        //     X=np.array([[0,0],[0.5,0.2],[0.2,0.5],[1,1],[4,4],[4.5,4.2],[4.2,4.5],[5,5],\
        //       [0,5],[0.5,5.2],[0.2,4.8],[1,6]],dtype=float); \
        //     y=np.array([0,0,0,0,1,1,1,1,2,2,2,2]); \
        //     m=LinearSVC(multi_class='crammer_singer',C=1.0,fit_intercept=True, \
        //       max_iter=200000,tol=1e-10).fit(X,y); \
        //     print(m.coef_.tolist(), m.intercept_.tolist(), m.predict(X).tolist())"
        //   # coef [[-0.06761865903618029,-0.24341085269880958],
        //   #       [0.30047599619312043,0.021705426371714007],
        //   #       [-0.23285733712058276,0.22170542636345167]]
        //   # int [0.9107847137145293,-0.622059023505724,-0.2887256901997209]
        //   # pred [0,0,0,0,1,1,1,1,2,2,2,2]
        let x = array![
            [0.0, 0.0],
            [0.5, 0.2],
            [0.2, 0.5],
            [1.0, 1.0],
            [4.0, 4.0],
            [4.5, 4.2],
            [4.2, 4.5],
            [5.0, 5.0],
            [0.0, 5.0],
            [0.5, 5.2],
            [0.2, 4.8],
            [1.0, 6.0],
        ];
        let y = array![0usize, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2];

        // Per-class oracle coef_/intercept_ (rows aligned to classes_ [0,1,2]).
        let exp_coef: [[f64; 2]; 3] = [
            [-0.067_618_659_036_180_29, -0.243_410_852_698_809_58],
            [0.300_475_996_193_120_43, 0.021_705_426_371_714_007],
            [-0.232_857_337_120_582_76, 0.221_705_426_363_451_67],
        ];
        let exp_int: [f64; 3] = [
            0.910_784_713_714_529_3,
            -0.622_059_023_505_724,
            -0.288_725_690_199_720_9,
        ];

        let Ok(fitted) = LinearSVC::<f64>::new()
            .with_multi_class(MultiClass::CrammerSinger)
            .with_c(1.0)
            .with_max_iter(200_000)
            .with_tol(1e-10)
            .fit(&x, &y)
        else {
            assert!(false, "crammer_singer multiclass fit must succeed");
            return;
        };

        assert!(!fitted.is_binary, "3-class CS must be multiclass");
        let rows = fitted.weight_vectors();
        let ints = fitted.intercepts();
        assert_eq!(rows.len(), 3);
        for m in 0..3 {
            assert!(
                (rows[m][0] - exp_coef[m][0]).abs() < 1e-2,
                "CS coef[{m}][0] {} vs oracle {}",
                rows[m][0],
                exp_coef[m][0]
            );
            assert!(
                (rows[m][1] - exp_coef[m][1]).abs() < 1e-2,
                "CS coef[{m}][1] {} vs oracle {}",
                rows[m][1],
                exp_coef[m][1]
            );
            assert!(
                (ints[m] - exp_int[m]).abs() < 1e-2,
                "CS intercept[{m}] {} vs oracle {}",
                ints[m],
                exp_int[m]
            );
        }

        // predict is all-correct (argmax over the per-class scores).
        let Ok(preds) = fitted.predict(&x) else {
            assert!(false, "CS multiclass predict must succeed");
            return;
        };
        assert_eq!(
            preds.to_vec(),
            y.to_vec(),
            "CS multiclass predict must be all-correct"
        );

        // Binary collapse (`_classes.py:340-344`): with 2 classes the joint solve
        // collapses to a SINGLE weight vector (`coef_` shape (1,2)) +
        // `intercept_ = int_1 - int_0`.
        //   python3 -c "import numpy as np; from sklearn.svm import LinearSVC; \
        //     X=np.array([[1,1],[1,2],[2,1],[2,2],[8,8],[8,9],[9,8],[9,9]],dtype=float); \
        //     y=np.array([0,0,0,0,1,1,1,1]); \
        //     m=LinearSVC(multi_class='crammer_singer',C=1.0,fit_intercept=True, \
        //       max_iter=200000,tol=1e-10).fit(X,y); \
        //     print(m.coef_.tolist(), m.intercept_.tolist())"
        //   # coef [[0.15503875968992287,0.15503875968992295]] int [-1.4806201550387597]
        const BIN_COEF_0: f64 = 0.155_038_759_689_922_87;
        const BIN_COEF_1: f64 = 0.155_038_759_689_922_95;
        const BIN_INT: f64 = -1.480_620_155_038_759_7;

        let xb = array![
            [1.0, 1.0],
            [1.0, 2.0],
            [2.0, 1.0],
            [2.0, 2.0],
            [8.0, 8.0],
            [8.0, 9.0],
            [9.0, 8.0],
            [9.0, 9.0],
        ];
        let yb = array![0usize, 0, 0, 0, 1, 1, 1, 1];

        let Ok(fitted_b) = LinearSVC::<f64>::new()
            .with_multi_class(MultiClass::CrammerSinger)
            .with_c(1.0)
            .with_max_iter(200_000)
            .with_tol(1e-10)
            .fit(&xb, &yb)
        else {
            assert!(false, "crammer_singer binary fit must succeed");
            return;
        };

        assert!(fitted_b.is_binary, "2-class CS must collapse to binary");
        let coef = fitted_b.coefficients();
        assert_eq!(coef.len(), 2, "binary collapse → single (1,2) weight row");
        assert!(
            (coef[0] - BIN_COEF_0).abs() < 1e-2,
            "CS binary coef[0] {} vs oracle {BIN_COEF_0}",
            coef[0]
        );
        assert!(
            (coef[1] - BIN_COEF_1).abs() < 1e-2,
            "CS binary coef[1] {} vs oracle {BIN_COEF_1}",
            coef[1]
        );
        assert!(
            (fitted_b.intercept() - BIN_INT).abs() < 1e-2,
            "CS binary intercept {} vs oracle {BIN_INT}",
            fitted_b.intercept()
        );

        let Ok(preds_b) = fitted_b.predict(&xb) else {
            assert!(false, "CS binary predict must succeed");
            return;
        };
        assert_eq!(
            preds_b.to_vec(),
            yb.to_vec(),
            "CS binary predict all-correct"
        );
    }

    #[test]
    fn test_predict_feature_mismatch() {
        let x = Array2::from_shape_vec(
            (6, 2),
            vec![1.0, 1.0, 1.0, 2.0, 2.0, 1.0, 8.0, 8.0, 8.0, 9.0, 9.0, 8.0],
        )
        .unwrap();
        let y = array![0, 0, 0, 1, 1, 1];

        let fitted = LinearSVC::<f64>::new()
            .with_max_iter(5000)
            .fit(&x, &y)
            .unwrap();

        let x_bad = Array2::from_shape_vec((3, 1), vec![1.0, 2.0, 3.0]).unwrap();
        assert!(fitted.predict(&x_bad).is_err());
    }
}