ferrolearn_linear/

svm.rs

//! Support Vector Machine with kernel trick.
//!
//! This module provides [`SVC`] (classification) and [`SVR`] (regression)
//! support vector machines trained using the **Sequential Minimal Optimization
//! (SMO)** algorithm (Platt, 1998).
//!
//! # Kernels
//!
//! Four built-in kernels are provided:
//!
//! - [`LinearKernel`]: `K(x, y) = x . y`
//! - [`RbfKernel`]: `K(x, y) = exp(-gamma * ||x - y||^2)`
//! - [`PolynomialKernel`]: `K(x, y) = (gamma * x . y + coef0)^degree`
//! - [`SigmoidKernel`]: `K(x, y) = tanh(gamma * x . y + coef0)`
//!
//! Users can implement the [`Kernel`] trait for custom kernels.
//!
//! # Multiclass
//!
//! `SVC` uses a one-vs-one strategy for multiclass classification.
//!
//! # Examples
//!
//! ```
//! use ferrolearn_linear::svm::{SVC, LinearKernel};
//! use ferrolearn_core::{Fit, Predict};
//! use ndarray::{array, Array2};
//!
//! let x = Array2::from_shape_vec((6, 2), vec![
//!     1.0, 1.0,  2.0, 1.0,  1.0, 2.0,
//!     5.0, 5.0,  6.0, 5.0,  5.0, 6.0,
//! ]).unwrap();
//! let y = array![0usize, 0, 0, 1, 1, 1];
//!
//! let model = SVC::<f64, LinearKernel>::new(LinearKernel);
//! 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 `#`. `SVC`/`SVR`/
//! `FittedSVC`/`FittedSVR`/`Kernel` + the four kernels are boundary estimator
//! types re-exported at the crate root (`pub use svm::{…}` in `lib.rs`) and
//! consumed by `nu_svm.rs` (`NuSVC`/`NuSVR` delegate to `SVC`/`SVR`) and
//! `one_class_svm.rs` (uses `Kernel`) — non-test production consumers; under
//! S5/R-DEFER-1 the fitted-attribute accessors are part of that boundary public
//! API surface. See `.design/linear/svm.md`.
//!
//! | REQ | Status | Evidence |
//! |---|---|---|
//! | REQ-1 (kernels + gamma scale/auto/float) | SHIPPED | The four kernel formulas + the three-way `pub enum Gamma<F> { Scale, Auto, Value }` (default `Scale`) resolved at fit time by `fn resolve_gamma in svm.rs` + `fn resolved_for_fit in svm.rs`: `Scale`=`1/(n_features·X.var())` (`_base.py:238-239`), `Auto`=`1/n_features` (`_base.py:240-241`), `Value(v)`=verbatim (`_base.py:242-243`). Builders `RbfKernel::with_gamma`/`with_gamma_scale`/`with_gamma_auto`. Non-test consumer: the kernel `gamma` field is resolved in the production `fn fit in svm.rs` (`self.kernel.resolved_for_fit(x)`). Pinned: `divergence_pin2_rbf_default_scale_gamma` (scale, green) + in-module `test_svc_gamma_auto_decision_function in svm.rs` (`_gamma=0.5`, df `[-0.9996,-0.9999,-0.9999,0.9999,0.9999,0.9996]` vs live `SVC(kernel='rbf',gamma='auto')`, R-CHAR-3, 1e-2) + `test_svc_gamma_scale_still_default` (`_gamma=0.118421`). |
//! | REQ-2 (C-SVC SMO fit) | SHIPPED | `fn smo_binary in svm.rs` (Fan-Chen-Lin WSS) converges to libsvm's `α`; pinned by `divergence_pin5_binary_fitted_attributes in tests/divergence_svm_fit.rs` (`dual_coef_ [[-0.0408,-0.0408,0.0816]]`, `support_ [1,2,3]`, `intercept_ [-1.8565]` vs live `SVC(kernel='linear',C=1.0)`). |
//! | REQ-3 (fitted attrs + binary sign flip) | SHIPPED | `FittedSVC::{support,support_vectors,n_support,dual_coef,intercept,coef} in svm.rs` emit the libsvm layout with the binary sign flip (`_base.py:258-262`); `coef_` is linear-only (`_base.py:642-666`). Pinned by `divergence_pin5_*` (binary) + `divergence_pin6_multiclass_dual_coef_packing` (multiclass `(n_class-1,n_SV)` packing). |
//! | REQ-4 (decision_function shape/sign/ovr) | SHIPPED | `FittedSVC::decision_function in svm.rs` returns the `SvmScores<F>` enum: binary -> `SvmScores::Binary` 1-D `(n,)` = `-raw_ovo.ravel()` (positive -> `classes_[1]`, `_base.py:538-539`); multiclass -> `SvmScores::Multiclass` `(n, n_classes)` via `fn ovr_decision_function in svm.rs` (default `SvmDecisionShape::Ovr`, transcribed from `multiclass.py:520-562`) applied to `dec<0`/`-dec` (`_base.py:780`), or raw `(n, n·(n-1)/2)` for `SvmDecisionShape::Ovo`. `SVC::decision_function_shape` field + `with_decision_function_shape`. Sign normalized: `fn raw_ovo` negates `decision_value_binary` to restore libsvm's lower-index-class-`+1` ovo convention. Pinned by `divergence_pin8_multiclass_ovr_decision_function` (ovr `(9,3)` row0 `[2.2366,0.8167,-0.1833]`, row3 `[1.0606,2.2262,-0.2333]`), `divergence_pin9_multiclass_ovo_decision_function` (ovo `(9,3)` row0 `[1.2222,1.2222,0.0]`), `divergence_pin10_binary_shape_contract` (binary 1-D `(6,)`) in `tests/divergence_svm_fit.rs` (R-CHAR-3, 1e-2). Consumer: `FittedNuSVC::decision_function in nu_svm.rs` delegates (non-test, propagates `SvmScores`). |
//! | REQ-5 (predict + tie-break) | SHIPPED | `fn predict in svm.rs` (FittedSVC) does libsvm ovo voting and breaks vote ties toward the LOWER class index via a strictly-greater first-max scan (keeps the first/lowest-index maximum since `classes` is `np.unique(y)`-sorted), matching libsvm/sklearn `super().predict` (`_base.py:813-814`) instead of `max_by_key`'s last-maximum. Pinned by `divergence_pin3_predict_labels` (separable-set labels) + `divergence_pin11_ovo_vote_tie_break_lower_index` (4-class vote tie `(0,2,2,2)` at `q=(-0.21,-8.976)` -> class 1) in `tests/divergence_svm_fit.rs` vs live `SVC(kernel='linear',C=1.0)`. |
//! | REQ-6 (epsilon-SVR) | SHIPPED | `fn smo_svr in svm.rs` + `FittedSVR::{support,support_vectors,n_support,dual_coef,intercept}`; pinned by `divergence_pin4_svr_predict_values` (predict) + `divergence_pin7_svr_fitted_attributes` (`support_ [0,5]`, `dual_coef_ [[-0.392,0.392]]`, `intercept_ [0.14]` vs live `SVR(kernel='linear',C=100,epsilon=0.1)`). |
//! | REQ-7 (multiclass one-vs-one) | SHIPPED | `fn fit in svm.rs` (SVC) trains one `smo_binary` per class pair, `classes` = `np.unique(y)`; pinned by `divergence_pin6_multiclass_dual_coef_packing` (3-class `dual_coef_ (2,6)` libsvm packing, `support_ [1,2,3,5,6,7]`, `n_support_ [2,2,2]`, `intercept_ [1.2222,1.2222,0.0]`). |
//! | REQ-8 (constructor param surface + defaults) | SHIPPED | `shrinking` (`SVC`/`SVR`, default `true`, `with_shrinking`; accepted for API parity, shrinking-invariant optimum so DOES NOT alter results — R-DEV-7); `break_ties` (`SVC`, default `false`, `with_break_ties`; `fn predict in svm.rs` ovr-argmax branch for `break_ties=true`+ovr+`n_classes>2`, `InvalidParameter` for the ovo combo, `_base.py:801-814`); default alignment `cache_size=200`, `max_iter=0` (= sklearn `-1`, no iteration limit; the `smo_binary`/`smo_svr` loops treat `0` as unbounded); REQ-1's `gamma` enum (`scale`/`auto`/float); and now **`class_weight`** (`SVC`, `pub class_weight: ClassWeight<F>` default `None`, `with_class_weight`). `fn compute_class_weight in svm.rs` mirrors `sklearn.utils.compute_class_weight` as called by `BaseSVC._validate_targets` (`class_weight_ = compute_class_weight(class_weight, classes, y)`, `_base.py:740`): `None`→1.0, `Balanced`→`n_samples/(n_classes·count_c)` (`_classes.py:122-124`), `Explicit`→1.0 default overridden by map. `fn smo_binary in svm.rs` now takes per-class box bounds `(cp, cn)` (the `y=+1`/`y=-1` upper bounds) instead of a scalar `c`, applied in the WSS `in_up`/`in_low` tests, the analytic-update box clip, and the free-SV bias recovery (`0<alpha_i<C_i`); when `cp==cn` the math is identical to before (the 13 divergence pins stay green). `fn fit in svm.rs` (SVC) computes `weights = compute_class_weight(...)` ONCE over the full `y`, then per ovo pair `(ci,cj)`: `cp = C·weights[cj]`, `cn = C·weights[ci]` (libsvm `weighted_C`, `_base.py:740`). Non-test consumer: `fn fit in svm.rs` consumes `self.class_weight` (the boundary `SVC`/`FittedSVC` types are re-exported at the crate root + consumed by `nu_svm.rs`). Pinned: `test_svc_class_weight_smoke`/`test_compute_class_weight_balanced`/`test_svc_break_ties_changes_label`/`test_svc_break_ties_ovo_errors`/`test_svc_default_params in svm.rs` (live oracle on the imbalanced 8×2 set: None `dual_coef_ [[-0.5,-1,1,0.5]]`/`intercept_ [-2.0]`/`support_ [1,3,5,6]`; balanced `[[-0.8,-0.8,1.3333,0.2667]]`/`-1.6667`; `{0:1,1:5}` `support_ [1,3,4,5]`/`-2.0`; R-CHAR-3, 1e-2; None≠balanced intercept). **R-DEV-7 design difference (preserved contract, NOT a gap):** estimator-level `kernel`(string-select)/`degree`/`coef0` are the type parameter `K`, set by construction; `random_state` is unused (ferrolearn's SMO is deterministic). `class_weight` is SVC-only (sklearn SVR has no `class_weight`). |
//! | REQ-9 (probability / predict_proba) | SHIPPED | `pub probability: bool` field on `SVC` (default `false`, `with_probability`) + `prob_a`/`prob_b`/`probability` on `FittedSVC`. The DETERMINISTIC Platt machinery is transcribed from libsvm: `fn sigmoid_train in svm.rs` (Newton iteration + prior init + target smoothing + step-halving line search, `svm.cpp:1919-2030`), `fn sigmoid_predict in svm.rs` (overflow-safe form, `svm.cpp:2032-2040`), `fn multiclass_probability in svm.rs` (Wu-Lin-Weng 2004 coupling, `svm.cpp:2043-2104`). `fn platt_cv_sigmoid in svm.rs` runs a per-ovo-pair 5-fold CV at fit time when `probability=true` (`svm.cpp:2107-2203`). `FittedSVC::predict_proba in svm.rs` builds the pairwise matrix via `sigmoid_predict` (clamped `[1e-7,1-1e-7]`, `svm.cpp:2937`) -> `multiclass_probability` -> `(n,n_classes)`; binary -> `[P(classes[0]),P(classes[1])]`; rows sum to 1. `FittedSVC::predict_log_proba` = `predict_proba.ln()` (`_base.py:866-894`). `probability=false` -> `InvalidParameter` carrying sklearn's `NotFittedError` text "predict_proba is not available when fitted with probability=False" (`_base.py:856-860`; no `NotFitted` variant by R-DEV-4 typestate). **RNG-CV boundary (documented divergence, NOT a gap):** libsvm's CV fold permutation is RNG-seeded, so sklearn's `probA_`/`probB_`/`predict_proba` are NON-DETERMINISTIC across `random_state` (`probA_` = -0.7749 at rs=0 vs -1.0541 at rs=1; the docstring admits CV-dependence). ferrolearn uses a DETERMINISTIC contiguous 5-fold split (analogous to the documented SGD shuffle boundary), so it does NOT bit-match sklearn's predict_proba VALUES — only the deterministic machinery + structural invariants + the raise contract are verified (R-CHAR-3: the asserted invariants are sklearn's DOCUMENTED contract, not copied values). Pinned by `test_svc_predict_proba_raises_when_probability_false`/`test_svc_predict_proba_binary_rows_sum_to_one`/`test_svc_predict_proba_binary_monotone_in_decision`/`test_svc_predict_log_proba_equals_log_of_proba`/`test_svc_predict_proba_multiclass_rows_sum_to_one`/`test_sigmoid_predict_overflow_safe`/`test_multiclass_probability_binary_reduces_to_pairwise in svm.rs`. Non-test consumer: `fn fit in svm.rs` (SVC) consumes `self.probability` (the boundary `SVC`/`FittedSVC` types are re-exported at the crate root + consumed by `nu_svm.rs`). |
//! | REQ-10 (ferray substrate) | NOT-STARTED | open #643. `svm.rs` imports `ndarray::{Array1, Array2, ScalarOperand}`, not `ferray-core`/`ferray::linalg` (R-SUBSTRATE). |
//! | REQ-11 (non-finite input rejected) | SHIPPED | Both fit entries reject any NaN/+/-inf BEFORE the SMO solve with `FerroError::InvalidParameter`, mirroring sklearn's `BaseLibSVM.fit` -> `_validate_data(X, y, …)` (`_base.py:190-197`, default `force_all_finite=True`) -> `ValueError("Input X contains NaN.")` / `"Input y contains NaN."` / `"… contains infinity …"`. `SVC::fit in svm.rs` checks `X` (`y` is `Array1<usize>` labels, finite by type); `SVR::fit in svm.rs` checks `X` AND the float target `y`. ferrolearn's `Fit::fit` signature has no `sample_weight` argument, so the sklearn `sample_weight`-finiteness raise has no fit-entry counterpart here. `.iter().any(|v| !v.is_finite())` catches both NaN and Inf; the finite path is byte-identical (the guard never fires on finite input — the 13+ SVC/SVR divergence pins stay green). Verified vs the live sklearn 1.5.2 oracle (R-CHAR-3): NaN/+inf/-inf in X for both, NaN/inf in y for SVR, all raise `ValueError` (`tests/divergence_svm_nonfinite.rs::{svc_*,svr_*}`). Non-test consumer: the existing `Fit::fit` consumers + the crate-root `pub use svm::{SVC, SVR, …}` re-exports. (#2269) |

use std::collections::HashMap;

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

// ---------------------------------------------------------------------------
// Kernel trait and built-in kernels
// ---------------------------------------------------------------------------

/// The `gamma` coefficient for the RBF / polynomial / sigmoid kernels,
/// mirroring scikit-learn's three-way `gamma` parameter
/// (`sklearn/svm/_base.py:235-243`,
/// `StrOptions({"scale", "auto"}) | Interval(Real, 0.0, None)`).
///
/// Resolved at fit time against the training matrix `X`
/// ([`Kernel::resolved_for_fit`]):
///
/// - [`Gamma::Scale`] (default): `1 / (n_features · X.var())` where `X.var()`
///   is the population variance (ddof=0) of the whole flattened `X`
///   (`_base.py:238-239`). When `X.var() == 0` sklearn falls back to `1.0`.
/// - [`Gamma::Auto`]: `1 / n_features` (`_base.py:240-241`).
/// - [`Gamma::Value`]: the float verbatim (`_base.py:242-243`).
///
/// The default is [`Gamma::Scale`], matching sklearn's `gamma='scale'`.
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum Gamma<F> {
    /// `gamma='scale'` (sklearn default): `1 / (n_features · X.var())`.
    Scale,
    /// `gamma='auto'`: `1 / n_features`.
    Auto,
    /// An explicit float gamma, used verbatim.
    Value(F),
}

impl<F> Default for Gamma<F> {
    /// sklearn's default is `gamma='scale'`.
    fn default() -> Self {
        Gamma::Scale
    }
}

/// Per-class scaling of the regularization parameter `C` for [`SVC`].
///
/// Mirrors `sklearn.svm.SVC`'s `class_weight` parameter
/// (`sklearn/svm/_classes.py:118-124`, constraint `{None, dict, 'balanced'}`):
/// it sets the C of class `i` to `class_weight[i]·C` (libsvm's per-class
/// `weighted_C[i] = C·class_weight_[i]`). The expanded per-class weights are
/// computed by [`compute_class_weight`] following
/// `sklearn.utils.compute_class_weight` semantics, as called from
/// `BaseSVC._validate_targets`
/// (`self.class_weight_ = compute_class_weight(self.class_weight, classes=cls,
/// y=y_)`, `sklearn/svm/_base.py:740`).
///
/// This mirrors [`crate::linear_svc::ClassWeight`] for cross-estimator
/// consistency, but is defined locally (no cross-import of `linear_svc`
/// internals).
#[derive(Debug, Clone, Default)]
pub enum ClassWeight<F> {
    /// Uniform weights (all classes weighted `1.0`). The default
    /// (`class_weight=None`).
    #[default]
    None,
    /// Balanced weights `n_samples / (n_classes · count_c)` per class `c`,
    /// matching `sklearn.utils.compute_class_weight("balanced", ...)`
    /// (`_classes.py:122-124`: `n_samples / (n_classes * np.bincount(y))`).
    Balanced,
    /// Explicit class-label -> weight map. Classes absent from the map default
    /// to `1.0`, matching the dict branch of `compute_class_weight`.
    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`, as called by
/// `BaseSVC._validate_targets`
/// (`compute_class_weight(self.class_weight, classes=cls, y=y_)`,
/// `sklearn/svm/_base.py:740`):
/// - `None` -> all `1.0`.
/// - `Balanced` -> `n_samples / (n_classes · count_c)` per class `c`,
///   where `count_c` is the number of samples with label `c`
///   (`_classes.py:122-124`).
/// - `Explicit(map)` -> `1.0` default, overridden by the map entries matched by
///   class label.
///
/// `classes` is the sorted unique label set; `y` is the per-sample label array.
/// Mirrors `ferrolearn_linear::linear_svc::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))`, indexed per
            // class (`_classes.py:124`).
            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(),
    }
}

/// A kernel function for SVM.
///
/// Computes the inner product of two vectors in a (possibly implicit)
/// higher-dimensional feature space.
pub trait Kernel<F: Float>: Clone + Send + Sync {
    /// Compute the kernel value between two vectors.
    fn compute(&self, x: &[F], y: &[F]) -> F;

    /// Resolve any data-dependent kernel parameters against the training data
    /// at fit time, returning a copy of the kernel with those parameters fixed.
    ///
    /// For kernels with a [`Gamma<F>`] parameter, the three-way `gamma`
    /// resolution mirrors scikit-learn (`sklearn/svm/_base.py:235-243`):
    /// [`Gamma::Scale`] (default) -> `1 / (n_features * X.var())` where
    /// `X.var()` is the population variance (ddof=0) over the whole flattened
    /// training matrix; [`Gamma::Auto`] -> `1 / n_features`; [`Gamma::Value`]
    /// is left verbatim. After resolution the stored `gamma` is always a
    /// concrete [`Gamma::Value`].
    ///
    /// The default implementation is a no-op (returns `self.clone()`), which is
    /// correct for parameter-free kernels such as [`LinearKernel`].
    #[must_use]
    fn resolved_for_fit(&self, _x: &Array2<F>) -> Self
    where
        Self: Sized,
    {
        self.clone()
    }

    /// Whether this is the linear kernel `K(x, y) = x . y`.
    ///
    /// sklearn exposes `coef_` (the primal weight vector
    /// `dual_coef_ @ support_vectors_`) ONLY for the linear kernel and raises
    /// `AttributeError` otherwise (`sklearn/svm/_base.py:650-651`). The default
    /// is `false`; [`LinearKernel`] overrides it to `true`.
    #[must_use]
    fn is_linear(&self) -> bool {
        false
    }
}

/// Compute the population variance (ddof=0) of all elements of `x`, mirroring
/// numpy's `X.var()` (`mean((x - mean)^2)`). Returns `None` when `x` is empty.
fn population_variance<F: Float>(x: &Array2<F>) -> Option<F> {
    let n = x.len();
    if n == 0 {
        return None;
    }
    let count = F::from(n)?;
    let sum = x.iter().fold(F::zero(), |acc, &v| acc + v);
    let mean = sum / count;
    let sq = x
        .iter()
        .fold(F::zero(), |acc, &v| acc + (v - mean) * (v - mean));
    Some(sq / count)
}

/// Extract the concrete float from a [`Gamma<F>`] for a direct `compute` call
/// without training data. After [`Kernel::resolved_for_fit`] the gamma is
/// always a [`Gamma::Value`], so this is the live path; an unresolved
/// `Scale`/`Auto` (e.g. a kernel used standalone outside a fit) has no `X` to
/// resolve against and falls back to `1.0`, matching the prior default-gamma
/// behavior of a directly-evaluated kernel.
fn gamma_value_or_one<F: Float>(gamma: Gamma<F>) -> F {
    match gamma {
        Gamma::Value(v) => v,
        Gamma::Scale | Gamma::Auto => F::one(),
    }
}

/// Resolve a [`Gamma<F>`] spec against the training matrix `X`, returning the
/// concrete float gamma, mirroring scikit-learn (`sklearn/svm/_base.py:235-243`):
///
/// - [`Gamma::Scale`] -> `1 / (n_features * X.var())` (`_base.py:238-239`).
///   When `X.var() == 0` (constant `X`) or `X` is empty, sklearn falls back to
///   `1.0` (`_base.py:239`: `if X_var != 0 else 1.0`), so we do the same
///   (avoiding a non-finite gamma).
/// - [`Gamma::Auto`] -> `1 / n_features` (`_base.py:240-241`).
/// - [`Gamma::Value`] -> the float verbatim (`_base.py:242-243`).
fn resolve_gamma<F: Float>(gamma: Gamma<F>, x: &Array2<F>) -> F {
    match gamma {
        Gamma::Value(v) => v,
        Gamma::Auto => match F::from(x.ncols()) {
            Some(nf) if nf > F::zero() => F::one() / nf,
            _ => F::one(),
        },
        Gamma::Scale => {
            let n_features = match F::from(x.ncols()) {
                Some(nf) if nf > F::zero() => nf,
                _ => return F::one(),
            };
            match population_variance(x) {
                Some(var) if var > F::zero() => F::one() / (n_features * var),
                // var == 0 (constant X) or empty: sklearn falls back to 1.0.
                _ => F::one(),
            }
        }
    }
}

/// Linear kernel: `K(x, y) = x . y`.
#[derive(Debug, Clone, Copy)]
pub struct LinearKernel;

impl<F: Float> Kernel<F> for LinearKernel {
    fn compute(&self, x: &[F], y: &[F]) -> F {
        x.iter()
            .zip(y.iter())
            .fold(F::zero(), |acc, (&a, &b)| acc + a * b)
    }

    fn is_linear(&self) -> bool {
        true
    }
}

/// Radial Basis Function (Gaussian) kernel.
///
/// `K(x, y) = exp(-gamma * ||x - y||^2)`
#[derive(Debug, Clone, Copy)]
pub struct RbfKernel<F> {
    /// The gamma parameter, a three-way [`Gamma<F>`] spec resolved at fit time
    /// (`sklearn/svm/_base.py:235-243`). Default [`Gamma::Scale`]
    /// (= `1 / (n_features * X.var())`); [`Gamma::Auto`] = `1 / n_features`;
    /// [`Gamma::Value`] is used verbatim.
    pub gamma: Gamma<F>,
}

impl<F: Float> RbfKernel<F> {
    /// Create a new RBF kernel with the default `gamma='scale'`.
    #[must_use]
    pub fn new() -> Self {
        Self {
            gamma: Gamma::Scale,
        }
    }

    /// Create a new RBF kernel with an explicit float gamma
    /// (`gamma=<float>`, [`Gamma::Value`]).
    #[must_use]
    pub fn with_gamma(gamma: F) -> Self {
        Self {
            gamma: Gamma::Value(gamma),
        }
    }

    /// Create a new RBF kernel with `gamma='scale'` ([`Gamma::Scale`],
    /// sklearn's default = `1 / (n_features * X.var())`).
    #[must_use]
    pub fn with_gamma_scale() -> Self {
        Self {
            gamma: Gamma::Scale,
        }
    }

    /// Create a new RBF kernel with `gamma='auto'` ([`Gamma::Auto`]
    /// = `1 / n_features`).
    #[must_use]
    pub fn with_gamma_auto() -> Self {
        Self { gamma: Gamma::Auto }
    }
}

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

impl<F: Float + Send + Sync> Kernel<F> for RbfKernel<F> {
    fn compute(&self, x: &[F], y: &[F]) -> F {
        let gamma = gamma_value_or_one(self.gamma);
        let sq_dist = x.iter().zip(y.iter()).fold(F::zero(), |acc, (&a, &b)| {
            let d = a - b;
            acc + d * d
        });
        (-gamma * sq_dist).exp()
    }

    fn resolved_for_fit(&self, x: &Array2<F>) -> Self {
        Self {
            gamma: Gamma::Value(resolve_gamma(self.gamma, x)),
        }
    }
}

/// Polynomial kernel: `K(x, y) = (gamma * x . y + coef0)^degree`.
#[derive(Debug, Clone, Copy)]
pub struct PolynomialKernel<F> {
    /// The gamma parameter, a three-way [`Gamma<F>`] spec resolved at fit time
    /// (`sklearn/svm/_base.py:235-243`). Default [`Gamma::Scale`].
    pub gamma: Gamma<F>,
    /// Polynomial degree.
    pub degree: usize,
    /// Independent term.
    pub coef0: F,
}

impl<F: Float> PolynomialKernel<F> {
    /// Create a new polynomial kernel with defaults (`gamma='scale'`,
    /// `degree=3`, `coef0=0`).
    #[must_use]
    pub fn new() -> Self {
        Self {
            gamma: Gamma::Scale,
            degree: 3,
            coef0: F::zero(),
        }
    }
}

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

impl<F: Float + Send + Sync> Kernel<F> for PolynomialKernel<F> {
    fn compute(&self, x: &[F], y: &[F]) -> F {
        let gamma = gamma_value_or_one(self.gamma);
        let dot: F = x
            .iter()
            .zip(y.iter())
            .fold(F::zero(), |acc, (&a, &b)| acc + a * b);
        let val = gamma * dot + self.coef0;
        let mut result = F::one();
        for _ in 0..self.degree {
            result = result * val;
        }
        result
    }

    fn resolved_for_fit(&self, x: &Array2<F>) -> Self {
        Self {
            gamma: Gamma::Value(resolve_gamma(self.gamma, x)),
            degree: self.degree,
            coef0: self.coef0,
        }
    }
}

/// Sigmoid kernel: `K(x, y) = tanh(gamma * x . y + coef0)`.
#[derive(Debug, Clone, Copy)]
pub struct SigmoidKernel<F> {
    /// The gamma parameter, a three-way [`Gamma<F>`] spec resolved at fit time
    /// (`sklearn/svm/_base.py:235-243`). Default [`Gamma::Scale`].
    pub gamma: Gamma<F>,
    /// Independent term.
    pub coef0: F,
}

impl<F: Float> SigmoidKernel<F> {
    /// Create a new sigmoid kernel with defaults (`gamma='scale'`, `coef0=0`).
    #[must_use]
    pub fn new() -> Self {
        Self {
            gamma: Gamma::Scale,
            coef0: F::zero(),
        }
    }
}

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

impl<F: Float + Send + Sync> Kernel<F> for SigmoidKernel<F> {
    fn compute(&self, x: &[F], y: &[F]) -> F {
        let gamma = gamma_value_or_one(self.gamma);
        let dot: F = x
            .iter()
            .zip(y.iter())
            .fold(F::zero(), |acc, (&a, &b)| acc + a * b);
        (gamma * dot + self.coef0).tanh()
    }

    fn resolved_for_fit(&self, x: &Array2<F>) -> Self {
        Self {
            gamma: Gamma::Value(resolve_gamma(self.gamma, x)),
            coef0: self.coef0,
        }
    }
}

// ---------------------------------------------------------------------------
// Kernel cache (LRU)
// ---------------------------------------------------------------------------

/// Simple LRU cache for kernel evaluations.
struct KernelCache<F> {
    cache: HashMap<(usize, usize), F>,
    order: Vec<(usize, usize)>,
    capacity: usize,
}

impl<F: Float> KernelCache<F> {
    fn new(capacity: usize) -> Self {
        Self {
            cache: HashMap::with_capacity(capacity),
            order: Vec::with_capacity(capacity),
            capacity,
        }
    }

    fn get_or_compute<K: Kernel<F>>(
        &mut self,
        i: usize,
        j: usize,
        kernel: &K,
        data: &[Vec<F>],
    ) -> F {
        let key = if i <= j { (i, j) } else { (j, i) };
        if let Some(&val) = self.cache.get(&key) {
            return val;
        }
        let val = kernel.compute(&data[i], &data[j]);
        if self.order.len() >= self.capacity
            && let Some(old_key) = self.order.first().copied()
        {
            self.cache.remove(&old_key);
            self.order.remove(0);
        }
        self.cache.insert(key, val);
        self.order.push(key);
        val
    }
}

// ---------------------------------------------------------------------------
// SMO solver for binary SVM
// ---------------------------------------------------------------------------

/// Result of a binary SMO solve.
struct SmoResult<F> {
    alphas: Vec<F>,
    bias: F,
}

/// SMO implementation (Platt 1998, Fan-Chen-Lin 2005 WSS).
///
/// Uses the dual gradient `grad_i = (Q * alpha)_i - 1` where
/// `Q_{ij} = y_i * y_j * K(x_i, x_j)`. Bias is computed after
/// convergence from the KKT conditions.
#[allow(
    clippy::too_many_arguments,
    reason = "the per-class box bounds (cp, cn) are separate args mirroring \
              libsvm's per-sample upper bound C_i (Cp for y=+1, Cn for y=-1)"
)]
fn smo_binary<F: Float, K: Kernel<F>>(
    data: &[Vec<F>],
    labels: &[F],
    kernel: &K,
    cp: F,
    cn: F,
    tol: F,
    max_iter: usize,
    cache_size: usize,
) -> Result<SmoResult<F>, FerroError> {
    let n = data.len();
    let mut alphas = vec![F::zero(); n];
    let mut cache = KernelCache::new(cache_size);

    // Per-sample box upper bound `C_i = (y_i > 0 ? Cp : Cn)` (libsvm `GETI`):
    // `class_weight` scales C per class so the +1 group (class_pos) gets `Cp`
    // and the -1 group (class_neg) gets `Cn`. When `cp == cn` the box is the
    // uniform `[0, C]` of the no-class-weight case.
    let c_of = |i: usize| -> F { if labels[i] > F::zero() { cp } else { cn } };

    // Gradient of the dual objective: grad_i = (Q*alpha)_i - 1
    // where Q_{ij} = y_i * y_j * K(x_i, x_j).
    // Initially alpha = 0, so grad_i = -1 for all i.
    let mut grad: Vec<F> = vec![-F::one(); n];

    let two = F::one() + F::one();
    let eps = F::from(1e-12).unwrap_or_else(F::epsilon);

    // `max_iter == 0` is the sklearn `max_iter=-1` ("no iteration limit",
    // libsvm runs to convergence) sentinel — the SMO loop then runs until the
    // KKT gap closes. A non-zero `max_iter` caps the iteration count.
    let mut iter = 0usize;
    loop {
        if max_iter != 0 && iter >= max_iter {
            break;
        }
        iter += 1;
        // Working set selection (Fan-Chen-Lin 2005):
        // I_up  = {i : (y_i=+1 and alpha_i < C) or (y_i=-1 and alpha_i > 0)}
        // I_low = {j : (y_j=+1 and alpha_j > 0) or (y_j=-1 and alpha_j < C)}
        // Select i = argmax_{t in I_up}  -y_t * grad_t
        // Select j = argmin_{t in I_low} -y_t * grad_t

        let mut i_up = None;
        let mut max_val = F::neg_infinity();
        let mut j_low = None;
        let mut min_val = F::infinity();

        for t in 0..n {
            let val = -labels[t] * grad[t];
            let c_t = c_of(t);

            let in_up = (labels[t] > F::zero() && alphas[t] < c_t - eps)
                || (labels[t] < F::zero() && alphas[t] > eps);

            let in_low = (labels[t] > F::zero() && alphas[t] > eps)
                || (labels[t] < F::zero() && alphas[t] < c_t - eps);

            if in_up && val > max_val {
                max_val = val;
                i_up = Some(t);
            }
            if in_low && val < min_val {
                min_val = val;
                j_low = Some(t);
            }
        }

        // Stopping criterion: KKT gap < tol
        if i_up.is_none() || j_low.is_none() || max_val - min_val < tol {
            break;
        }

        let i = i_up.unwrap();
        let j = j_low.unwrap();

        if i == j {
            break;
        }

        // Compute second-order info
        let kii = cache.get_or_compute(i, i, kernel, data);
        let kjj = cache.get_or_compute(j, j, kernel, data);
        let kij = cache.get_or_compute(i, j, kernel, data);
        let eta = kii + kjj - two * kij;

        if eta <= eps {
            continue;
        }

        // Bounds for alpha_j, respecting the per-sample box bounds
        // `0 <= alpha_i <= C_i` and `0 <= alpha_j <= C_j` (libsvm allows a
        // different upper bound per sample under `class_weight`).
        let old_ai = alphas[i];
        let old_aj = alphas[j];
        let ci = c_of(i);
        let cj = c_of(j);

        let (lo, hi) = if labels[i] == labels[j] {
            // alpha_i + alpha_j = sum (const): alpha_j in
            // [max(0, sum - C_i), min(C_j, sum)].
            let sum = old_ai + old_aj;
            ((sum - ci).max(F::zero()), sum.min(cj))
        } else {
            // alpha_i = alpha_j - diff (const diff): alpha_j in
            // [max(0, diff), min(C_j, C_i + diff)].
            let diff = old_aj - old_ai;
            (diff.max(F::zero()), (ci + diff).min(cj))
        };

        if (hi - lo).abs() < eps {
            continue;
        }

        // Analytic update for alpha_j (Platt 1998).
        // E_k = y_k * grad_k (dual error, where grad = Q*alpha - e).
        // alpha_j_new = alpha_j + y_j * (E_i - E_j) / eta
        //             = alpha_j + y_j * (y_i * grad_i - y_j * grad_j) / eta
        let mut new_aj = old_aj + labels[j] * (labels[i] * grad[i] - labels[j] * grad[j]) / eta;

        // Clip to [lo, hi]
        if new_aj > hi {
            new_aj = hi;
        }
        if new_aj < lo {
            new_aj = lo;
        }

        if (new_aj - old_aj).abs() < eps {
            continue;
        }

        let new_ai = old_ai + labels[i] * labels[j] * (old_aj - new_aj);

        alphas[i] = new_ai;
        alphas[j] = new_aj;

        // Update dual gradient: grad_k += delta_alpha_i * Q_{k,i} + delta_alpha_j * Q_{k,j}
        // where Q_{k,t} = y_k * y_t * K(k,t)
        let delta_ai = new_ai - old_ai;
        let delta_aj = new_aj - old_aj;

        for (k, grad_k) in grad.iter_mut().enumerate() {
            let kki = cache.get_or_compute(k, i, kernel, data);
            let kkj = cache.get_or_compute(k, j, kernel, data);
            *grad_k = *grad_k
                + delta_ai * labels[k] * labels[i] * kki
                + delta_aj * labels[k] * labels[j] * kkj;
        }
    }

    // Compute bias from KKT conditions.
    // For support vectors with 0 < alpha_i < C:
    //   y_i * (sum_j alpha_j * y_j * K(i,j) + b) = 1
    //   b = y_i - sum_j alpha_j * y_j * K(i,j)
    // (since y_i^2 = 1, y_i * (y_i * f) = f, so b = 1/y_i - sum = y_i - sum)
    let mut b_sum = F::zero();
    let mut b_count = 0usize;

    for i in 0..n {
        if alphas[i] > eps && alphas[i] < c_of(i) - eps {
            // This is a free support vector (`0 < alpha_i < C_i`).
            let mut f_no_b = F::zero();
            for j in 0..n {
                if alphas[j] > eps {
                    f_no_b =
                        f_no_b + alphas[j] * labels[j] * cache.get_or_compute(i, j, kernel, data);
                }
            }
            b_sum = b_sum + labels[i] - f_no_b;
            b_count += 1;
        }
    }

    let bias = if b_count > 0 {
        b_sum / F::from(b_count).unwrap()
    } else {
        // Fallback: use all support vectors (bounded ones too)
        let mut b_sum_all = F::zero();
        let mut b_count_all = 0usize;
        for i in 0..n {
            if alphas[i] > eps {
                let mut f_no_b = F::zero();
                for j in 0..n {
                    if alphas[j] > eps {
                        f_no_b = f_no_b
                            + alphas[j] * labels[j] * cache.get_or_compute(i, j, kernel, data);
                    }
                }
                b_sum_all = b_sum_all + labels[i] - f_no_b;
                b_count_all += 1;
            }
        }
        if b_count_all > 0 {
            b_sum_all / F::from(b_count_all).unwrap()
        } else {
            F::zero()
        }
    };

    Ok(SmoResult { alphas, bias })
}

// ---------------------------------------------------------------------------
// Platt scaling (probability estimates)
// ---------------------------------------------------------------------------

/// Fit the Platt sigmoid `P(y=+1 | f) = 1 / (1 + exp(A·f + B))` to a set of
/// decision values `dec_values` with binary labels `labels` (`+1` / `-1`),
/// returning the `(A, B)` parameters.
///
/// A faithful transcription of libsvm's `sigmoid_train`
/// (`sklearn/svm/src/libsvm/svm.cpp:1919-2030`): the prior-based initial point
/// (`A=0`, `B=log((prior0+1)/(prior1+1))`), the `t` target smoothing
/// (`hiTarget=(prior1+1)/(prior1+2)`, `loTarget=1/(prior0+2)`), the Newton
/// iteration with the regularized Hessian (`H' = H + sigma·I`,
/// `sigma=1e-12`), the gradient/Hessian accumulation, the step-halving line
/// search (`min_step=1e-10`, sufficient-decrease constant `0.0001`),
/// `max_iter=100`, and the `eps=1e-5` gradient stopping criterion. The
/// overflow-safe `fApB>=0` branching matches the C code exactly.
#[allow(
    clippy::too_many_lines,
    reason = "a faithful one-to-one transcription of libsvm's sigmoid_train \
              Newton loop (svm.cpp:1919-2030); splitting it would obscure the \
              line-by-line correspondence to the C oracle"
)]
fn sigmoid_train<F: Float>(dec_values: &[F], labels: &[F]) -> (F, F) {
    let l = dec_values.len();
    let zero = F::zero();
    let one = F::one();
    let two = one + one;

    let mut prior1 = zero;
    let mut prior0 = zero;
    for &lab in labels {
        if lab > zero {
            prior1 = prior1 + one;
        } else {
            prior0 = prior0 + one;
        }
    }

    let max_iter = 100usize;
    let min_step = F::from(1e-10).unwrap_or_else(F::epsilon);
    let sigma = F::from(1e-12).unwrap_or_else(F::epsilon);
    let eps = F::from(1e-5).unwrap_or_else(F::epsilon);
    let suff = F::from(0.0001).unwrap_or_else(F::epsilon);

    let hi_target = (prior1 + one) / (prior1 + two);
    let lo_target = one / (prior0 + two);

    // Per-sample target smoothed labels `t`.
    let t: Vec<F> = labels
        .iter()
        .map(|&lab| if lab > zero { hi_target } else { lo_target })
        .collect();

    // Initial point and initial function value.
    let mut a = zero;
    let mut b = ((prior0 + one) / (prior1 + one)).ln();

    let funcval = |a: F, b: F| -> F {
        let mut fval = zero;
        for i in 0..l {
            let f_ap_b = dec_values[i] * a + b;
            if f_ap_b >= zero {
                fval = fval + t[i] * f_ap_b + (one + (-f_ap_b).exp()).ln();
            } else {
                fval = fval + (t[i] - one) * f_ap_b + (one + f_ap_b.exp()).ln();
            }
        }
        fval
    };

    let mut fval = funcval(a, b);

    for _iter in 0..max_iter {
        // Update gradient and Hessian (H' = H + sigma·I).
        let mut h11 = sigma;
        let mut h22 = sigma;
        let mut h21 = zero;
        let mut g1 = zero;
        let mut g2 = zero;
        for i in 0..l {
            let f_ap_b = dec_values[i] * a + b;
            let (p, q) = if f_ap_b >= zero {
                let e = (-f_ap_b).exp();
                (e / (one + e), one / (one + e))
            } else {
                let e = f_ap_b.exp();
                (one / (one + e), e / (one + e))
            };
            let d2 = p * q;
            h11 = h11 + dec_values[i] * dec_values[i] * d2;
            h22 = h22 + d2;
            h21 = h21 + dec_values[i] * d2;
            let d1 = t[i] - p;
            g1 = g1 + dec_values[i] * d1;
            g2 = g2 + d1;
        }

        // Stopping criterion.
        if g1.abs() < eps && g2.abs() < eps {
            break;
        }

        // Newton direction: -inv(H')·g.
        let det = h11 * h22 - h21 * h21;
        let d_a = -(h22 * g1 - h21 * g2) / det;
        let d_b = -(-h21 * g1 + h11 * g2) / det;
        let gd = g1 * d_a + g2 * d_b;

        // Line search (step halving).
        let mut stepsize = one;
        while stepsize >= min_step {
            let new_a = a + stepsize * d_a;
            let new_b = b + stepsize * d_b;
            let newf = funcval(new_a, new_b);
            if newf < fval + suff * stepsize * gd {
                a = new_a;
                b = new_b;
                fval = newf;
                break;
            }
            stepsize = stepsize / two;
        }

        if stepsize < min_step {
            // Line search failed — libsvm bails out of the Newton loop.
            break;
        }
    }

    (a, b)
}

/// Evaluate the Platt sigmoid `P(y=+1 | f) = 1 / (1 + exp(A·f + B))` at a single
/// decision value, in the overflow-safe form of libsvm's `sigmoid_predict`
/// (`sklearn/svm/src/libsvm/svm.cpp:2032-2040`):
/// `fApB = decision·A + B`; if `fApB >= 0` return `exp(-fApB)/(1+exp(-fApB))`,
/// else `1/(1+exp(fApB))` (avoiding `exp` overflow / catastrophic
/// cancellation).
fn sigmoid_predict<F: Float>(decision: F, a: F, b: F) -> F {
    let f_ap_b = decision * a + b;
    if f_ap_b >= F::zero() {
        let e = (-f_ap_b).exp();
        e / (F::one() + e)
    } else {
        F::one() / (F::one() + f_ap_b.exp())
    }
}

/// Wu-Lin-Weng (2004) pairwise coupling ("Method 2"): given the `k×k` pairwise
/// probability matrix `r` (where `r[i][j] = P(class i | class i or j)`),
/// produce the `k` coupled class probabilities `p`.
///
/// A faithful transcription of libsvm's `multiclass_probability`
/// (`sklearn/svm/src/libsvm/svm.cpp:2043-2104`): build the `Q` matrix from the
/// pairwise probabilities, then run the fixed-point iteration
/// (`max_iter = max(100, k)`, `eps = 0.005/k`) that minimizes the coupling
/// objective, normalized so the returned probabilities sum to 1.
fn multiclass_probability<F: Float>(k: usize, r: &Array2<F>) -> Vec<F> {
    let zero = F::zero();
    let one = F::one();
    let k_f = F::from(k).unwrap_or(one);

    let mut p = vec![one / k_f; k];
    // Q[t][j].
    let mut q = Array2::<F>::zeros((k, k));
    for t in 0..k {
        for j in 0..t {
            q[[t, t]] = q[[t, t]] + r[[j, t]] * r[[j, t]];
            q[[t, j]] = q[[j, t]];
        }
        for j in (t + 1)..k {
            q[[t, t]] = q[[t, t]] + r[[j, t]] * r[[j, t]];
            q[[t, j]] = -r[[j, t]] * r[[t, j]];
        }
    }

    let max_iter = 100.max(k);
    let eps = F::from(0.005).unwrap_or_else(F::epsilon) / k_f;
    let mut qp = vec![zero; k];

    for _iter in 0..max_iter {
        // Recompute Qp, pQp for numerical accuracy.
        let mut p_qp = zero;
        for t in 0..k {
            qp[t] = zero;
            for j in 0..k {
                qp[t] = qp[t] + q[[t, j]] * p[j];
            }
            p_qp = p_qp + p[t] * qp[t];
        }
        let mut max_error = zero;
        for &qpt in qp.iter().take(k) {
            let error = (qpt - p_qp).abs();
            if error > max_error {
                max_error = error;
            }
        }
        if max_error < eps {
            break;
        }

        for t in 0..k {
            let qtt = q[[t, t]];
            if qtt == zero {
                continue;
            }
            let diff = (-qp[t] + p_qp) / qtt;
            p[t] = p[t] + diff;
            p_qp = (p_qp + diff * (diff * qtt + two_qp(qp[t]))) / (one + diff) / (one + diff);
            for j in 0..k {
                qp[j] = (qp[j] + diff * q[[t, j]]) / (one + diff);
                p[j] = p[j] / (one + diff);
            }
        }
    }

    p
}

/// `2·x` helper for [`multiclass_probability`] (libsvm `2*Qp[t]`).
#[inline]
fn two_qp<F: Float>(x: F) -> F {
    x + x
}

/// Decision value of a freshly-trained binary SMO sub-model on a query sample,
/// in ferrolearn's sign convention (positive favors the `+1` label, i.e. the
/// higher-index `class_pos` group).
fn sub_decision_value<F: Float, K: Kernel<F>>(
    sv_data: &[Vec<F>],
    sv_coefs: &[F],
    bias: F,
    kernel: &K,
    q: &[F],
) -> F {
    let mut val = bias;
    for (sv, &coef) in sv_data.iter().zip(sv_coefs.iter()) {
        val = val + coef * kernel.compute(sv, q);
    }
    val
}

/// A freshly-trained binary sub-model in this crate's (ferrolearn) sign
/// convention: support-vector feature rows, their coefficients
/// (`alpha_i·y_i`, `class_pos = +1` side), and the decision bias such that
/// [`sub_decision_value`] is positive favoring `class_pos`. Returned by the
/// per-solver TRAINER closure that [`platt_cv_sigmoid`] invokes on each CV
/// training fold.
pub(crate) type SubModel<F> = (Vec<Vec<F>>, Vec<F>, F);

/// Fit the per-ovo-pair Platt sigmoid `(A, B)` via a DETERMINISTIC 5-fold CV
/// over the pair's samples, mirroring libsvm's `svm_binary_svc_probability`
/// (`sklearn/svm/src/libsvm/svm.cpp:2107-2203`) EXCEPT for the fold
/// permutation.
///
/// libsvm shuffles the fold assignment with an RNG seeded by `random_state`
/// (`svm.cpp:2116-2122`), which makes the resulting `(A, B)` (sklearn's
/// `probA_`/`probB_`) and thus `predict_proba` NON-DETERMINISTIC across
/// `random_state`. To keep ferrolearn deterministic (it has no libsvm RNG
/// seed; cf. the documented SGD shuffle boundary, R-DEV-4), the folds here use a
/// DETERMINISTIC CLASS-STRATIFIED assignment instead of libsvm's random shuffle:
/// each sample's fold is its WITHIN-CLASS running index modulo `nr_fold`. Because
/// the per-ovo-pair samples arrive GROUPED by class, a naive contiguous
/// `[i·l/5, (i+1)·l/5)` split would make whole folds single-class — so the 4-fold
/// training set could miss a class entirely and `sigmoid_train` would collapse to
/// the degenerate `(A, B) = (0, 0)`. Stratifying within class keeps both classes
/// in every training set (when each class has ≥2 samples), restoring libsvm's
/// structural contract (a non-degenerate sigmoid) without its randomness. The
/// rest is a faithful transcription: train a binary sub-model on the 4 training
/// folds,
/// `predict_values` the held-out fold (in libsvm sign), with the degenerate
/// one-class-fold fallbacks (`+1` / `-1` / `0`, `svm.cpp:2161-2169`), then
/// [`sigmoid_train`] over all out-of-fold decisions.
///
/// # The `train_fold` trainer abstraction
///
/// libsvm's `svm_binary_svc_probability` trains each CV sub-model with the
/// SAME `svm_type` as the outer model (`svm.cpp:2147-2150`, a copy of the
/// outer `svm_parameter` with `probability=0`): C-SVC sub-models for `SVC`,
/// NU-SVC sub-models for `NuSVC`. ferrolearn threads that choice through a
/// `train_fold` closure: given the training-fold `(data, labels)` (in
/// ferrolearn sign, `class_pos = +1`), it returns the fitted [`SubModel`]
/// (`Some`) or `None` on a degenerate/failed sub-solve. `SVC` passes a
/// closure wrapping [`smo_binary`] (C-SVC); `NuSVC` passes a closure wrapping
/// [`solve_nu_svc`] (the genuine `Solver_NU`). The CV split, degenerate-fold
/// fallbacks, held-out scoring via [`sub_decision_value`], and the final
/// [`sigmoid_train`] are SOLVER-AGNOSTIC, so SVC's `(A, B)` is byte-identical
/// to the pre-refactor inline-`smo_binary` path.
///
/// `sub_labels` is ferrolearn's sign (`+1` = higher-index `class_pos`,
/// `-1` = lower-index `class_neg`). The decision values and labels passed to
/// [`sigmoid_train`] are converted to libsvm sign (`+1` = lower-index
/// `class_neg`, matching `raw_ovo`) so the fitted `(A, B)` is consistent with
/// the raw ovo decision used by [`FittedSVC::predict_proba`].
pub(crate) fn platt_cv_sigmoid<F: Float, K: Kernel<F>>(
    sub_data: &[Vec<F>],
    sub_labels: &[F],
    kernel: &K,
    train_fold: impl Fn(&[Vec<F>], &[F]) -> Option<SubModel<F>>,
) -> (F, F) {
    let l = sub_data.len();
    let nr_fold = 5usize;
    // Out-of-fold decision value per sample, in libsvm sign (+1 = class_neg).
    let mut dec_values = vec![F::zero(); l];

    // DETERMINISTIC class-stratified fold assignment. libsvm shuffles the fold
    // permutation with an RNG (`svm.cpp:2116-2122`) so each fold mixes both
    // classes; ferrolearn stays deterministic (no libsvm RNG seed; cf. the
    // sanctioned SGD-shuffle boundary, R-DEV-4) by instead assigning each sample
    // to a fold by its WITHIN-CLASS running index modulo `nr_fold`. The
    // per-ovo-pair samples arrive GROUPED by class (`[class_neg..., class_pos...]`,
    // built by the `FittedSVC::fit` loop), so a CONTIGUOUS `[i·l/5, (i+1)·l/5)`
    // split would make whole folds single-class and the 4-fold training set could
    // MISS a class entirely → a trivial sub-model → constant held-out decisions →
    // `sigmoid_train` returns the degenerate `(A, B) = (0, 0)`. Spreading each
    // class proportionally across all folds keeps BOTH classes in every training
    // set whenever each class has ≥2 samples, restoring the structural contract
    // (a non-degenerate sigmoid) at every input — matching libsvm's intent
    // without its randomness.
    let mut pos_seen = 0usize;
    let mut neg_seen = 0usize;
    let mut fold_of = vec![0usize; l];
    for (j, &lab) in sub_labels.iter().enumerate() {
        if lab > F::zero() {
            fold_of[j] = pos_seen % nr_fold;
            pos_seen += 1;
        } else {
            fold_of[j] = neg_seen % nr_fold;
            neg_seen += 1;
        }
    }

    for fold in 0..nr_fold {
        // Training set = all samples NOT assigned to this fold.
        let mut tr_data: Vec<Vec<F>> = Vec::with_capacity(l);
        let mut tr_labels: Vec<F> = Vec::with_capacity(l);
        for (j, row) in sub_data.iter().enumerate() {
            if fold_of[j] != fold {
                tr_data.push(row.clone());
                tr_labels.push(sub_labels[j]);
            }
        }

        // Count classes in the training folds (ferrolearn sign).
        let mut p_count = 0usize;
        let mut n_count = 0usize;
        for &lab in &tr_labels {
            if lab > F::zero() {
                p_count += 1;
            } else {
                n_count += 1;
            }
        }

        // Degenerate folds: libsvm assigns a constant decision
        // (`svm.cpp:2161-2169`). In ferrolearn sign a held-out sample gets
        // +1 (all-positive train), -1 (all-negative train), or 0 (empty); we
        // store the libsvm-sign value = negation. The held-out fold is now the
        // (non-contiguous) set `{ j : fold_of[j] == fold }`, not a slice.
        let held_out = (0..l).filter(|&j| fold_of[j] == fold);
        if p_count == 0 && n_count == 0 {
            for j in held_out {
                dec_values[j] = F::zero();
            }
            continue;
        } else if n_count == 0 {
            // train all +1 (class_pos) -> ferrolearn dec +1 -> libsvm -1.
            for j in held_out {
                dec_values[j] = -F::one();
            }
            continue;
        } else if p_count == 0 {
            for j in held_out {
                dec_values[j] = F::one();
            }
            continue;
        }

        // Train a probability-free sub-model on the training folds via the
        // per-solver trainer (C-SVC for SVC, NU-SVC for NuSVC).
        let Some((sv_data, sv_coefs, bias)) = train_fold(&tr_data, &tr_labels) else {
            // A failed/degenerate sub-solve falls back to a neutral 0 decision.
            for j in held_out {
                dec_values[j] = F::zero();
            }
            continue;
        };

        // Score the held-out fold; store in libsvm sign (negate ferrolearn).
        for j in held_out {
            let dec_ferro = sub_decision_value(&sv_data, &sv_coefs, bias, kernel, &sub_data[j]);
            dec_values[j] = -dec_ferro;
        }
    }

    // libsvm labels for sigmoid_train: +1 = lower-index class_neg, matching
    // the libsvm-sign decision values (so `-sub_labels`).
    let libsvm_labels: Vec<F> = sub_labels.iter().map(|&lab| -lab).collect();
    sigmoid_train(&dec_values, &libsvm_labels)
}

// ---------------------------------------------------------------------------
// decision_function shape + scores
// ---------------------------------------------------------------------------

/// The shape convention for [`FittedSVC::decision_function`] in the multiclass
/// case, mirroring scikit-learn's `SVC.decision_function_shape`
/// (`sklearn/svm/_base.py:778-781`).
///
/// - [`SvmDecisionShape::Ovr`] (default): one-vs-rest scores, shape
///   `(n_samples, n_classes)`, produced by the `_ovr_decision_function`
///   transform (`sklearn/utils/multiclass.py:520-562`).
/// - [`SvmDecisionShape::Ovo`]: the raw one-vs-one decision values, shape
///   `(n_samples, n_class·(n_class-1)/2)`.
///
/// The binary case is unaffected (it always collapses to a 1-D `(n_samples,)`
/// score, `_base.py:538-539`).
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum SvmDecisionShape {
    /// One-vs-rest: `(n_samples, n_classes)` via `_ovr_decision_function`
    /// (sklearn's default).
    #[default]
    Ovr,
    /// One-vs-one: raw `(n_samples, n_class·(n_class-1)/2)` decision values.
    Ovo,
}

/// The result of [`FittedSVC::decision_function`].
///
/// Mirrors scikit-learn's polymorphic `SVC.decision_function` return
/// (`sklearn/svm/_base.py:536-541, 778-781`): the binary case collapses the
/// single ovo column to a 1-D `(n_samples,)` array (`-dec.ravel()`,
/// `_base.py:538-539`), while the multiclass case returns
/// `(n_samples, n_classes)` (ovr, default) or
/// `(n_samples, n·(n-1)/2)` (ovo). Structurally parallels
/// [`crate::linear_svc::DecisionScores`] for cross-estimator consistency.
#[derive(Debug, Clone, PartialEq)]
pub enum SvmScores<F> {
    /// Binary decision values, shape `(n_samples,)`. A POSITIVE value predicts
    /// `classes_[1]` (`-dec.ravel()`, `_base.py:538-539`).
    Binary(Array1<F>),
    /// Multiclass decision values: `(n_samples, n_classes)` for
    /// [`SvmDecisionShape::Ovr`] or `(n_samples, n·(n-1)/2)` for
    /// [`SvmDecisionShape::Ovo`].
    Multiclass(Array2<F>),
}

impl<F: Clone> SvmScores<F> {
    /// Number of samples scored (the leading axis length in both variants).
    #[must_use]
    pub fn n_samples(&self) -> usize {
        match self {
            SvmScores::Binary(v) => v.len(),
            SvmScores::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 {
            SvmScores::Binary(v) => Some(v),
            SvmScores::Multiclass(_) => None,
        }
    }

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

// ---------------------------------------------------------------------------
// SVC (Support Vector Classifier)
// ---------------------------------------------------------------------------

/// Support Vector Classifier.
///
/// Uses Sequential Minimal Optimization (SMO) to solve the dual QP.
/// Supports multiclass via one-vs-one strategy.
///
/// # Type Parameters
///
/// - `F`: The floating-point type (`f32` or `f64`).
/// - `K`: The kernel type (e.g., [`LinearKernel`], [`RbfKernel`]).
#[derive(Debug, Clone)]
pub struct SVC<F, K> {
    /// The kernel function.
    pub kernel: K,
    /// Regularization parameter (penalty for misclassification).
    pub c: F,
    /// Convergence tolerance.
    pub tol: F,
    /// Maximum number of SMO iterations. `0` is the sklearn `max_iter=-1`
    /// sentinel meaning **no iteration limit** (the SMO runs to convergence);
    /// a non-zero value caps the iteration count
    /// (`sklearn/svm/_classes.py`, `max_iter` default `-1`).
    pub max_iter: usize,
    /// Size of the kernel evaluation LRU cache (perf-only; default `200` to
    /// match sklearn's `cache_size=200`).
    pub cache_size: usize,
    /// Whether to use libsvm's shrinking heuristic
    /// (`sklearn/svm/_base.py:339`, `_classes.py` `shrinking=True`).
    ///
    /// ferrolearn's SMO has no shrinking heuristic: shrinking is a libsvm
    /// performance optimization that does NOT change the converged optimum
    /// (R-DEV-7). This flag is accepted for API parity (default `true`,
    /// matching sklearn) but DOES NOT alter the fitted result — the converged
    /// `α`/`dual_coef_`/`intercept_` are shrinking-invariant.
    pub shrinking: bool,
    /// The multiclass `decision_function` shape convention
    /// (`sklearn/svm/_base.py:778-781`); default
    /// [`SvmDecisionShape::Ovr`] (sklearn's `decision_function_shape='ovr'`).
    pub decision_function_shape: SvmDecisionShape,
    /// Whether `predict` breaks ties by the one-vs-rest decision confidence
    /// instead of the libsvm vote (`break_ties`, `sklearn/svm/_classes.py`
    /// default `False`; semantics in `BaseSVC.predict`,
    /// `sklearn/svm/_base.py:801-814`).
    ///
    /// When `true` AND [`SvmDecisionShape::Ovr`] AND `n_classes > 2`,
    /// `predict = argmax(decision_function(X))` (the ovr decision, which breaks
    /// ties by confidence); otherwise the libsvm ovo vote (with lower-index
    /// tie-break) is used. `break_ties=true` with [`SvmDecisionShape::Ovo`] is
    /// rejected at predict time (`InvalidParameter`), matching sklearn
    /// (`_base.py:801-804`).
    pub break_ties: bool,
    /// Per-class scaling of `C` (`class_weight`, `sklearn/svm/_classes.py:118-124`).
    /// Default [`ClassWeight::None`] (all classes weighted `1.0`). For an ovo
    /// pair `(a, b)` with `a < b`, the C of the `y=+1` group (class `b`) is
    /// `C·class_weight_[b]` and the C of the `y=-1` group (class `a`) is
    /// `C·class_weight_[a]`; the weights are computed ONCE over the full `y`
    /// by [`compute_class_weight`] (`_base.py:740`).
    pub class_weight: ClassWeight<F>,
    /// Whether to enable Platt-scaling probability estimates
    /// (`probability`, `sklearn/svm/_classes.py`, default `False`).
    ///
    /// When `true`, [`Fit::fit`] runs an internal 5-fold cross-validation per
    /// one-vs-one pair, fits a sigmoid `1/(1+exp(A·f+B))` over the out-of-fold
    /// decision values ([`sigmoid_train`], libsvm `svm_binary_svc_probability`,
    /// `svm.cpp:2107-2203`), and stores the per-pair `(A, B)` so
    /// [`FittedSVC::predict_proba`]/[`FittedSVC::predict_log_proba`] are
    /// available. When `false` (the default) `predict_proba` returns an error
    /// (`_base.py:820-827`).
    ///
    /// **RNG boundary (documented divergence).** libsvm's
    /// `svm_binary_svc_probability` shuffles the CV fold assignment with an
    /// RNG seeded by `random_state`, so sklearn's `probA_`/`probB_` (and hence
    /// the exact `predict_proba` values) are NON-DETERMINISTIC across
    /// `random_state` — the docstring itself warns "the results can be slightly
    /// different than those obtained by predict". ferrolearn instead uses a
    /// DETERMINISTIC 5-fold split (contiguous folds, no RNG shuffle), so it
    /// CANNOT and DOES NOT bit-match sklearn's `predict_proba` values. What is
    /// reproduced exactly is the DETERMINISTIC machinery ([`sigmoid_train`],
    /// [`sigmoid_predict`], [`multiclass_probability`]) and the STRUCTURAL
    /// contract (rows sum to 1, entries in `[0, 1]`, monotone in the binary
    /// decision value, the raise-when-`probability=false`). This is analogous
    /// to the SGD shuffle boundary already documented in this codebase.
    pub probability: bool,
}

impl<F: Float, K: Kernel<F>> SVC<F, K> {
    /// Create a new `SVC` with the given kernel and default hyperparameters
    /// matching sklearn (`sklearn/svm/_classes.py` `SVC.__init__`).
    ///
    /// Defaults: `C = 1.0`, `tol = 1e-3`, `max_iter = 0` (= sklearn `-1`, no
    /// iteration limit), `cache_size = 200`, `shrinking = true`,
    /// `decision_function_shape = Ovr`, `break_ties = false`,
    /// `class_weight = None`, `probability = false`.
    #[must_use]
    pub fn new(kernel: K) -> Self {
        Self {
            kernel,
            c: F::one(),
            tol: F::from(1e-3).unwrap_or_else(F::epsilon),
            max_iter: 0,
            cache_size: 200,
            shrinking: true,
            decision_function_shape: SvmDecisionShape::Ovr,
            break_ties: false,
            class_weight: ClassWeight::None,
            probability: false,
        }
    }

    /// Enable/disable Platt-scaling probability estimates (`sklearn`
    /// `probability`, default `false`). When `true`, [`Fit::fit`] runs the
    /// internal per-pair 5-fold CV + [`sigmoid_train`] so
    /// [`FittedSVC::predict_proba`]/[`FittedSVC::predict_log_proba`] are
    /// available; when `false` they return an error.
    ///
    /// See the [`SVC::probability`] field doc for the documented RNG-CV
    /// exact-value boundary (sklearn is non-deterministic across
    /// `random_state`; only the deterministic machinery + structural
    /// invariants + the raise contract are reproduced).
    #[must_use]
    pub fn with_probability(mut self, probability: bool) -> Self {
        self.probability = probability;
        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 `shrinking` flag (`sklearn` `shrinking`, default `true`).
    ///
    /// Accepted for API parity; does NOT alter the converged optimum
    /// (ferrolearn's SMO has no shrinking heuristic — R-DEV-7).
    #[must_use]
    pub fn with_shrinking(mut self, shrinking: bool) -> Self {
        self.shrinking = shrinking;
        self
    }

    /// Set the `break_ties` flag (`sklearn` `break_ties`, default `false`,
    /// `sklearn/svm/_base.py:801-814`).
    #[must_use]
    pub fn with_break_ties(mut self, break_ties: bool) -> Self {
        self.break_ties = break_ties;
        self
    }

    /// Set the multiclass `decision_function` shape convention
    /// (`'ovr'` default / `'ovo'`, `sklearn/svm/_base.py:778-781`).
    #[must_use]
    pub fn with_decision_function_shape(mut self, shape: SvmDecisionShape) -> Self {
        self.decision_function_shape = shape;
        self
    }

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

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

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

    /// Set the kernel cache size.
    #[must_use]
    pub fn with_cache_size(mut self, cache_size: usize) -> Self {
        self.cache_size = cache_size;
        self
    }
}

/// A single binary SVM model (one pair of classes in one-vs-one).
#[derive(Debug, Clone)]
struct BinarySvm<F> {
    /// Support vectors (stored as rows).
    support_vectors: Vec<Vec<F>>,
    /// Original training-row index of each support vector (parallel to
    /// `support_vectors`/`dual_coefs`). Used to build the global, per-class
    /// grouped `support_` set (`sklearn/svm/_base.py:318-410`).
    sv_indices: Vec<usize>,
    /// Dual coefficients: `alpha_i * y_i` for each support vector, where this
    /// crate maps the lower-index class (`class_neg`) to `y = -1` and the
    /// higher-index class (`class_pos`) to `y = +1`. NOTE this is the OPPOSITE
    /// sign convention to libsvm internally (libsvm gives the lower-index class
    /// `+1`); the public-attribute layout compensates in
    /// [`FittedSVC::dual_coef`].
    dual_coefs: Vec<F>,
    /// Bias term.
    bias: F,
    /// The two class labels: (negative_class, positive_class). `class_neg` is
    /// the lower class index and `class_pos` the higher (the ovo pair `(a, b)`
    /// with `a < b`).
    class_neg: usize,
    class_pos: usize,
}

/// Fitted Support Vector Classifier.
///
/// Stores one binary SVM per pair of classes (one-vs-one). Implements
/// [`Predict`] to produce class labels.
#[derive(Debug, Clone)]
pub struct FittedSVC<F, K> {
    /// The kernel used for predictions.
    kernel: K,
    /// One binary SVM per class pair, in libsvm ovo pair order
    /// `(0,1),(0,2),...,(0,k-1),(1,2),...` (the `(ci,cj)` double loop).
    binary_models: Vec<BinarySvm<F>>,
    /// Sorted unique classes (`classes_ = np.unique(y)`).
    classes: Vec<usize>,
    /// The training feature matrix, retained so the libsvm-layout fitted
    /// attributes (`support_`, `support_vectors_`) can index back into the
    /// original rows (`sklearn/svm/_base.py:318-410`).
    x_train: Array2<F>,
    /// The training labels (class index per row), retained so `support_` can
    /// be grouped by class.
    y_train: Vec<usize>,
    /// The multiclass `decision_function` shape convention carried over from
    /// the unfitted [`SVC`] (`sklearn/svm/_base.py:778-781`).
    decision_function_shape: SvmDecisionShape,
    /// The `break_ties` flag carried over from the unfitted [`SVC`]
    /// (`sklearn/svm/_base.py:801-814`).
    break_ties: bool,
    /// Whether probability estimates were fitted (`probability`,
    /// `sklearn/svm/_classes.py`). When `false`, [`Self::predict_proba`]
    /// returns an error (`_base.py:820-827`).
    probability: bool,
    /// Per-ovo-pair Platt sigmoid `A` parameter (`probA_`,
    /// `sklearn/svm/src/libsvm/svm.cpp:2200` via `sigmoid_train`), parallel to
    /// `binary_models`. Empty when `probability == false`.
    prob_a: Vec<F>,
    /// Per-ovo-pair Platt sigmoid `B` parameter (`probB_`), parallel to
    /// `binary_models`. Empty when `probability == false`.
    prob_b: Vec<F>,
}

/// One ovo binary sub-model in **this crate's sign convention** (higher-index
/// `class_pos` is the `+1` side, matching [`BinarySvm`] and
/// [`FittedSVC::decision_value_binary`]). Used by [`FittedSVC::from_nu_ovo`] to
/// assemble a nu-SVC fitted model that reuses all of [`FittedSVC`]'s accessors
/// / `decision_function` / `predict`.
///
/// The nu-SVC solver ([`solve_nu_svc`]) is fed the per-pair labels in this same
/// convention (`class_pos = +1`), so `sv_coefs`/`bias_internal` are already in
/// this-crate sign and `from_nu_ovo` stores them verbatim.
pub(crate) struct NuOvoPair<F> {
    /// Support-vector feature rows for this pair.
    pub sv_data: Vec<Vec<F>>,
    /// Per-SV coefficient `alpha·y/r` (this-crate sign, `class_pos = +1`),
    /// equal to the public binary `dual_coef_` value (the nu_svc binary flip
    /// `public = -internal` cancels with `internal = -stored`).
    pub sv_coefs: Vec<F>,
    /// Original training-row index of each support vector.
    pub sv_indices: Vec<usize>,
    /// Decision bias for the `+1`-side (`class_pos`) in this crate's
    /// convention (`f(x) = Σ sv_coef·K + bias_internal`).
    pub bias_internal: F,
    /// Lower-index class label (this crate's `-1` side).
    pub class_neg: usize,
    /// Higher-index class label (this crate's `+1` side).
    pub class_pos: usize,
}

impl<F: Float + Send + Sync + ScalarOperand + 'static, K: Kernel<F> + 'static> FittedSVC<F, K> {
    /// Assemble a [`FittedSVC`] from per-ovo-pair nu-SVC sub-models (in libsvm
    /// sign convention) so that [`NuSVC`](crate::nu_svm::NuSVC) reuses the full
    /// libsvm-layout fitted-attribute machinery (`support_`/`dual_coef_`/
    /// `intercept_`/`coef_`/`decision_function`/`predict`) without duplicating
    /// it (`sklearn/svm/_base.py:318-410`).
    ///
    /// Each [`NuOvoPair`] is already in this crate's [`BinarySvm`] sign
    /// convention (higher-index `class_pos` as the `+1` side, because
    /// [`solve_nu_svc`] is fed labels in that convention), so the coefficients
    /// and bias are stored verbatim. The resulting public `dual_coef_`/
    /// `intercept_` then carry the binary nu_svc sign flip exactly as `c_svc`
    /// does (`_base.py:258-262`, predicate `_impl in ["c_svc","nu_svc"]`).
    ///
    /// When `probability` is `true`, `prob_a`/`prob_b` are the per-ovo-pair
    /// sigmoid `(A, B)` parameters fitted by [`platt_cv_sigmoid`] with the
    /// NU-SVC sub-solver ([`solve_nu_svc`]) ([`NuSVC`](crate::nu_svm::NuSVC)
    /// REQ-9), so the assembled [`FittedSVC`]'s
    /// [`Self::predict_proba`]/[`Self::predict_log_proba`] work identically to
    /// a `probability=true` C-SVC fit — the coupling/`sigmoid_predict` path is
    /// solver-agnostic (it consumes only `raw_ovo` + `prob_a`/`prob_b`). When
    /// `probability` is `false`, `prob_a`/`prob_b` are empty.
    #[allow(
        clippy::too_many_arguments,
        reason = "carries the full nu-ovo assembly plus the probability state \
                  (probability flag + per-pair probA/probB) in one constructor"
    )]
    pub(crate) fn from_nu_ovo(
        kernel: K,
        pairs: Vec<NuOvoPair<F>>,
        classes: Vec<usize>,
        x_train: Array2<F>,
        y_train: Vec<usize>,
        decision_function_shape: SvmDecisionShape,
        break_ties: bool,
        probability: bool,
        prob_a: Vec<F>,
        prob_b: Vec<F>,
    ) -> Self {
        let binary_models = pairs
            .into_iter()
            .map(|pair| BinarySvm {
                support_vectors: pair.sv_data,
                sv_indices: pair.sv_indices,
                dual_coefs: pair.sv_coefs,
                bias: pair.bias_internal,
                class_neg: pair.class_neg,
                class_pos: pair.class_pos,
            })
            .collect();

        FittedSVC {
            kernel,
            binary_models,
            classes,
            x_train,
            y_train,
            decision_function_shape,
            break_ties,
            probability,
            prob_a,
            prob_b,
        }
    }
}

impl<F: Float, K: Kernel<F>> FittedSVC<F, K> {
    /// Compute the decision function value for a single sample against a
    /// binary model.
    fn decision_value_binary(&self, x: &[F], model: &BinarySvm<F>) -> F {
        let mut val = model.bias;
        for (sv, &coef) in model.support_vectors.iter().zip(model.dual_coefs.iter()) {
            val = val + coef * self.kernel.compute(sv, x);
        }
        val
    }

    /// Raw one-vs-one decision values in **libsvm sign convention**, shape
    /// `(n_samples, n·(n-1)/2)`, columns in ovo pair order
    /// `(0,1),(0,2),...,(1,2),...` (the `(ci,cj)` double loop).
    ///
    /// libsvm/sklearn use the LOWER-index class as the `+1` side, so a POSITIVE
    /// value means the lower-index class wins (`sklearn/svm/_base.py:520-524`).
    /// This crate's [`Self::decision_value_binary`] uses the HIGHER-index class
    /// (`class_pos`) as `+1`, the opposite sign — so the raw ovo value is the
    /// negation of `decision_value_binary` to restore libsvm's convention.
    fn raw_ovo(&self, x: &Array2<F>) -> Array2<F> {
        let n_samples = x.nrows();
        let n_models = self.binary_models.len();
        let mut result = Array2::<F>::zeros((n_samples, n_models));

        for s in 0..n_samples {
            let xi: Vec<F> = x.row(s).to_vec();
            for (m, model) in self.binary_models.iter().enumerate() {
                // Negate to match libsvm's "lower-index class = +1" sign.
                result[[s, m]] = self.decision_value_binary(&xi, model).neg();
            }
        }

        result
    }

    /// The continuous one-vs-rest decision function transformed from the
    /// one-vs-one scores, mirroring `_ovr_decision_function`
    /// (`sklearn/utils/multiclass.py:520-562`), shape `(n_samples, n_classes)`.
    ///
    /// `predictions[s,k] = if raw_ovo[s,k] < 0 { 1 } else { 0 }` and
    /// `confidences[s,k] = -raw_ovo[s,k]`, matching sklearn's call
    /// `_ovr_decision_function(dec < 0, -dec, n_classes)`
    /// (`sklearn/svm/_base.py:780`). The ovo pair iteration `(i,j)` with `i<j`,
    /// `k = 0,1,...`, is the SAME order as the `raw_ovo` columns.
    fn ovr_decision_function(&self, raw_ovo: &Array2<F>) -> Array2<F> {
        let n_samples = raw_ovo.nrows();
        let n_classes = self.classes.len();
        let mut votes = Array2::<F>::zeros((n_samples, n_classes));
        let mut sum_of_confidences = Array2::<F>::zeros((n_samples, n_classes));
        let one = F::one();

        let mut k = 0usize;
        for i in 0..n_classes {
            for j in (i + 1)..n_classes {
                for s in 0..n_samples {
                    let dec = raw_ovo[[s, k]];
                    let confidence = dec.neg(); // -dec
                    // sum_of_confidences[:, i] -= confidences[:, k]
                    sum_of_confidences[[s, i]] = sum_of_confidences[[s, i]] - confidence;
                    // sum_of_confidences[:, j] += confidences[:, k]
                    sum_of_confidences[[s, j]] = sum_of_confidences[[s, j]] + confidence;
                    // predictions[s,k] = (dec < 0) ? 1 : 0
                    // votes[predictions==0, i] += 1; votes[predictions==1, j] += 1
                    if dec < F::zero() {
                        votes[[s, j]] = votes[[s, j]] + one;
                    } else {
                        votes[[s, i]] = votes[[s, i]] + one;
                    }
                }
                k += 1;
            }
        }

        // transformed = sum_of_confidences / (3 * (|sum_of_confidences| + 1))
        // return votes + transformed.
        let three = match F::from(3.0) {
            Some(v) => v,
            None => one + one + one,
        };
        let mut out = votes;
        for s in 0..n_samples {
            for c in 0..n_classes {
                let soc = sum_of_confidences[[s, c]];
                let transformed = soc / (three * (soc.abs() + one));
                out[[s, c]] = out[[s, c]] + transformed;
            }
        }
        out
    }

    /// The decision function for the samples in `x`
    /// (`sklearn/svm/_base.py:536-541, 778-781`).
    ///
    /// - **Binary** (`n_classes == 2`): [`SvmScores::Binary`], shape
    ///   `(n_samples,)` = `-raw_ovo.ravel()` (sklearn flips the sign for
    ///   `c_svc`/`nu_svc` binary, `_base.py:538-539`), so a POSITIVE value
    ///   predicts `classes_[1]`. Because this crate's `decision_value_binary`
    ///   already uses the higher-index class as `+1`, `-raw_ovo` equals
    ///   `decision_value_binary` directly.
    /// - **Multiclass [`SvmDecisionShape::Ovr`]** (default):
    ///   [`SvmScores::Multiclass`], shape `(n_samples, n_classes)` =
    ///   `_ovr_decision_function(raw_ovo)` (`_base.py:780`).
    /// - **Multiclass [`SvmDecisionShape::Ovo`]**: [`SvmScores::Multiclass`],
    ///   shape `(n_samples, n·(n-1)/2)` = the raw ovo values.
    ///
    /// # Errors
    ///
    /// Returns `Ok` for valid input.
    pub fn decision_function(&self, x: &Array2<F>) -> Result<SvmScores<F>, FerroError> {
        let raw_ovo = self.raw_ovo(x);

        if self.classes.len() == 2 {
            // Binary: -raw_ovo.ravel() = +decision_value_binary (1-D).
            let n_samples = raw_ovo.nrows();
            let mut binary = Array1::<F>::zeros(n_samples);
            for s in 0..n_samples {
                binary[s] = raw_ovo[[s, 0]].neg();
            }
            return Ok(SvmScores::Binary(binary));
        }

        match self.decision_function_shape {
            SvmDecisionShape::Ovo => Ok(SvmScores::Multiclass(raw_ovo)),
            SvmDecisionShape::Ovr => {
                Ok(SvmScores::Multiclass(self.ovr_decision_function(&raw_ovo)))
            }
        }
    }

    /// Class probability estimates, shape `(n_samples, n_classes)`; columns
    /// correspond to `classes_` in sorted order
    /// (`sklearn/svm/_base.py:829-864`, `libsvm.predict_probability`,
    /// `svm.cpp:2918-2955`).
    ///
    /// Built from the per-pair Platt sigmoids fitted at `fit` time when
    /// `probability=true`: the raw one-vs-one decision values (libsvm sign,
    /// lower-index class `+1`) are mapped to pairwise probabilities via
    /// [`sigmoid_predict`] (clamped to `[1e-7, 1-1e-7]`, `svm.cpp:2929-2938`),
    /// then coupled by [`multiclass_probability`] (Wu-Lin-Weng 2004,
    /// `svm.cpp:2941`). For the binary case `multiclass_probability` reduces to
    /// `[sigmoid_predict(dec), 1 - sigmoid_predict(dec)]` =
    /// `[P(classes_[0]), P(classes_[1])]`. Each row sums to 1.
    ///
    /// **RNG-CV exact-value boundary (documented divergence).** Because the
    /// underlying `(A, B)` come from a cross-validation whose fold assignment
    /// is RNG-dependent in libsvm/sklearn (sklearn's `probA_`/`probB_` and
    /// `predict_proba` values change with `random_state`), ferrolearn uses a
    /// DETERMINISTIC contiguous 5-fold split instead and therefore does NOT
    /// bit-match sklearn's `predict_proba` values. The reproduced contract is
    /// structural: rows sum to 1, entries in `[0, 1]`, and (binary) the
    /// `classes_[1]` column is monotone non-decreasing in the
    /// `decision_function` value. See [`SVC::probability`].
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::InvalidParameter`] when the model was fitted with
    /// `probability=false`, with the message mirroring sklearn's
    /// `NotFittedError` text "predict_proba is not available when fitted with
    /// probability=False" (`_base.py:856-860`). (This crate has no `NotFitted`
    /// variant — predict-before-fit is a compile error via the typestate,
    /// R-DEV-4 — so the "fitted-without-probability" runtime condition maps to
    /// `InvalidParameter`; the binding maps it to the matching `PyErr`.)
    pub fn predict_proba(&self, x: &Array2<F>) -> Result<Array2<F>, FerroError> {
        if !self.probability {
            return Err(FerroError::InvalidParameter {
                name: "probability".into(),
                reason: "predict_proba is not available when fitted with probability=False".into(),
            });
        }

        let raw_ovo = self.raw_ovo(x);
        let n_samples = raw_ovo.nrows();
        let n_classes = self.classes.len();
        let min_prob = F::from(1e-7).unwrap_or_else(F::epsilon);
        let max_prob = F::one() - min_prob;

        let mut out = Array2::<F>::zeros((n_samples, n_classes));

        for s in 0..n_samples {
            // Build the k×k pairwise probability matrix for this sample.
            let mut pairwise = Array2::<F>::zeros((n_classes, n_classes));
            let mut k = 0usize;
            for i in 0..n_classes {
                for j in (i + 1)..n_classes {
                    // dec_values[k] is the raw ovo value (libsvm sign: positive
                    // favors the lower-index class i = classes_[i]).
                    let dec = raw_ovo[[s, k]];
                    let (a, b) = (self.prob_a[k], self.prob_b[k]);
                    let mut pij = sigmoid_predict(dec, a, b);
                    // Clamp to [min_prob, 1-min_prob] (`svm.cpp:2937`).
                    if pij < min_prob {
                        pij = min_prob;
                    }
                    if pij > max_prob {
                        pij = max_prob;
                    }
                    pairwise[[i, j]] = pij;
                    pairwise[[j, i]] = F::one() - pij;
                    k += 1;
                }
            }
            let probs = multiclass_probability(n_classes, &pairwise);
            for (c, &pc) in probs.iter().enumerate() {
                out[[s, c]] = pc;
            }
        }

        Ok(out)
    }

    /// Natural-log class probability estimates, shape `(n_samples, n_classes)`
    /// = `predict_proba(x).ln()` elementwise (`sklearn/svm/_base.py:866-894`:
    /// `np.log(self.predict_proba(X))`).
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::NotFitted`] when the model was fitted with
    /// `probability=false` (delegated from [`Self::predict_proba`]).
    pub fn predict_log_proba(&self, x: &Array2<F>) -> Result<Array2<F>, FerroError> {
        self.predict_proba(x).map(|p| p.mapv(F::ln))
    }

    /// Whether Platt-scaling probability estimates were fitted
    /// (`probability`, `sklearn/svm/_classes.py`); when `false`,
    /// [`Self::predict_proba`]/[`Self::predict_log_proba`] raise.
    #[must_use]
    pub fn probability(&self) -> bool {
        self.probability
    }

    /// The per-ovo-pair Platt sigmoid `A` parameters (`probA_`,
    /// `sklearn/svm/_base.py`), length `n_class·(n_class-1)/2`. Empty when
    /// fitted with `probability=false`.
    #[must_use]
    pub fn prob_a(&self) -> Array1<F> {
        Array1::from_vec(self.prob_a.clone())
    }

    /// The per-ovo-pair Platt sigmoid `B` parameters (`probB_`,
    /// `sklearn/svm/_base.py`), length `n_class·(n_class-1)/2`. Empty when
    /// fitted with `probability=false`.
    #[must_use]
    pub fn prob_b(&self) -> Array1<F> {
        Array1::from_vec(self.prob_b.clone())
    }
}

impl<F: Float + ScalarOperand + 'static, K: Kernel<F>> FittedSVC<F, K> {
    /// Build the global, per-class-grouped support-vector index set, mirroring
    /// libsvm's `support_` layout (`sklearn/svm/_base.py:318-410`): the indices
    /// of the training rows that are a support vector in AT LEAST ONE ovo
    /// binary problem, deduplicated, grouped by class (all of class
    /// `classes_[0]` first, then `classes_[1]`, ...), ascending within a class.
    ///
    /// Returns `(support, per_class_indices)` where `support` is the flat
    /// grouped index vector and `per_class_indices[c]` is the (ascending)
    /// list of global training-row indices that are SVs for class
    /// `classes_[c]`.
    fn build_support(&self) -> (Vec<usize>, Vec<Vec<usize>>) {
        let n_classes = self.classes.len();
        // Per-class set of training-row indices that are an SV anywhere.
        let mut per_class: Vec<Vec<usize>> = vec![Vec::new(); n_classes];
        let mut seen: Vec<bool> = vec![false; self.y_train.len()];

        for model in &self.binary_models {
            for &idx in &model.sv_indices {
                if !seen[idx] {
                    seen[idx] = true;
                    let cls = self.y_train[idx];
                    if let Some(ci) = self.classes.iter().position(|&c| c == cls) {
                        per_class[ci].push(idx);
                    }
                }
            }
        }

        for group in &mut per_class {
            group.sort_unstable();
        }

        let support: Vec<usize> = per_class.iter().flatten().copied().collect();
        (support, per_class)
    }

    /// Indices of the support vectors into the training set, **grouped by
    /// class** (all class-`classes_[0]` SVs first, then `classes_[1]`, ...),
    /// ascending within each class.
    ///
    /// Mirrors `SVC.support_` (`sklearn/svm/_base.py:318-410`).
    #[must_use]
    pub fn support(&self) -> Array1<usize> {
        let (support, _) = self.build_support();
        Array1::from_vec(support)
    }

    /// The support vectors `X[support_]`, shape `(n_SV, n_features)`.
    ///
    /// Mirrors `SVC.support_vectors_` (`sklearn/svm/_base.py:318-410`).
    #[must_use]
    pub fn support_vectors(&self) -> Array2<F> {
        let (support, _) = self.build_support();
        let n_features = self.x_train.ncols();
        let mut out = Array2::<F>::zeros((support.len(), n_features));
        for (row, &idx) in support.iter().enumerate() {
            out.row_mut(row).assign(&self.x_train.row(idx));
        }
        out
    }

    /// Number of support vectors per class (`n_support_`,
    /// `sklearn/svm/_base.py:668-682`), parallel to `classes_`.
    #[must_use]
    pub fn n_support(&self) -> Vec<usize> {
        let (_, per_class) = self.build_support();
        per_class.iter().map(Vec::len).collect()
    }

    /// Dual coefficients in the libsvm public layout, shape
    /// `(n_class - 1, n_SV)` (`sklearn/svm/_base.py:318-410`, the `dual_coef_`
    /// attribute), columns in `support_` (per-class-grouped) order.
    ///
    /// For an SV belonging to class `i`, row `m` holds its coefficient in the
    /// binary classifier between class `i` and the `m`-th OTHER class (the
    /// other classes in increasing index order, skipping `i`). In the ovo pair
    /// `(a, b)` with `a < b`, libsvm uses class `a` as the `+1` side and `b` as
    /// `-1`; the stored coefficient is `alpha * y_libsvm`.
    ///
    /// This crate stores `alpha * y` per pair with the OPPOSITE sign
    /// (`class_neg = a` mapped to `-1`), so the libsvm-internal coefficient is
    /// the negation of the stored value. For `n_class == 2` sklearn negates the
    /// internal coefficient again to form the PUBLIC binary attribute
    /// (`sklearn/svm/_base.py:258-262`: `dual_coef_ = -dual_coef_`), which
    /// leaves the public binary value equal to this crate's stored value; for
    /// `n_class > 2` the public value IS the libsvm-internal value (no flip).
    #[must_use]
    pub fn dual_coef(&self) -> Array2<F> {
        let n_classes = self.classes.len();
        let (support, _per_class) = self.build_support();
        let n_sv = support.len();

        if n_classes == 2 {
            // Binary: public dual_coef_ = -internal, and internal = -stored,
            // so public = stored. The single ovo model holds one stored coef
            // per SV keyed by training index; map them into support_ column
            // order.
            let mut out = Array2::<F>::zeros((1, n_sv));
            if let Some(model) = self.binary_models.first() {
                for (sv_idx, &coef) in model.sv_indices.iter().zip(model.dual_coefs.iter()) {
                    if let Some(col) = support.iter().position(|&s| s == *sv_idx) {
                        out[[0, col]] = coef;
                    }
                }
            }
            return out;
        }

        // Multiclass: public dual_coef_ = libsvm-internal = -(stored). Row m
        // for an SV of class i is its coefficient in the pair (i, m-th other).
        let mut out = Array2::<F>::zeros((n_classes - 1, n_sv));

        // Column index in `support_` for a given training-row index.
        let col_of: HashMap<usize, usize> =
            support.iter().enumerate().map(|(c, &i)| (i, c)).collect();

        for model in &self.binary_models {
            let a = model.class_neg; // lower-index class in the pair
            let b = model.class_pos; // higher-index class
            let ai = self.classes.iter().position(|&c| c == a);
            let bi = self.classes.iter().position(|&c| c == b);
            let (ai, bi) = match (ai, bi) {
                (Some(ai), Some(bi)) => (ai, bi),
                _ => continue,
            };

            for (sv_idx, &stored) in model.sv_indices.iter().zip(model.dual_coefs.iter()) {
                let Some(&col) = col_of.get(sv_idx) else {
                    continue;
                };
                let cls = self.y_train[*sv_idx];
                let internal = stored.neg(); // libsvm internal = -(stored)
                // Determine which row this pair occupies for class `cls`:
                // the count of OTHER classes with index < the partner's index.
                let (own_class_index, partner_class_index) =
                    if cls == a { (ai, bi) } else { (bi, ai) };
                // Row m = number of other classes (excluding own) with class
                // index < partner_class_index.
                let mut row = 0usize;
                for ci in 0..n_classes {
                    if ci == own_class_index {
                        continue;
                    }
                    if ci < partner_class_index {
                        row += 1;
                    }
                }
                out[[row, col]] = internal;
            }
        }

        out
    }

    /// The per-ovo-problem intercepts, length `n_class * (n_class - 1) / 2`,
    /// in pair order `(0,1),(0,2),(1,2),...` (`intercept_`,
    /// `sklearn/svm/_base.py:318-410`). For `n_class == 2` sklearn negates the
    /// internal bias to form the public attribute
    /// (`sklearn/svm/_base.py:258-262`: `intercept_ *= -1`).
    ///
    /// This crate's per-pair `bias` is recovered for the decision function
    /// `sum coef*K + bias` with `class_pos` (the higher index) as the `+1`
    /// side. libsvm/sklearn use the lower-index class as `+1`, so the
    /// libsvm-internal intercept is the negation of this crate's `bias`. For
    /// binary, the public attribute negates the internal again, leaving the
    /// public value equal to this crate's stored `bias`; for multiclass the
    /// public value IS the internal `-bias`.
    #[must_use]
    pub fn intercept(&self) -> Array1<F> {
        let n_classes = self.classes.len();
        let vals: Vec<F> = if n_classes == 2 {
            self.binary_models.iter().map(|m| m.bias).collect()
        } else {
            self.binary_models.iter().map(|m| m.bias.neg()).collect()
        };
        Array1::from_vec(vals)
    }

    /// Primal weight vector `coef_ = dual_coef_ @ support_vectors_`, shape
    /// `(n_class - 1, n_features)` — available ONLY for the linear kernel
    /// (`sklearn/svm/_base.py:650-666`). Returns `None` for any other kernel
    /// (sklearn raises `AttributeError`).
    #[must_use]
    pub fn coef(&self) -> Option<Array2<F>> {
        if !self.kernel.is_linear() {
            return None;
        }
        let dual = self.dual_coef(); // (n_class-1, n_SV)
        let svs = self.support_vectors(); // (n_SV, n_features)
        Some(dual.dot(&svs))
    }
}

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

    /// Fit the SVC model using SMO.
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if `x` and `y` have different
    /// sample counts.
    /// Returns [`FerroError::InvalidParameter`] if `C` is not positive.
    /// Returns [`FerroError::InsufficientSamples`] if fewer than 2 classes.
    fn fit(&self, x: &Array2<F>, y: &Array1<usize>) -> Result<FittedSVC<F, K>, FerroError> {
        let (n_samples, _n_features) = x.dim();

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

        // Reject non-finite input (NaN / +/-inf) in X BEFORE the X/y length
        // (`ShapeMismatch`) check, mirroring sklearn's `BaseLibSVM.fit` ->
        // `_validate_data(X, y, ...)` (`sklearn/svm/_base.py:190`) which routes to
        // `check_X_y`: `check_array(X, force_all_finite=True)` runs BEFORE
        // `check_consistent_length(X, y)` (`_base.py:208`), so on an input that is
        // BOTH non-finite AND length-mismatched sklearn raises the finiteness
        // `ValueError("Input X contains NaN.")` / `"... contains infinity ..."`,
        // NOT a consistency error (#2270). `y` is class labels
        // (`Array1<usize>`), already finite by type, so only `X` needs the float
        // finiteness check. `.iter().any(|v| !v.is_finite())` catches both NaN and
        // +/-inf; on finite input the guard never fires, so the fitted SVC
        // attributes (`support_`/`dual_coef_`/`intercept_`/`coef_`) are
        // byte-identical and a finite length-mismatch still yields `ShapeMismatch`.
        if x.iter().any(|v| !v.is_finite()) {
            return Err(FerroError::InvalidParameter {
                name: "X".into(),
                reason: "Input X contains NaN or infinity.".into(),
            });
        }

        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(),
            });
        }

        // Determine unique classes.
        let mut classes: Vec<usize> = y.to_vec();
        classes.sort_unstable();
        classes.dedup();

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

        // Resolve any data-dependent kernel parameters (e.g. a `None` gamma ->
        // sklearn's default `gamma='scale'` = 1/(n_features * X.var()),
        // `_base.py:236-239`) against the training data BEFORE fitting, and use
        // this resolved kernel for both fitting and prediction.
        let kernel = self.kernel.resolved_for_fit(x);

        // Convert data to Vec<Vec<F>> for kernel cache.
        let data: Vec<Vec<F>> = (0..n_samples).map(|i| x.row(i).to_vec()).collect();

        // Per-class weights computed ONCE over the FULL y (libsvm's
        // `class_weight_ = compute_class_weight(class_weight, classes, y)`,
        // `_base.py:740`); `weighted_C[c] = C·class_weight_[c]`.
        let y_vec: Vec<usize> = y.to_vec();
        let weights = compute_class_weight(&self.class_weight, &classes, &y_vec);

        // One-vs-one: train one binary SVM per pair.
        let n_classes = classes.len();
        let mut binary_models = Vec::new();
        // Per-pair Platt sigmoid params (only filled when `probability`).
        let mut prob_a: Vec<F> = Vec::new();
        let mut prob_b: Vec<F> = Vec::new();

        for ci in 0..n_classes {
            for cj in (ci + 1)..n_classes {
                let class_neg = classes[ci];
                let class_pos = classes[cj];

                // Per-class box bounds for this ovo pair `(class_neg, class_pos)`:
                // the `y=+1` group (class_pos = classes[cj]) gets `Cp = C·w[cj]`
                // and the `y=-1` group (class_neg = classes[ci]) gets
                // `Cn = C·w[ci]` (`weighted_C`, `_base.py:740`). The `weights`
                // vector is aligned to `classes`, so the class-index = the
                // position in `classes` (`ci`/`cj`).
                let cp = self.c * weights[cj];
                let cn = self.c * weights[ci];

                // Extract samples for these two classes.
                let mut sub_data = Vec::new();
                let mut sub_labels = Vec::new();
                let mut sub_indices = Vec::new();

                for s in 0..n_samples {
                    let label = y[s];
                    if label == class_neg {
                        sub_data.push(data[s].clone());
                        sub_labels.push(-F::one());
                        sub_indices.push(s);
                    } else if label == class_pos {
                        sub_data.push(data[s].clone());
                        sub_labels.push(F::one());
                        sub_indices.push(s);
                    }
                }

                let result = smo_binary(
                    &sub_data,
                    &sub_labels,
                    &kernel,
                    cp,
                    cn,
                    self.tol,
                    self.max_iter,
                    self.cache_size,
                )?;

                // Extract support vectors (non-zero alphas).
                let eps = F::from(1e-8).unwrap_or_else(F::epsilon);
                let mut sv_data = Vec::new();
                let mut sv_coefs = Vec::new();
                let mut sv_idx = Vec::new();

                for (k, &alpha) in result.alphas.iter().enumerate() {
                    if alpha > eps {
                        sv_data.push(sub_data[k].clone());
                        sv_coefs.push(alpha * sub_labels[k]);
                        // Record the ORIGINAL training-row index of this SV
                        // (sub_indices maps the per-pair row k back to X).
                        sv_idx.push(sub_indices[k]);
                    }
                }

                // Platt-scaling CV for this ovo pair (only when probability).
                // The CV sub-models are C-SVC (the SAME svm_type as the outer
                // model, libsvm `svm.cpp:2147-2150`): the `train_fold` closure
                // wraps `smo_binary` + SV extraction, returning the fitted
                // sub-model in this crate's sign (`class_pos = +1`).
                if self.probability {
                    let (tol, max_iter, cache_size) = (self.tol, self.max_iter, self.cache_size);
                    let sub_eps = F::from(1e-8).unwrap_or_else(F::epsilon);
                    let (a, b) = platt_cv_sigmoid(
                        &sub_data,
                        &sub_labels,
                        &kernel,
                        |tr_data: &[Vec<F>], tr_labels: &[F]| {
                            let sub = smo_binary(
                                tr_data, tr_labels, &kernel, cp, cn, tol, max_iter, cache_size,
                            )
                            .ok()?;
                            let mut sv_d: Vec<Vec<F>> = Vec::new();
                            let mut sv_c: Vec<F> = Vec::new();
                            for (k, &alpha) in sub.alphas.iter().enumerate() {
                                if alpha > sub_eps {
                                    sv_d.push(tr_data[k].clone());
                                    sv_c.push(alpha * tr_labels[k]);
                                }
                            }
                            Some((sv_d, sv_c, sub.bias))
                        },
                    );
                    prob_a.push(a);
                    prob_b.push(b);
                }

                binary_models.push(BinarySvm {
                    support_vectors: sv_data,
                    sv_indices: sv_idx,
                    dual_coefs: sv_coefs,
                    bias: result.bias,
                    class_neg,
                    class_pos,
                });
            }
        }

        Ok(FittedSVC {
            kernel,
            binary_models,
            classes,
            x_train: x.clone(),
            y_train: y.to_vec(),
            decision_function_shape: self.decision_function_shape,
            break_ties: self.break_ties,
            probability: self.probability,
            prob_a,
            prob_b,
        })
    }
}

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

    /// Predict class labels for the given feature matrix.
    ///
    /// Uses one-vs-one voting (each binary model casts a vote, the most-voted
    /// class wins, ties broken toward the lower class index), matching libsvm's
    /// `super().predict` (`sklearn/svm/_base.py:813-814`).
    ///
    /// When `break_ties == true` AND [`SvmDecisionShape::Ovr`] AND
    /// `n_classes > 2`, ties are instead broken by the one-vs-rest decision
    /// confidence: `predict = argmax(decision_function(X))`
    /// (`sklearn/svm/_base.py:806-811`).
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if the number of features
    /// does not match the training data. Returns
    /// [`FerroError::InvalidParameter`] when `break_ties == true` and the
    /// decision-function shape is [`SvmDecisionShape::Ovo`]
    /// (`sklearn/svm/_base.py:801-804`).
    fn predict(&self, x: &Array2<F>) -> Result<Array1<usize>, FerroError> {
        let n_samples = x.nrows();
        let n_classes = self.classes.len();

        // sklearn raises when break_ties=True and decision_function_shape='ovo'
        // (`_base.py:801-804`), regardless of n_classes.
        if self.break_ties && self.decision_function_shape == SvmDecisionShape::Ovo {
            return Err(FerroError::InvalidParameter {
                name: "break_ties".into(),
                reason: "break_ties must be False when decision_function_shape is 'ovo'".into(),
            });
        }

        // break_ties=True, ovr, multiclass: predict = argmax(decision_function)
        // (`_base.py:806-811`). The ovr decision breaks vote ties by confidence.
        if self.break_ties && self.decision_function_shape == SvmDecisionShape::Ovr && n_classes > 2
        {
            let scores = self.decision_function(x)?;
            let mc = match scores {
                SvmScores::Multiclass(m) => m,
                // n_classes > 2 always yields the multiclass variant.
                SvmScores::Binary(_) => {
                    return Err(FerroError::InvalidParameter {
                        name: "break_ties".into(),
                        reason: "ovr decision function unavailable for break-tie predict".into(),
                    });
                }
            };
            let mut predictions = Array1::<usize>::zeros(n_samples);
            for s in 0..n_samples {
                // argmax over the row; ties keep the first (lowest) index via a
                // strictly-greater scan, matching numpy's argmax.
                let mut best_idx = 0usize;
                let mut best_val = mc[[s, 0]];
                for c in 1..n_classes {
                    if mc[[s, c]] > best_val {
                        best_val = mc[[s, c]];
                        best_idx = c;
                    }
                }
                predictions[s] = self.classes[best_idx];
            }
            return Ok(predictions);
        }

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

        for s in 0..n_samples {
            let xi: Vec<F> = x.row(s).to_vec();
            let mut votes = vec![0usize; n_classes];

            for model in &self.binary_models {
                let val = self.decision_value_binary(&xi, model);
                let winner = if val >= F::zero() {
                    model.class_pos
                } else {
                    model.class_neg
                };
                if let Some(idx) = self.classes.iter().position(|&c| c == winner) {
                    votes[idx] += 1;
                }
            }

            // libsvm/sklearn break ovo vote ties toward the LOWER class index
            // (`sklearn/svm/_base.py:813-814`: `super().predict` -> libsvm
            // `svm_predict` keeps the first argmax). `classes` is ascending
            // (`np.unique(y)`), so a strictly-greater scan keeps the
            // first/lowest-index maximum — unlike `max_by_key`, which returns
            // the LAST maximum (the higher index).
            let mut best_class_idx = 0usize;
            let mut best_votes = 0usize;
            for (i, &v) in votes.iter().enumerate() {
                if v > best_votes {
                    best_votes = v;
                    best_class_idx = i;
                }
            }

            predictions[s] = self.classes[best_class_idx];
        }

        Ok(predictions)
    }
}

// ---------------------------------------------------------------------------
// SVR (Support Vector Regressor)
// ---------------------------------------------------------------------------

/// Support Vector Regressor.
///
/// Uses SMO to solve the epsilon-SVR dual problem.
///
/// # Type Parameters
///
/// - `F`: The floating-point type (`f32` or `f64`).
/// - `K`: The kernel type.
#[derive(Debug, Clone)]
pub struct SVR<F, K> {
    /// The kernel function.
    pub kernel: K,
    /// Regularization parameter.
    pub c: F,
    /// Epsilon tube width (insensitive loss zone).
    pub epsilon: F,
    /// Convergence tolerance.
    pub tol: F,
    /// Maximum number of SMO iterations. `0` is the sklearn `max_iter=-1`
    /// sentinel meaning **no iteration limit** (the SMO runs to convergence).
    pub max_iter: usize,
    /// Size of the kernel evaluation LRU cache (perf-only; default `200`).
    pub cache_size: usize,
    /// Whether to use libsvm's shrinking heuristic (`sklearn` `shrinking`,
    /// default `true`). Accepted for API parity; does NOT alter the converged
    /// optimum (ferrolearn's SMO has no shrinking heuristic — R-DEV-7).
    pub shrinking: bool,
}

impl<F: Float, K: Kernel<F>> SVR<F, K> {
    /// Create a new `SVR` with the given kernel and default hyperparameters
    /// matching sklearn (`sklearn/svm/_classes.py` `SVR.__init__`).
    ///
    /// Defaults: `C = 1.0`, `epsilon = 0.1`, `tol = 1e-3`, `max_iter = 0`
    /// (= sklearn `-1`, no iteration limit), `cache_size = 200`,
    /// `shrinking = true`.
    #[must_use]
    pub fn new(kernel: K) -> Self {
        Self {
            kernel,
            c: F::one(),
            epsilon: F::from(0.1).unwrap_or_else(F::epsilon),
            tol: F::from(1e-3).unwrap_or_else(F::epsilon),
            max_iter: 0,
            cache_size: 200,
            shrinking: true,
        }
    }

    /// Set the `shrinking` flag (`sklearn` `shrinking`, default `true`).
    ///
    /// Accepted for API parity; does NOT alter the converged optimum
    /// (ferrolearn's SMO has no shrinking heuristic — R-DEV-7).
    #[must_use]
    pub fn with_shrinking(mut self, shrinking: bool) -> Self {
        self.shrinking = shrinking;
        self
    }

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

    /// Set the epsilon tube width.
    #[must_use]
    pub fn with_epsilon(mut self, epsilon: F) -> Self {
        self.epsilon = epsilon;
        self
    }

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

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

    /// Set the kernel cache size.
    #[must_use]
    pub fn with_cache_size(mut self, cache_size: usize) -> Self {
        self.cache_size = cache_size;
        self
    }
}

/// Fitted Support Vector Regressor.
///
/// Stores the support vectors, dual coefficients, and bias.
#[derive(Debug, Clone)]
pub struct FittedSVR<F, K> {
    /// The kernel used for predictions.
    kernel: K,
    /// Support vectors.
    support_vectors: Vec<Vec<F>>,
    /// Original training-row index of each support vector (parallel to
    /// `support_vectors`/`dual_coefs`), for the `support_` attribute.
    sv_indices: Vec<usize>,
    /// Dual coefficients (alpha_i* - alpha_i) for each support vector.
    dual_coefs: Vec<F>,
    /// Bias term.
    bias: F,
}

impl<F: Float, K: Kernel<F>> FittedSVR<F, K> {
    /// Compute the decision function value for a single sample.
    fn decision_value(&self, x: &[F]) -> F {
        let mut val = self.bias;
        for (sv, &coef) in self.support_vectors.iter().zip(self.dual_coefs.iter()) {
            val = val + coef * self.kernel.compute(sv, x);
        }
        val
    }

    /// Compute the raw decision function values for each sample.
    ///
    /// # Errors
    ///
    /// Returns `Ok` always (provided for API symmetry).
    pub fn decision_function(&self, x: &Array2<F>) -> Result<Array1<F>, FerroError> {
        let n_samples = x.nrows();
        let mut result = Array1::<F>::zeros(n_samples);
        for s in 0..n_samples {
            let xi: Vec<F> = x.row(s).to_vec();
            result[s] = self.decision_value(&xi);
        }
        Ok(result)
    }

    /// Indices of the support vectors into the training set, ascending
    /// (`SVR.support_`, `sklearn/svm/_base.py:318-410`). SVR has a single
    /// "class", so there is no per-class grouping; the SVs are kept in
    /// training-row order.
    #[must_use]
    pub fn support(&self) -> Array1<usize> {
        Array1::from_vec(self.sv_indices.clone())
    }

    /// The support vectors, shape `(n_SV, n_features)`
    /// (`SVR.support_vectors_`).
    #[must_use]
    pub fn support_vectors(&self) -> Array2<F> {
        let n_sv = self.support_vectors.len();
        let n_features = self.support_vectors.first().map_or(0, Vec::len);
        let mut out = Array2::<F>::zeros((n_sv, n_features));
        for (r, sv) in self.support_vectors.iter().enumerate() {
            for (c, &v) in sv.iter().enumerate() {
                out[[r, c]] = v;
            }
        }
        out
    }

    /// Number of support vectors. For SVR `n_support_` has size 1
    /// (`sklearn/svm/_base.py:680-682`).
    #[must_use]
    pub fn n_support(&self) -> Vec<usize> {
        vec![self.support_vectors.len()]
    }

    /// Dual coefficients `alpha*_i - alpha_i`, shape `(1, n_SV)`
    /// (`SVR.dual_coef_`). No sign flip applies to SVR
    /// (`sklearn/svm/_base.py:260` restricts the flip to `c_svc`/`nu_svc`).
    #[must_use]
    pub fn dual_coef(&self) -> Array2<F> {
        let n_sv = self.dual_coefs.len();
        let mut out = Array2::<F>::zeros((1, n_sv));
        for (c, &v) in self.dual_coefs.iter().enumerate() {
            out[[0, c]] = v;
        }
        out
    }

    /// The intercept, length 1 (`SVR.intercept_`). The SVR decision function is
    /// `sum coef*K + bias`, matching libsvm's `f(x) = ... + rho`, so the public
    /// intercept equals this crate's stored `bias` (no sign flip).
    #[must_use]
    pub fn intercept(&self) -> Array1<F> {
        Array1::from_vec(vec![self.bias])
    }
}

/// Solve epsilon-SVR using SMO.
///
/// Reformulates the epsilon-SVR dual as a standard 2n-variable QP and
/// solves it with the Fan-Chen-Lin WSS approach, analogous to `smo_binary`.
///
/// The 2n variables are indexed 0..2n:
/// - Index `k` (k < n)  corresponds to alpha\*\_k  with label +1
/// - Index `k` (k >= n) corresponds to alpha\_{k-n} with label -1
///
/// The Q matrix is: `Q_{ij} = s_i * s_j * K(p_i, p_j)` where `s` is the
/// sign (+1 or -1) and `p` maps to the original sample index.
///
/// The linear term is: `q_k = epsilon - y_{p_k} * s_k`.
#[allow(clippy::too_many_arguments)]
fn smo_svr<F: Float, K: Kernel<F>>(
    data: &[Vec<F>],
    targets: &[F],
    kernel: &K,
    c: F,
    epsilon: F,
    tol: F,
    max_iter: usize,
    cache_size: usize,
) -> Result<(Vec<F>, F), FerroError> {
    let n = data.len();
    let m = 2 * n; // Total number of dual variables.

    // Encoding: variable k in [0, m)
    //   k < n  => alpha*_k   (sign = +1, sample index = k)
    //   k >= n => alpha_{k-n} (sign = -1, sample index = k - n)
    //
    // The dual is: min 1/2 * beta^T Q beta + q^T beta
    //   s.t. 0 <= beta_k <= C, sum_k s_k * beta_k = 0
    // where beta_k = alpha*_k or alpha_{k-n},
    //       Q_{ij} = s_i * s_j * K(p_i, p_j),
    //       q_k    = epsilon - y_{p_k} * s_k.
    //
    // This has the same structure as the SVC dual.

    let sign = |k: usize| -> F { if k < n { F::one() } else { -F::one() } };
    let sample = |k: usize| -> usize { if k < n { k } else { k - n } };

    let mut beta = vec![F::zero(); m];
    let mut cache = KernelCache::new(cache_size);

    // Gradient: grad_k = (Q * beta)_k + q_k.  Initially beta=0 so grad_k = q_k.
    // q_k = epsilon - y_{p_k} * s_k
    let mut grad: Vec<F> = (0..m)
        .map(|k| epsilon - targets[sample(k)] * sign(k))
        .collect();

    let two = F::one() + F::one();
    let eps_num = F::from(1e-12).unwrap_or_else(F::epsilon);

    // `max_iter == 0` is the sklearn `max_iter=-1` ("no iteration limit")
    // sentinel — run until the KKT gap closes; non-zero caps the count.
    let mut iter = 0usize;
    loop {
        if max_iter != 0 && iter >= max_iter {
            break;
        }
        iter += 1;
        // WSS: same as SVC but with the extended variables.
        // I_up  = {k : (s_k=+1 and beta_k < C) or (s_k=-1 and beta_k > 0)}
        // I_low = {k : (s_k=+1 and beta_k > 0) or (s_k=-1 and beta_k < C)}
        // Select i = argmax_{k in I_up}  -s_k * grad_k
        // Select j = argmin_{k in I_low} -s_k * grad_k

        let mut i_up = None;
        let mut max_val = F::neg_infinity();
        let mut j_low = None;
        let mut min_val = F::infinity();

        for k in 0..m {
            let sk = sign(k);
            let val = -sk * grad[k];

            let in_up =
                (sk > F::zero() && beta[k] < c - eps_num) || (sk < F::zero() && beta[k] > eps_num);
            let in_low =
                (sk > F::zero() && beta[k] > eps_num) || (sk < F::zero() && beta[k] < c - eps_num);

            if in_up && val > max_val {
                max_val = val;
                i_up = Some(k);
            }
            if in_low && val < min_val {
                min_val = val;
                j_low = Some(k);
            }
        }

        if i_up.is_none() || j_low.is_none() || max_val - min_val < tol {
            break;
        }

        let i = i_up.unwrap();
        let j = j_low.unwrap();

        if i == j {
            break;
        }

        let si = sign(i);
        let sj = sign(j);
        let pi = sample(i);
        let pj = sample(j);

        // Q_{ii} = si*si*K(pi,pi) = K(pi,pi),  similarly for jj and ij
        let kii = cache.get_or_compute(pi, pi, kernel, data);
        let kjj = cache.get_or_compute(pj, pj, kernel, data);
        let kij = cache.get_or_compute(pi, pj, kernel, data);

        // eta = Q_{ii} + Q_{jj} - 2*Q_{ij} = K(pi,pi) + K(pj,pj) - 2*si*sj*K(pi,pj)
        let eta = kii + kjj - two * si * sj * kij;

        if eta <= eps_num {
            continue;
        }

        // Bounds for beta_j
        let old_bi = beta[i];
        let old_bj = beta[j];

        let (lo, hi) = if si == sj {
            let sum = old_bi + old_bj;
            ((sum - c).max(F::zero()), sum.min(c))
        } else {
            let diff = old_bj - old_bi;
            (diff.max(F::zero()), (c + diff).min(c))
        };

        if (hi - lo).abs() < eps_num {
            continue;
        }

        // Analytic update: beta_j += s_j * (E_i - E_j) / eta
        // where E_k = s_k * grad_k
        let mut new_bj = old_bj + sj * (si * grad[i] - sj * grad[j]) / eta;

        if new_bj > hi {
            new_bj = hi;
        }
        if new_bj < lo {
            new_bj = lo;
        }

        if (new_bj - old_bj).abs() < eps_num {
            continue;
        }

        let new_bi = old_bi + si * sj * (old_bj - new_bj);

        beta[i] = new_bi;
        beta[j] = new_bj;

        // Update gradient: grad_k += delta_bi * Q_{k,i} + delta_bj * Q_{k,j}
        // Q_{k,t} = s_k * s_t * K(p_k, p_t)
        let delta_bi = new_bi - old_bi;
        let delta_bj = new_bj - old_bj;

        for (k, grad_k) in grad.iter_mut().enumerate() {
            let sk = sign(k);
            let pk = sample(k);
            let kki = cache.get_or_compute(pk, pi, kernel, data);
            let kkj = cache.get_or_compute(pk, pj, kernel, data);
            *grad_k = *grad_k + delta_bi * sk * si * kki + delta_bj * sk * sj * kkj;
        }
    }

    // Recover alpha*_i = beta_i (i < n) and alpha_i = beta_{i+n} (i >= n).
    // Coefficient for prediction: coef_i = alpha*_i - alpha_i.
    let coefs: Vec<F> = (0..n).map(|i| beta[i] - beta[i + n]).collect();

    // Compute bias from KKT conditions on free support vectors.
    // For k where 0 < beta_k < C:
    //   grad_k = 0 at optimality => (Q*beta)_k + q_k = 0
    //   sum_t beta_t * s_k * s_t * K(p_k, p_t) + epsilon - y_{p_k} * s_k = 0
    //   s_k * sum_t (beta_t * s_t) * K(p_k, p_t) = y_{p_k} * s_k - epsilon
    //   sum_t coef_t_effective * K(p_k, p_t) = y_{p_k} - epsilon / s_k  (nah, let's use f directly)
    //
    // Prediction: f(x) = sum_i coef_i * K(x_i, x) + b
    // For free alpha*_i (0 < alpha*_i < C): y_i - f(x_i) = epsilon  => b = y_i - epsilon - sum coef_j * K(i,j)
    // For free alpha_i  (0 < alpha_i  < C): f(x_i) - y_i = epsilon  => b = y_i + epsilon - sum coef_j * K(i,j)

    let mut b_sum = F::zero();
    let mut b_count = 0usize;

    for i in 0..n {
        let mut kernel_sum = F::zero();
        let has_free = (beta[i] > eps_num && beta[i] < c - eps_num)
            || (beta[i + n] > eps_num && beta[i + n] < c - eps_num);

        if !has_free {
            continue;
        }

        for (j, &cj) in coefs.iter().enumerate() {
            if cj.abs() > eps_num {
                kernel_sum = kernel_sum + cj * cache.get_or_compute(i, j, kernel, data);
            }
        }

        if beta[i] > eps_num && beta[i] < c - eps_num {
            // alpha*_i is free: y_i - f(x_i) = epsilon => b = y_i - epsilon - kernel_sum
            b_sum = b_sum + targets[i] - epsilon - kernel_sum;
            b_count += 1;
        }
        if beta[i + n] > eps_num && beta[i + n] < c - eps_num {
            // alpha_i is free: f(x_i) - y_i = epsilon => b = y_i + epsilon - kernel_sum
            b_sum = b_sum + targets[i] + epsilon - kernel_sum;
            b_count += 1;
        }
    }

    let bias = if b_count > 0 {
        b_sum / F::from(b_count).unwrap()
    } else {
        F::zero()
    };

    Ok((coefs, bias))
}

// ---------------------------------------------------------------------------
// Solver_NU — the libsvm nu-parameterized solver (nu-SVC / nu-SVR)
// ---------------------------------------------------------------------------

/// Output of the generic [`solver_nu_core`] solve.
struct NuResult<F> {
    /// The dual variables `alpha` (length `l`, the solver-internal variables).
    alpha: Vec<F>,
    /// `rho = (r1 - r2) / 2` (the per-pair bias term, `Solver_NU::calculate_rho`
    /// returns this; `sklearn/svm/src/libsvm/svm.cpp:1417`).
    rho: F,
    /// `r = (r1 + r2) / 2` (`si->r`, the nu-SVC `/r` rescale factor /
    /// the nu-SVR `-epsilon`, `svm.cpp:1416`).
    r: F,
}

/// The generic libsvm `Solver_NU` core: solves
/// `min 0.5 αᵀQα + pᵀα  s.t.  0≤α_k≤C_k,  yᵀα=0,  eᵀα=const`
/// where `Q[k][t] = y_k·y_t·K(sample(k), sample(t))`, with the nu-specific
/// second-order working-set selection (separately over the `y=+1` / `y=-1`
/// groups, four running maxima `Gmaxp/Gmaxp2/Gmaxn/Gmaxn2`) and the
/// two-group `rho`/`r` recovery.
///
/// A faithful transcription of libsvm's `Solver_NU` + the shared
/// `Solver::Solve` update step (`sklearn/svm/src/libsvm/svm.cpp:1166-1418`
/// for the nu-specific `select_working_set`/`calculate_rho`, and `:663-940`
/// for the gradient init + analytic 2-variable update + objective). Like the
/// existing [`smo_binary`]/[`smo_svr`] solvers this is a NATURAL-ORDER,
/// no-shrinking, DETERMINISTIC variant: libsvm's shrinking heuristic
/// (`Solver_NU::do_shrinking`, `svm.cpp:1318`) is a performance optimization
/// that does NOT change the converged optimum (R-DEV-7), so it is omitted; the
/// `active_size` always equals `l`. The result (`alpha`, `rho`, `r`) is in
/// libsvm's convention so the caller can reconstruct `dual_coef_`/`intercept_`
/// that match the live `NuSVC`/`NuSVR` oracle after the public sign flip.
///
/// Arguments:
/// - `data`: the ORIGINAL per-sample feature rows (length `n`), keyed by the
///   [`KernelCache`] via `sample(i)` for the kernel evaluation.
/// - `m`: number of solver variables (`l` in libsvm; `= n` for nu-SVC,
///   `= 2n` for nu-SVR).
/// - `sample`: maps a solver variable index `k` to the original sample index
///   (into `data`) used for the kernel evaluation `K(sample(i), sample(j))`.
/// - `y`: the per-variable sign (`+1` / `-1`).
/// - `p`: the linear term (`p[k]`; `0` for nu-SVC, `∓prob.y` for nu-SVR).
/// - `c`: the per-variable upper bound `C_k`.
/// - `alpha`: the (greedily-initialized) starting dual variables, modified
///   in place into the solution.
#[allow(
    clippy::too_many_arguments,
    reason = "a faithful transcription of libsvm's Solver_NU::Solve threads the \
              problem (n, m, sample, y, p, c, alpha) + hyperparameters through \
              one call; splitting them would obscure the C oracle correspondence"
)]
#[allow(
    clippy::too_many_lines,
    reason = "a one-to-one transcription of libsvm's Solver_NU select_working_set \
              + the shared Solver analytic 2-variable update + calculate_rho \
              (svm.cpp:663-940, 1186-1418); kept inline to preserve the \
              line-by-line correspondence to the C oracle"
)]
fn solver_nu_core<F: Float, K: Kernel<F>>(
    data: &[Vec<F>],
    m: usize,
    sample: &dyn Fn(usize) -> usize,
    y: &[F],
    p: &[F],
    c: &[F],
    mut alpha: Vec<F>,
    kernel: &K,
    tol: F,
    max_iter: usize,
    cache_size: usize,
) -> NuResult<F> {
    let zero = F::zero();
    let two = F::one() + F::one();
    // libsvm's TAU (`svm.cpp:316`): a tiny positive floor for the quadratic
    // coefficient when the kernel is not strictly positive-definite.
    let tau = F::from(1e-12).unwrap_or_else(F::epsilon);
    let eps_bound = F::from(1e-12).unwrap_or_else(F::epsilon);

    let mut cache = KernelCache::new(cache_size);
    // `QD[i] = K(sample(i), sample(i))` (since `y_i^2 = 1`, `Q_ii = K(i,i)`).
    // The cache keys on the ORIGINAL-sample index `sample(i)` into `data`.
    let qd: Vec<F> = (0..m)
        .map(|i| {
            let si = sample(i);
            cache.get_or_compute(si, si, kernel, data)
        })
        .collect();

    // `q(i, j) = y_i·y_j·K(sample(i), sample(j))` (libsvm `SVC_Q::get_Q`,
    // `svm.cpp:1436-1446`). The cache is keyed by original-sample index.
    let q_entry = |i: usize, j: usize, cache: &mut KernelCache<F>| -> F {
        y[i] * y[j] * cache.get_or_compute(sample(i), sample(j), kernel, data)
    };

    // Bound predicates (`Solver::update_alpha_status`, `svm.cpp:588-598`).
    let is_upper = |k: usize, a: &[F]| -> bool { a[k] >= c[k] - eps_bound };
    let is_lower = |k: usize, a: &[F]| -> bool { a[k] <= eps_bound };

    // Gradient `G[k] = (Q·alpha)_k + p[k]` (`svm.cpp:693-715`). Since the
    // greedy init has some `alpha_k > 0`, accumulate the full product.
    let mut grad: Vec<F> = p.to_vec();
    #[allow(
        clippy::needless_range_loop,
        reason = "i indexes alpha while the inner loop forms Q[i][j]·alpha[i]"
    )]
    for i in 0..m {
        if alpha[i] > eps_bound {
            let ai = alpha[i];
            for (j, g) in grad.iter_mut().enumerate() {
                *g = *g + ai * q_entry(i, j, &mut cache);
            }
        }
    }

    let mut iter = 0usize;
    loop {
        if max_iter != 0 && iter >= max_iter {
            break;
        }
        iter += 1;

        // ---- Solver_NU::select_working_set (svm.cpp:1186-1296) ----
        let mut gmaxp = F::neg_infinity();
        let mut gmaxp2 = F::neg_infinity();
        let mut gmaxp_idx: isize = -1;
        let mut gmaxn = F::neg_infinity();
        let mut gmaxn2 = F::neg_infinity();
        let mut gmaxn_idx: isize = -1;

        for t in 0..m {
            if y[t] > zero {
                if !is_upper(t, &alpha) && -grad[t] >= gmaxp {
                    gmaxp = -grad[t];
                    gmaxp_idx = t as isize;
                }
            } else if !is_lower(t, &alpha) && grad[t] >= gmaxn {
                gmaxn = grad[t];
                gmaxn_idx = t as isize;
            }
        }

        let ip = gmaxp_idx;
        let in_ = gmaxn_idx;

        let mut gmin_idx: isize = -1;
        let mut obj_diff_min = F::infinity();

        for j in 0..m {
            if y[j] > zero {
                if !is_lower(j, &alpha) {
                    let grad_diff = gmaxp + grad[j];
                    if grad[j] >= gmaxp2 {
                        gmaxp2 = grad[j];
                    }
                    if grad_diff > zero && ip != -1 {
                        let ipi = ip as usize;
                        let quad_coef = qd[ipi] + qd[j] - two * q_entry(ipi, j, &mut cache);
                        let obj_diff = if quad_coef > zero {
                            -(grad_diff * grad_diff) / quad_coef
                        } else {
                            -(grad_diff * grad_diff) / tau
                        };
                        if obj_diff <= obj_diff_min {
                            gmin_idx = j as isize;
                            obj_diff_min = obj_diff;
                        }
                    }
                }
            } else if !is_upper(j, &alpha) {
                let grad_diff = gmaxn - grad[j];
                if -grad[j] >= gmaxn2 {
                    gmaxn2 = -grad[j];
                }
                if grad_diff > zero && in_ != -1 {
                    let ini = in_ as usize;
                    let quad_coef = qd[ini] + qd[j] - two * q_entry(ini, j, &mut cache);
                    let obj_diff = if quad_coef > zero {
                        -(grad_diff * grad_diff) / quad_coef
                    } else {
                        -(grad_diff * grad_diff) / tau
                    };
                    if obj_diff <= obj_diff_min {
                        gmin_idx = j as isize;
                        obj_diff_min = obj_diff;
                    }
                }
            }
        }

        // Stopping criterion (`svm.cpp:1286`).
        if (gmaxp + gmaxp2).max(gmaxn + gmaxn2) < tol || gmin_idx == -1 {
            break;
        }

        let i = if y[gmin_idx as usize] > zero {
            gmaxp_idx
        } else {
            gmaxn_idx
        };
        if i == -1 {
            break;
        }
        let i = i as usize;
        let j = gmin_idx as usize;
        if i == j {
            break;
        }

        // ---- Solver::Solve analytic 2-variable update (svm.cpp:756-852) ----
        let q_ij = q_entry(i, j, &mut cache);
        let c_i = c[i];
        let c_j = c[j];
        let old_alpha_i = alpha[i];
        let old_alpha_j = alpha[j];

        if y[i] != y[j] {
            // `Q_i[j] = y_i·y_j·K = -K(i,j)` here, so `+2·Q_i[j]` in libsvm.
            let mut quad_coef = qd[i] + qd[j] + two * q_ij;
            if quad_coef <= zero {
                quad_coef = tau;
            }
            let delta = (-grad[i] - grad[j]) / quad_coef;
            let diff = alpha[i] - alpha[j];
            alpha[i] = alpha[i] + delta;
            alpha[j] = alpha[j] + delta;

            if diff > zero {
                if alpha[j] < zero {
                    alpha[j] = zero;
                    alpha[i] = diff;
                }
            } else if alpha[i] < zero {
                alpha[i] = zero;
                alpha[j] = -diff;
            }
            if diff > c_i - c_j {
                if alpha[i] > c_i {
                    alpha[i] = c_i;
                    alpha[j] = c_i - diff;
                }
            } else if alpha[j] > c_j {
                alpha[j] = c_j;
                alpha[i] = c_j + diff;
            }
        } else {
            let mut quad_coef = qd[i] + qd[j] - two * q_ij;
            if quad_coef <= zero {
                quad_coef = tau;
            }
            let delta = (grad[i] - grad[j]) / quad_coef;
            let sum = alpha[i] + alpha[j];
            alpha[i] = alpha[i] - delta;
            alpha[j] = alpha[j] + delta;

            if sum > c_i {
                if alpha[i] > c_i {
                    alpha[i] = c_i;
                    alpha[j] = sum - c_i;
                }
            } else if alpha[j] < zero {
                alpha[j] = zero;
                alpha[i] = sum;
            }
            if sum > c_j {
                if alpha[j] > c_j {
                    alpha[j] = c_j;
                    alpha[i] = sum - c_j;
                }
            } else if alpha[i] < zero {
                alpha[i] = zero;
                alpha[j] = sum;
            }
        }

        // ---- Update gradient (svm.cpp:856-862) ----
        let delta_alpha_i = alpha[i] - old_alpha_i;
        let delta_alpha_j = alpha[j] - old_alpha_j;
        #[allow(
            clippy::needless_range_loop,
            reason = "k indexes grad and is also the kernel column Q[k][i]/Q[k][j]"
        )]
        for k in 0..m {
            let q_ki = q_entry(k, i, &mut cache);
            let q_kj = q_entry(k, j, &mut cache);
            grad[k] = grad[k] + q_ki * delta_alpha_i + q_kj * delta_alpha_j;
        }
    }

    // ---- Solver_NU::calculate_rho (svm.cpp:1370-1418) ----
    let mut nr_free1 = 0usize;
    let mut nr_free2 = 0usize;
    let mut ub1 = F::infinity();
    let mut ub2 = F::infinity();
    let mut lb1 = F::neg_infinity();
    let mut lb2 = F::neg_infinity();
    let mut sum_free1 = zero;
    let mut sum_free2 = zero;

    for i in 0..m {
        let upper = is_upper(i, &alpha);
        let lower = is_lower(i, &alpha);
        if y[i] > zero {
            if upper {
                lb1 = lb1.max(grad[i]);
            } else if lower {
                ub1 = ub1.min(grad[i]);
            } else {
                nr_free1 += 1;
                sum_free1 = sum_free1 + grad[i];
            }
        } else if upper {
            lb2 = lb2.max(grad[i]);
        } else if lower {
            ub2 = ub2.min(grad[i]);
        } else {
            nr_free2 += 1;
            sum_free2 = sum_free2 + grad[i];
        }
    }

    let r1 = if nr_free1 > 0 {
        sum_free1 / F::from(nr_free1).unwrap_or_else(F::one)
    } else {
        (ub1 + lb1) / two
    };
    let r2 = if nr_free2 > 0 {
        sum_free2 / F::from(nr_free2).unwrap_or_else(F::one)
    } else {
        (ub2 + lb2) / two
    };

    NuResult {
        alpha,
        rho: (r1 - r2) / two,
        r: (r1 + r2) / two,
    }
}

/// The recovered nu-SVC binary sub-model: support vectors, their `dual_coef`
/// values (`alpha_i·y_i/r`, libsvm's `sv_coef`), original indices, and the
/// libsvm-internal bias `-rho/r` (so the decision function is
/// `f(x) = Σ sv_coef·K(sv, x) - rho/r`).
pub(crate) struct NuSvcModel<F> {
    pub sv_data: Vec<Vec<F>>,
    pub sv_coefs: Vec<F>,
    pub sv_indices: Vec<usize>,
    /// libsvm-internal bias term: `-rho/r`. Used directly as the `+1`-side
    /// decision bias for the LOWER-index class (libsvm convention).
    pub bias_internal: F,
}

/// Solve the libsvm **nu-SVC** dual for a single binary sub-problem
/// (`solve_nu_svc`, `sklearn/svm/src/libsvm/svm.cpp:1646-1708`):
/// `min 0.5 αᵀQα  s.t.  yᵀα=0, eᵀα=ν·l, 0≤α_i≤1`, `Q_ij=y_iy_jK`.
///
/// `data`/`labels` are the per-pair samples and signs (`+1` for the higher-index
/// `class_pos`, `-1` for the lower-index `class_neg`, matching this crate's
/// [`BinarySvm`] convention). The greedy `alpha` init (`svm.cpp:1667-1682`)
/// fills each class up to `min(C_i, nu·l/2 − running_sum)`. After
/// [`solver_nu_core`], libsvm rescales `alpha_i ← alpha_i·y_i/r` and
/// `rho ← rho/r` (`svm.cpp:1696-1702`); the support-vector coefficient is
/// `sv_coef = alpha·y/r` and the decision bias is `−rho/r`.
///
/// Returns `None` when `r ≈ 0` (degenerate — no usable rescale, e.g. a
/// pathological all-bound solution), letting the caller surface a clean error.
#[allow(
    clippy::too_many_arguments,
    reason = "the solver + per-pair samples/labels + hyperparameters thread \
              through one call mirroring libsvm's solve_nu_svc"
)]
pub(crate) fn solve_nu_svc<F: Float, K: Kernel<F>>(
    data: &[Vec<F>],
    labels: &[F],
    kernel: &K,
    nu: F,
    tol: F,
    max_iter: usize,
    cache_size: usize,
) -> Option<NuSvcModel<F>> {
    let l = data.len();
    let zero = F::zero();
    let one = F::one();
    let two = one + one;
    // `C[i] = prob->W[i] = 1` (unit instance weights, `svm.cpp:1664`).
    let c = vec![one; l];

    // Greedy alpha init (`svm.cpp:1667-1682`): `nu_l = Σ nu·C[i] = nu·l`,
    // `sum_pos = sum_neg = nu_l/2`, fill each class greedily up to `C[i]`.
    let nu_l = nu * F::from(l).unwrap_or_else(F::zero);
    let mut sum_pos = nu_l / two;
    let mut sum_neg = nu_l / two;
    let mut alpha = vec![zero; l];
    for i in 0..l {
        if labels[i] > zero {
            alpha[i] = c[i].min(sum_pos);
            sum_pos = sum_pos - alpha[i];
        } else {
            alpha[i] = c[i].min(sum_neg);
            sum_neg = sum_neg - alpha[i];
        }
    }

    let p = vec![zero; l]; // nu-SVC linear term is 0 (`zeros`, svm.cpp:1684).
    let sample = |k: usize| k;
    let res = solver_nu_core(
        data, l, &sample, labels, &p, &c, alpha, kernel, tol, max_iter, cache_size,
    );

    let r = res.r;
    if r.abs() <= F::from(1e-12).unwrap_or_else(F::epsilon) {
        return None;
    }

    // libsvm: `alpha_i *= y_i / r`, `rho /= r` (`svm.cpp:1696-1702`).
    // `sv_coef = alpha·y/r`; decision bias = `-rho/r`.
    let eps_sv = F::from(1e-8).unwrap_or_else(F::epsilon);
    let mut sv_data = Vec::new();
    let mut sv_coefs = Vec::new();
    let mut sv_indices = Vec::new();
    #[allow(
        clippy::needless_range_loop,
        reason = "i indexes alpha, the sign labels[i], and the sample data[i] together"
    )]
    for i in 0..l {
        let coef = res.alpha[i] * labels[i] / r;
        if coef.abs() > eps_sv {
            sv_data.push(data[i].clone());
            sv_coefs.push(coef);
            sv_indices.push(i);
        }
    }

    Some(NuSvcModel {
        sv_data,
        sv_coefs,
        sv_indices,
        bias_internal: (res.rho / r).neg(),
    })
}

/// The recovered nu-SVR model: prediction coefficients `α*−α` per sample,
/// support indices, and the bias.
pub(crate) struct NuSvrModel<F> {
    pub sv_data: Vec<Vec<F>>,
    pub sv_coefs: Vec<F>,
    pub sv_indices: Vec<usize>,
    /// Prediction bias: `f(x) = Σ coef·K(sv, x) + bias`. libsvm's decision
    /// function is `Σ coef·K − rho`, so `bias = −rho`.
    pub bias: F,
}

/// Solve the libsvm **nu-SVR** dual (`solve_nu_svr`,
/// `sklearn/svm/src/libsvm/svm.cpp:1795-1839`): a `2l`-variable
/// `(α, α*)` dual with the learned-tube `nu` constraint, both `nu` AND `C`
/// used (`epsilon` is replaced by `nu`).
///
/// Variable layout (libsvm): `k < l` is `α*_k` with sign `+1` and linear term
/// `−y_k`; `k ≥ l` is `α_{k−l}` with sign `−1` and linear term `+y_{k−l}`;
/// `C_k = W·C` for all (`svm.cpp:1806-1824`). The greedy init fills both halves
/// to `min(sum, C)` where `sum = (Σ C·nu)/2` (`svm.cpp:1814-1817`). After
/// [`solver_nu_core`] the prediction coefficient is `coef_k = α*_k − α_k`
/// (`svm.cpp:1832-1833`) and the bias is `−rho` (libsvm `f = Σ coef·K − rho`).
#[allow(
    clippy::too_many_arguments,
    reason = "the solver + samples/targets + (nu, C) + hyperparameters thread \
              through one call mirroring libsvm's solve_nu_svr"
)]
pub(crate) fn solve_nu_svr<F: Float, K: Kernel<F>>(
    data: &[Vec<F>],
    targets: &[F],
    kernel: &K,
    nu: F,
    c_param: F,
    tol: F,
    max_iter: usize,
    cache_size: usize,
) -> NuSvrModel<F> {
    let l = data.len();
    let m = 2 * l;
    let zero = F::zero();
    let one = F::one();
    let two = one + one;

    // `C[i] = C[i+l] = W·C` (`svm.cpp:1809`); `sum = (Σ C·nu)/2`.
    let c = vec![c_param; m];
    let mut sum = c_param * nu * F::from(l).unwrap_or_else(F::zero) / two;

    let mut alpha = vec![zero; m];
    let mut y = vec![zero; m];
    let mut p = vec![zero; m];
    for i in 0..l {
        let a = sum.min(c[i]); // alpha2[i] = alpha2[i+l] = min(sum, C[i])
        alpha[i] = a;
        alpha[i + l] = a;
        sum = sum - a;
        p[i] = targets[i].neg(); // linear_term[i]   = -y_i
        y[i] = one;
        p[i + l] = targets[i]; // linear_term[i+l] = +y_i
        y[i + l] = -one;
    }

    let sample = |k: usize| if k < l { k } else { k - l };
    let res = solver_nu_core(
        data, m, &sample, &y, &p, &c, alpha, kernel, tol, max_iter, cache_size,
    );

    // coef_i = alpha2[i] - alpha2[i+l] = α*_i - α_i (`svm.cpp:1832-1833`).
    let eps_sv = F::from(1e-8).unwrap_or_else(F::epsilon);
    let mut sv_data = Vec::new();
    let mut sv_coefs = Vec::new();
    let mut sv_indices = Vec::new();
    #[allow(
        clippy::needless_range_loop,
        reason = "i pairs alpha2[i] with alpha2[i+l] and indexes the sample data[i]"
    )]
    for i in 0..l {
        let coef = res.alpha[i] - res.alpha[i + l];
        if coef.abs() > eps_sv {
            sv_data.push(data[i].clone());
            sv_coefs.push(coef);
            sv_indices.push(i);
        }
    }

    NuSvrModel {
        sv_data,
        sv_coefs,
        sv_indices,
        bias: res.rho.neg(),
    }
}

impl<F: Float + Send + Sync + ScalarOperand + 'static, K: Kernel<F> + 'static>
    Fit<Array2<F>, Array1<F>> for SVR<F, K>
{
    type Fitted = FittedSVR<F, K>;
    type Error = FerroError;

    /// Fit the SVR model using SMO.
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if `x` and `y` have different
    /// sample counts.
    /// Returns [`FerroError::InvalidParameter`] if `C` is not positive.
    fn fit(&self, x: &Array2<F>, y: &Array1<F>) -> Result<FittedSVR<F, K>, FerroError> {
        let (n_samples, _n_features) = x.dim();

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

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

        // Reject non-finite input (NaN / +/-inf) in X or the float target y
        // BEFORE the X/y length (`ShapeMismatch`) check, mirroring sklearn's
        // `BaseLibSVM.fit` -> `_validate_data(X, y, ...)`
        // (`sklearn/svm/_base.py:190`) which routes to `check_X_y`:
        // `check_array(X, force_all_finite=True)` then `check_array(y)` both run
        // BEFORE `check_consistent_length(X, y)` (`_base.py:208`), so on an input
        // that is BOTH non-finite AND length-mismatched sklearn raises the
        // finiteness `ValueError("Input X contains NaN.")` /
        // `"Input y contains NaN."` / `"... contains infinity ..."`, NOT a
        // consistency error (#2270). SVR's `y` is float regression targets, so
        // both X and y are finiteness-checked. `.iter().any(|v| !v.is_finite())`
        // catches NaN and +/-inf; on finite input the guard never fires, so the
        // fitted SVR attributes are byte-identical and a finite length-mismatch
        // still yields `ShapeMismatch`.
        if x.iter().any(|v| !v.is_finite()) {
            return Err(FerroError::InvalidParameter {
                name: "X".into(),
                reason: "Input X contains NaN or infinity.".into(),
            });
        }
        if y.iter().any(|v| !v.is_finite()) {
            return Err(FerroError::InvalidParameter {
                name: "y".into(),
                reason: "Input y contains NaN or infinity.".into(),
            });
        }

        if 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(),
            });
        }

        // Resolve any data-dependent kernel parameters (e.g. a `None` gamma ->
        // sklearn's default `gamma='scale'` = 1/(n_features * X.var()),
        // `_base.py:236-239`) against the training data BEFORE fitting, and use
        // this resolved kernel for both fitting and prediction.
        let kernel = self.kernel.resolved_for_fit(x);

        let data: Vec<Vec<F>> = (0..n_samples).map(|i| x.row(i).to_vec()).collect();
        let targets: Vec<F> = y.to_vec();

        let (coefs, bias) = smo_svr(
            &data,
            &targets,
            &kernel,
            self.c,
            self.epsilon,
            self.tol,
            self.max_iter,
            self.cache_size,
        )?;

        // Extract support vectors (non-zero coefficients).
        let eps = F::from(1e-8).unwrap_or_else(F::epsilon);
        let mut sv_data = Vec::new();
        let mut sv_coefs = Vec::new();
        let mut sv_idx = Vec::new();

        for (i, &coef) in coefs.iter().enumerate() {
            if coef.abs() > eps {
                sv_data.push(data[i].clone());
                sv_coefs.push(coef);
                sv_idx.push(i);
            }
        }

        Ok(FittedSVR {
            kernel,
            support_vectors: sv_data,
            sv_indices: sv_idx,
            dual_coefs: sv_coefs,
            bias,
        })
    }
}

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

    /// Predict target values for the given feature matrix.
    ///
    /// # Errors
    ///
    /// Returns `Ok` always for valid input.
    fn predict(&self, x: &Array2<F>) -> Result<Array1<F>, FerroError> {
        self.decision_function(x)
    }
}

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

    #[test]
    fn test_linear_kernel() {
        let k = LinearKernel;
        let x = vec![1.0, 2.0, 3.0];
        let y = vec![4.0, 5.0, 6.0];
        assert_relative_eq!(k.compute(&x, &y), 32.0, epsilon = 1e-10);
    }

    #[test]
    fn test_rbf_kernel() {
        let k = RbfKernel::with_gamma(1.0);
        let x = vec![0.0, 0.0];
        let y = vec![0.0, 0.0];
        assert_relative_eq!(k.compute(&x, &y), 1.0, epsilon = 1e-10);

        // Different points should give value < 1
        let y2 = vec![1.0, 0.0];
        let val: f64 = k.compute(&x, &y2);
        assert!(val < 1.0);
        assert!(val > 0.0);
    }

    #[test]
    fn test_polynomial_kernel() {
        let k = PolynomialKernel {
            gamma: Gamma::Value(1.0),
            degree: 2,
            coef0: 1.0,
        };
        let x = vec![1.0f64, 0.0];
        let y = vec![1.0, 0.0];
        // (1.0 * 1.0 + 1.0)^2 = 4.0
        assert_relative_eq!(k.compute(&x, &y), 4.0, epsilon = 1e-10);
    }

    #[test]
    fn test_sigmoid_kernel() {
        let k = SigmoidKernel {
            gamma: Gamma::Value(1.0),
            coef0: 0.0,
        };
        let x = vec![0.0f64];
        let y = vec![0.0];
        // tanh(0) = 0
        assert_relative_eq!(k.compute(&x, &y), 0.0, epsilon = 1e-10);
    }

    #[test]
    fn test_svc_linear_separable() {
        // Two well-separated clusters.
        let x = Array2::from_shape_vec(
            (8, 2),
            vec![
                1.0, 1.0, 1.5, 1.0, 1.0, 1.5, 1.5, 1.5, // class 0
                5.0, 5.0, 5.5, 5.0, 5.0, 5.5, 5.5, 5.5, // class 1
            ],
        )
        .unwrap();
        let y = array![0usize, 0, 0, 0, 1, 1, 1, 1];

        let model = SVC::<f64, LinearKernel>::new(LinearKernel).with_c(10.0);
        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_svc_rbf_xor() {
        // XOR problem: not linearly separable, needs RBF kernel.
        let x = Array2::from_shape_vec(
            (8, 2),
            vec![
                0.0, 0.0, 0.1, 0.1, // class 0
                1.0, 1.0, 1.1, 1.1, // class 0
                1.0, 0.0, 1.1, 0.1, // class 1
                0.0, 1.0, 0.1, 1.1, // class 1
            ],
        )
        .unwrap();
        let y = array![0usize, 0, 0, 0, 1, 1, 1, 1];

        let kernel = RbfKernel::with_gamma(5.0);
        let model = SVC::new(kernel).with_c(100.0).with_max_iter(50000);
        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 for XOR, got {correct}"
        );
    }

    #[test]
    fn test_svc_multiclass() {
        // Three linearly separable clusters.
        let x = Array2::from_shape_vec(
            (9, 2),
            vec![
                0.0, 0.0, 0.5, 0.0, 0.0, 0.5, // class 0
                5.0, 0.0, 5.5, 0.0, 5.0, 0.5, // class 1
                0.0, 5.0, 0.5, 5.0, 0.0, 5.5, // class 2
            ],
        )
        .unwrap();
        let y = array![0usize, 0, 0, 1, 1, 1, 2, 2, 2];

        let model = SVC::new(LinearKernel).with_c(10.0);
        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 >= 7,
            "Expected at least 7 correct for multiclass, got {correct}"
        );
    }

    #[test]
    fn test_svc_decision_function() -> TestResult {
        let x = Array2::from_shape_vec(
            (6, 2),
            vec![
                1.0, 1.0, 1.5, 1.0, 1.0, 1.5, // class 0
                5.0, 5.0, 5.5, 5.0, 5.0, 5.5, // class 1
            ],
        )
        .map_err(|_| err("shape"))?;
        let y = array![0usize, 0, 0, 1, 1, 1];

        let model = SVC::new(LinearKernel).with_c(10.0);
        let fitted = model.fit(&x, &y)?;

        let df = fitted.decision_function(&x)?;
        // Binary: 1-D (n,) score (sklearn ravels the single ovo column,
        // `_base.py:538-539`); the multiclass borrow is None.
        assert!(df.as_multiclass().is_none());
        let bin = df.as_binary().ok_or_else(|| err("binary 1-D"))?;
        assert_eq!(bin.len(), 6);

        // Class 0 samples should have negative decision values,
        // class 1 should have positive (positive -> classes_[1]).
        for (i, &v) in bin.iter().enumerate().take(3) {
            assert!(v < 0.5, "Class 0 sample {i} has decision value {v}");
        }
        Ok(())
    }

    // -----------------------------------------------------------------------
    // REQ-4 smoke tests: decision_function shape + sign + ovr/ovo transform.
    //
    // Expected values from the LIVE sklearn 1.5.2 oracle (R-CHAR-3):
    //
    //   import numpy as np; from sklearn.svm import SVC
    //   X=np.array([[1.,1.],[2.,1.],[1.,2.],[5.,5.],[6.,5.],[5.,6.]])
    //   y=np.array([0,0,0,1,1,1]); m=SVC(kernel='linear',C=1.0).fit(X,y)
    //   m.decision_function(X)  # (6,) [-1.2853,-0.9997,-0.9997,0.9995,1.2851,1.2851]
    //
    //   X3=[[0,0],[.5,0],[0,.5],[5,0],[5.5,0],[5,.5],[0,5],[.5,5],[0,5.5]]
    //   y3=[0,0,0,1,1,1,2,2,2]; m3=SVC(kernel='linear',C=1.0).fit(X3,y3)
    //   m3.decision_function(X3)            # ovr (9,3) row0 [2.2366,0.8167,-0.1833]
    //                                       #         row3 [1.0606,2.2262,-0.2333]
    //   SVC(...,decision_function_shape='ovo').fit(X3,y3).decision_function(X3)
    //                                       # ovo (9,3) row0 [1.2222,1.2222,0.0]
    // -----------------------------------------------------------------------

    fn three_class_9x2() -> Result<(Array2<f64>, Array1<usize>), FerroError> {
        let x = Array2::from_shape_vec(
            (9, 2),
            vec![
                0.0, 0.0, 0.5, 0.0, 0.0, 0.5, 5.0, 0.0, 5.5, 0.0, 5.0, 0.5, 0.0, 5.0, 0.5, 5.0,
                0.0, 5.5,
            ],
        )
        .map_err(|_| err("shape"))?;
        let y = array![0usize, 0, 0, 1, 1, 1, 2, 2, 2];
        Ok((x, y))
    }

    #[test]
    fn test_svc_decision_function_binary_values() -> TestResult {
        let m = binary_fit()?;
        let x = Array2::from_shape_vec(
            (6, 2),
            vec![1.0, 1.0, 2.0, 1.0, 1.0, 2.0, 5.0, 5.0, 6.0, 5.0, 5.0, 6.0],
        )
        .map_err(|_| err("shape"))?;
        let df = m.decision_function(&x)?;
        let bin = df.as_binary().ok_or_else(|| err("binary"))?;
        assert_eq!(bin.len(), 6);
        let oracle = [-1.2853, -0.9997, -0.9997, 0.9995, 1.2851, 1.2851];
        for (i, &exp) in oracle.iter().enumerate() {
            assert!(
                (bin[i] - exp).abs() < 1e-2,
                "binary df[{i}] = {} vs oracle {exp}",
                bin[i]
            );
        }
        Ok(())
    }

    #[test]
    fn test_svc_decision_function_ovr() -> TestResult {
        let (x, y) = three_class_9x2()?;
        let m = SVC::new(LinearKernel)
            .with_c(1.0)
            .with_tol(1e-6)
            .with_max_iter(200_000)
            .fit(&x, &y)?;
        let df = m.decision_function(&x)?;
        let mc = df.as_multiclass().ok_or_else(|| err("multiclass"))?;
        assert_eq!(mc.dim(), (9, 3));
        // ovr (default): row0 [2.2366,0.8167,-0.1833], row3 [1.0606,2.2262,-0.2333].
        let row0 = [2.2366, 0.8167, -0.1833];
        let row3 = [1.0606, 2.2262, -0.2333];
        for (c, &v) in row0.iter().enumerate() {
            assert!(
                (mc[[0, c]] - v).abs() < 1e-2,
                "ovr row0[{c}] = {} vs oracle {v}",
                mc[[0, c]]
            );
        }
        for (c, &v) in row3.iter().enumerate() {
            assert!(
                (mc[[3, c]] - v).abs() < 1e-2,
                "ovr row3[{c}] = {} vs oracle {v}",
                mc[[3, c]]
            );
        }
        Ok(())
    }

    #[test]
    fn test_svc_decision_function_ovo() -> TestResult {
        let (x, y) = three_class_9x2()?;
        let m = SVC::new(LinearKernel)
            .with_c(1.0)
            .with_tol(1e-6)
            .with_max_iter(200_000)
            .with_decision_function_shape(SvmDecisionShape::Ovo)
            .fit(&x, &y)?;
        let df = m.decision_function(&x)?;
        let mc = df.as_multiclass().ok_or_else(|| err("multiclass"))?;
        assert_eq!(mc.dim(), (9, 3));
        // ovo: row0 [1.2222,1.2222,0.0].
        let row0 = [1.2222, 1.2222, 0.0];
        for (c, &v) in row0.iter().enumerate() {
            assert!(
                (mc[[0, c]] - v).abs() < 1e-2,
                "ovo row0[{c}] = {} vs oracle {v}",
                mc[[0, c]]
            );
        }
        Ok(())
    }

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

        let model = SVC::new(LinearKernel).with_c(0.0);
        assert!(model.fit(&x, &y).is_err());
    }

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

        let model = SVC::new(LinearKernel);
        assert!(model.fit(&x, &y).is_err());
    }

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

        let model = SVC::new(LinearKernel);
        assert!(model.fit(&x, &y).is_err());
    }

    #[test]
    fn test_svr_simple() {
        // Simple linear regression: y = 2x
        let x = Array2::from_shape_vec((6, 1), vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]).unwrap();
        let y = array![2.0, 4.0, 6.0, 8.0, 10.0, 12.0];

        let model = SVR::new(LinearKernel)
            .with_c(100.0)
            .with_epsilon(0.1)
            .with_max_iter(50000);
        let fitted = model.fit(&x, &y).unwrap();
        let preds = fitted.predict(&x).unwrap();

        // Check predictions are reasonably close.
        for (p, &actual) in preds.iter().zip(y.iter()) {
            assert!(
                (*p - actual).abs() < 2.0,
                "SVR prediction {p} too far from actual {actual}"
            );
        }
    }

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

        let model = SVR::new(LinearKernel).with_c(100.0).with_epsilon(0.1);
        let fitted = model.fit(&x, &y).unwrap();

        let df = fitted.decision_function(&x).unwrap();
        let preds = fitted.predict(&x).unwrap();

        // Decision function and predict should return the same values.
        for i in 0..4 {
            assert_relative_eq!(df[i], preds[i], epsilon = 1e-10);
        }
    }

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

        let model = SVR::new(LinearKernel).with_c(-1.0);
        assert!(model.fit(&x, &y).is_err());
    }

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

        let model = SVR::new(LinearKernel);
        assert!(model.fit(&x, &y).is_err());
    }

    // -----------------------------------------------------------------------
    // REQ-3 smoke tests: libsvm-layout fitted attributes (binary sign flip).
    //
    // Expected values from the LIVE sklearn 1.5.2 oracle (R-CHAR-3), never
    // copied from the ferrolearn side:
    //
    //   python3 -c "import numpy as np; from sklearn.svm import SVC, SVR
    //   X=np.array([[1.,1.],[2.,1.],[1.,2.],[5.,5.],[6.,5.],[5.,6.]])
    //   y=np.array([0,0,0,1,1,1]); m=SVC(kernel='linear',C=1.0).fit(X,y)
    //   print(m.support_.tolist(), m.n_support_.tolist(),
    //         m.dual_coef_.tolist(), m.intercept_.tolist(), m.coef_.tolist())"
    //   # [1, 2, 3] [2, 1] [[-0.0408,-0.0408,0.0816]] [-1.8565] [[0.2856,0.2856]]
    //
    //   X3=[[0,0],[.5,0],[0,.5],[5,0],[5.5,0],[5,.5],[0,5],[.5,5],[0,5.5]]
    //   y3=[0,0,0,1,1,1,2,2,2]; m3=SVC(kernel='linear',C=1.0).fit(X3,y3)
    //   # support_ [1,2,3,5,6,7] n_support_ [2,2,2]
    //   # dual_coef_ [[0.0988,0,-0.0988,0,-0.0988,0],[0,0.0988,0,0.0494,0,-0.0494]]
    //   # intercept_ [1.2222,1.2222,0.0]
    //
    //   Xr=[[1],[2],[3],[4],[5],[6]]; yr=[2,4,6,8,10,12]
    //   mr=SVR(kernel='linear',C=100,epsilon=0.1).fit(Xr,yr)
    //   # support_ [0,5] dual_coef_ [[-0.392,0.392]] intercept_ [0.14] n_support_ [2]
    //
    // The tests return `Result` and use `?`/`ok_or` (no unwrap/expect/panic).
    // -----------------------------------------------------------------------

    type TestResult = Result<(), FerroError>;

    fn err(msg: &str) -> FerroError {
        FerroError::InvalidParameter {
            name: "test".into(),
            reason: msg.into(),
        }
    }

    fn binary_fit() -> Result<FittedSVC<f64, LinearKernel>, FerroError> {
        let x = Array2::from_shape_vec(
            (6, 2),
            vec![1.0, 1.0, 2.0, 1.0, 1.0, 2.0, 5.0, 5.0, 6.0, 5.0, 5.0, 6.0],
        )
        .map_err(|_| err("shape"))?;
        let y = array![0usize, 0, 0, 1, 1, 1];
        SVC::new(LinearKernel)
            .with_c(1.0)
            .with_tol(1e-6)
            .with_max_iter(200_000)
            .fit(&x, &y)
    }

    #[test]
    fn test_svc_binary_support_attrs() -> TestResult {
        let m = binary_fit()?;
        // support_ [1,2,3], grouped by class (class0:[1,2], class1:[3]).
        assert_eq!(m.support().to_vec(), vec![1, 2, 3]);
        // n_support_ [2,1].
        assert_eq!(m.n_support(), vec![2, 1]);
        // support_vectors_ = X[support_].
        let svs = m.support_vectors();
        assert_eq!(svs.dim(), (3, 2));
        let expected = [[2.0, 1.0], [1.0, 2.0], [5.0, 5.0]];
        for (r, row) in expected.iter().enumerate() {
            for (c, &v) in row.iter().enumerate() {
                assert_relative_eq!(svs[[r, c]], v, epsilon = 1e-10);
            }
        }
        Ok(())
    }

    #[test]
    fn test_svc_binary_dual_coef_sign_flip() -> TestResult {
        let m = binary_fit()?;
        // dual_coef_ shape (1,3) = [[-0.0408,-0.0408,0.0816]] (binary sign flip).
        let dc = m.dual_coef();
        assert_eq!(dc.dim(), (1, 3));
        let oracle = [-0.0408, -0.0408, 0.0816];
        for (c, &v) in oracle.iter().enumerate() {
            assert!(
                (dc[[0, c]] - v).abs() < 1e-2,
                "dual_coef_[0,{c}] = {} vs oracle {v}",
                dc[[0, c]]
            );
        }
        Ok(())
    }

    #[test]
    fn test_svc_binary_intercept_and_coef() -> TestResult {
        let m = binary_fit()?;
        // intercept_ [-1.8565], length 1 (binary sign flip).
        let ic = m.intercept();
        assert_eq!(ic.len(), 1);
        assert!(
            (ic[0] - (-1.8565)).abs() < 1e-2,
            "intercept_ = {} vs oracle -1.8565",
            ic[0]
        );
        // coef_ [[0.2856,0.2856]] shape (1,2) for the linear kernel.
        let coef = m.coef().ok_or_else(|| err("linear kernel exposes coef_"))?;
        assert_eq!(coef.dim(), (1, 2));
        for c in 0..2 {
            assert!(
                (coef[[0, c]] - 0.2856).abs() < 1e-2,
                "coef_[0,{c}] = {} vs oracle 0.2856",
                coef[[0, c]]
            );
        }
        Ok(())
    }

    #[test]
    fn test_svc_coef_none_for_nonlinear() -> TestResult {
        // coef_ is only available for the linear kernel; RBF -> None
        // (sklearn raises AttributeError, _base.py:650-651).
        let x = Array2::from_shape_vec(
            (6, 2),
            vec![1.0, 1.0, 2.0, 1.0, 1.0, 2.0, 5.0, 5.0, 6.0, 5.0, 5.0, 6.0],
        )
        .map_err(|_| err("shape"))?;
        let y = array![0usize, 0, 0, 1, 1, 1];
        let m = SVC::new(RbfKernel::with_gamma(0.5)).fit(&x, &y)?;
        assert!(m.coef().is_none());
        Ok(())
    }

    fn multiclass_fit() -> Result<FittedSVC<f64, LinearKernel>, FerroError> {
        let x = Array2::from_shape_vec(
            (9, 2),
            vec![
                0.0, 0.0, 0.5, 0.0, 0.0, 0.5, 5.0, 0.0, 5.5, 0.0, 5.0, 0.5, 0.0, 5.0, 0.5, 5.0,
                0.0, 5.5,
            ],
        )
        .map_err(|_| err("shape"))?;
        let y = array![0usize, 0, 0, 1, 1, 1, 2, 2, 2];
        SVC::new(LinearKernel)
            .with_c(1.0)
            .with_tol(1e-6)
            .with_max_iter(200_000)
            .fit(&x, &y)
    }

    #[test]
    fn test_svc_multiclass_support_attrs() -> TestResult {
        let m = multiclass_fit()?;
        // support_ [1,2,3,5,6,7] grouped by class; n_support_ [2,2,2].
        assert_eq!(m.support().to_vec(), vec![1, 2, 3, 5, 6, 7]);
        assert_eq!(m.n_support(), vec![2, 2, 2]);
        // intercept_ [1.2222,1.2222,0.0] (no sign flip for multiclass).
        let ic = m.intercept();
        assert_eq!(ic.len(), 3);
        let oracle_ic = [1.2222, 1.2222, 0.0];
        for (i, &v) in oracle_ic.iter().enumerate() {
            assert!(
                (ic[i] - v).abs() < 1e-2,
                "intercept_[{i}] = {} vs oracle {v}",
                ic[i]
            );
        }
        Ok(())
    }

    #[test]
    fn test_svc_multiclass_dual_coef_packing() -> TestResult {
        let m = multiclass_fit()?;
        // dual_coef_ shape (2,6), libsvm packing (cols = SVs [1,2,3,5,6,7]):
        //   row0 = [0.0988, 0.0, -0.0988, 0.0, -0.0988, 0.0]
        //   row1 = [0.0, 0.0988, 0.0, 0.0494, 0.0, -0.0494]
        let dc = m.dual_coef();
        assert_eq!(dc.dim(), (2, 6));
        let oracle = [
            [0.0988, 0.0, -0.0988, 0.0, -0.0988, 0.0],
            [0.0, 0.0988, 0.0, 0.0494, 0.0, -0.0494],
        ];
        for (r, row) in oracle.iter().enumerate() {
            for (c, &v) in row.iter().enumerate() {
                assert!(
                    (dc[[r, c]] - v).abs() < 1e-2,
                    "dual_coef_[{r},{c}] = {} vs oracle {v}",
                    dc[[r, c]]
                );
            }
        }
        Ok(())
    }

    #[test]
    fn test_svr_linear_attrs() -> TestResult {
        // SVR(kernel='linear', C=100, epsilon=0.1) on the 6x1 set.
        let x = Array2::from_shape_vec((6, 1), vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0])
            .map_err(|_| err("shape"))?;
        let y = array![2.0, 4.0, 6.0, 8.0, 10.0, 12.0];
        let m = SVR::new(LinearKernel)
            .with_c(100.0)
            .with_epsilon(0.1)
            .with_tol(1e-6)
            .with_max_iter(200_000)
            .fit(&x, &y)?;
        // support_ [0,5]; n_support_ [2]; dual_coef_ (1,2) [[-0.392,0.392]];
        // intercept_ [0.14].
        assert_eq!(m.support().to_vec(), vec![0, 5]);
        assert_eq!(m.n_support(), vec![2]);
        let dc = m.dual_coef();
        assert_eq!(dc.dim(), (1, 2));
        assert!(
            (dc[[0, 0]] - (-0.392)).abs() < 1e-2,
            "dual_coef_[0,0] = {} vs oracle -0.392",
            dc[[0, 0]]
        );
        assert!(
            (dc[[0, 1]] - 0.392).abs() < 1e-2,
            "dual_coef_[0,1] = {} vs oracle 0.392",
            dc[[0, 1]]
        );
        let ic = m.intercept();
        assert_eq!(ic.len(), 1);
        assert!(
            (ic[0] - 0.14).abs() < 1e-2,
            "intercept_ = {} vs oracle 0.14",
            ic[0]
        );
        Ok(())
    }

    // -----------------------------------------------------------------------
    // REQ-1 (gamma='auto') + REQ-8 (break_ties) smoke tests.
    //
    // Expected values from the LIVE sklearn 1.5.2 oracle (R-CHAR-3):
    //
    //   import numpy as np; from sklearn.svm import SVC
    //   X=np.array([[1.,1.],[2.,1.],[1.,2.],[5.,5.],[6.,5.],[5.,6.]])
    //   y=np.array([0,0,0,1,1,1])
    //   m=SVC(kernel='rbf',C=1.0,gamma='auto').fit(X,y)
    //   m._gamma                       # 0.5  (= 1/n_features = 1/2)
    //   m.decision_function(X)         # [-0.9996,-0.9999,-0.9999,
    //                                  #   0.9999, 0.9999, 0.9996]
    //
    //   break_ties: a symmetric, cleanly-separable 3-class set so ferrolearn's
    //   SMO converges to libsvm's optimum (each class has 2 SVs). Near the
    //   centroid the three ovr scores are ~1.0 (a near 1-1-1 vote tie), so the
    //   libsvm vote breaks toward the LOWEST class index (0) while ovr-argmax
    //   breaks by confidence. Oracle (re-derived vs the live oracle):
    //   Xb=[[0,0],[.5,0],[0,.5],[10,0],[10.5,0],[10,.5],[5,9],[5.5,9],[5,9.5]]
    //   yb=[0,0,0,1,1,1,2,2,2]
    //   q1=[[5.19,3.342]]: vote -> 0, break_ties -> 2 (df [0.9942,0.999,1.0068])
    //   q2=[[5.19,3.241]]: vote -> 0, break_ties -> 1 (df [0.9999,1.0051,0.9949])
    //   SVC(...,decision_function_shape='ovo',break_ties=True) -> ValueError
    // -----------------------------------------------------------------------

    #[test]
    fn test_svc_gamma_auto_decision_function() -> TestResult {
        // gamma='auto' on the 6x2 set: _gamma = 1/n_features = 0.5.
        let x = Array2::from_shape_vec(
            (6, 2),
            vec![1.0, 1.0, 2.0, 1.0, 1.0, 2.0, 5.0, 5.0, 6.0, 5.0, 5.0, 6.0],
        )
        .map_err(|_| err("shape"))?;
        let y = array![0usize, 0, 0, 1, 1, 1];

        // Confirm the resolved gamma is 1/n_features (sklearn `_base.py:241`).
        let resolved = RbfKernel::<f64>::with_gamma_auto().resolved_for_fit(&x);
        assert!(
            (gamma_value_or_one(resolved.gamma) - 0.5).abs() < 1e-12,
            "gamma='auto' resolved to {} vs oracle 0.5",
            gamma_value_or_one(resolved.gamma)
        );

        let m = SVC::new(RbfKernel::<f64>::with_gamma_auto())
            .with_c(1.0)
            .with_tol(1e-6)
            .with_max_iter(200_000)
            .fit(&x, &y)?;
        let df = m.decision_function(&x)?;
        let bin = df.as_binary().ok_or_else(|| err("binary"))?;
        assert_eq!(bin.len(), 6);
        let oracle = [-0.9996, -0.9999, -0.9999, 0.9999, 0.9999, 0.9996];
        for (i, &exp) in oracle.iter().enumerate() {
            assert!(
                (bin[i] - exp).abs() < 1e-2,
                "gamma=auto df[{i}] = {} vs oracle {exp}",
                bin[i]
            );
        }
        Ok(())
    }

    #[test]
    fn test_svc_gamma_scale_still_default() -> TestResult {
        // gamma='scale' (the default) must STILL resolve to
        // 1/(n_features * X.var()); on the 6x2 set X.var()=4.2222 ->
        // _gamma = 1/(2*4.2222) = 0.11842 (sklearn `_base.py:238-239`).
        let x = Array2::from_shape_vec(
            (6, 2),
            vec![1.0, 1.0, 2.0, 1.0, 1.0, 2.0, 5.0, 5.0, 6.0, 5.0, 5.0, 6.0],
        )
        .map_err(|_| err("shape"))?;
        let resolved = RbfKernel::<f64>::new().resolved_for_fit(&x);
        assert!(
            (gamma_value_or_one(resolved.gamma) - 0.118_421).abs() < 1e-4,
            "gamma='scale' resolved to {} vs oracle 0.118421",
            gamma_value_or_one(resolved.gamma)
        );
        Ok(())
    }

    fn break_ties_set() -> Result<(Array2<f64>, Array1<usize>), FerroError> {
        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, 5.0, 9.0, 5.5, 9.0,
                5.0, 9.5,
            ],
        )
        .map_err(|_| err("shape"))?;
        let y = array![0usize, 0, 0, 1, 1, 1, 2, 2, 2];
        Ok((x, y))
    }

    fn break_ties_fit(
        break_ties: bool,
        shape: SvmDecisionShape,
    ) -> Result<FittedSVC<f64, LinearKernel>, FerroError> {
        let (x, y) = break_ties_set()?;
        SVC::new(LinearKernel)
            .with_c(1.0)
            .with_tol(1e-6)
            .with_max_iter(200_000)
            .with_break_ties(break_ties)
            .with_decision_function_shape(shape)
            .fit(&x, &y)
    }

    #[test]
    fn test_svc_break_ties_changes_label() -> TestResult {
        // q1=(5.19,3.342): vote -> 0, break_ties -> 2.
        let q1 = Array2::from_shape_vec((1, 2), vec![5.19, 3.342]).map_err(|_| err("shape"))?;
        // q2=(5.19,3.241): vote -> 0, break_ties -> 1.
        let q2 = Array2::from_shape_vec((1, 2), vec![5.19, 3.241]).map_err(|_| err("shape"))?;

        let vote = break_ties_fit(false, SvmDecisionShape::Ovr)?;
        let bt = break_ties_fit(true, SvmDecisionShape::Ovr)?;

        // break_ties=false (default): libsvm vote -> lowest-index class 0.
        assert_eq!(vote.predict(&q1)?[0], 0, "vote q1 should be 0");
        assert_eq!(vote.predict(&q2)?[0], 0, "vote q2 should be 0");

        // break_ties=true + ovr: argmax(decision_function).
        assert_eq!(
            bt.predict(&q1)?[0],
            2,
            "break_ties q1 (ovr-argmax) should be 2"
        );
        assert_eq!(
            bt.predict(&q2)?[0],
            1,
            "break_ties q2 (ovr-argmax) should be 1"
        );
        Ok(())
    }

    #[test]
    fn test_svc_break_ties_ovo_errors() -> TestResult {
        // sklearn raises when break_ties=True and decision_function_shape='ovo'
        // (`_base.py:801-804`).
        let m = break_ties_fit(true, SvmDecisionShape::Ovo)?;
        let q = Array2::from_shape_vec((1, 2), vec![5.19, 3.342]).map_err(|_| err("shape"))?;
        assert!(m.predict(&q).is_err());
        Ok(())
    }

    #[test]
    fn test_svc_default_params() {
        // sklearn defaults: cache_size=200, max_iter=-1 (= 0 sentinel),
        // shrinking=True, break_ties=False, decision_function_shape='ovr'.
        let m = SVC::<f64, LinearKernel>::new(LinearKernel);
        assert_eq!(m.cache_size, 200);
        assert_eq!(m.max_iter, 0);
        assert!(m.shrinking);
        assert!(!m.break_ties);
        assert_eq!(m.decision_function_shape, SvmDecisionShape::Ovr);
        let r = SVR::<f64, LinearKernel>::new(LinearKernel);
        assert_eq!(r.cache_size, 200);
        assert_eq!(r.max_iter, 0);
        assert!(r.shrinking);
    }

    /// The overlapping imbalanced binary set used to pin `class_weight`.
    fn class_weight_xy() -> Result<(Array2<f64>, Array1<usize>), FerroError> {
        let x = Array2::from_shape_vec(
            (8, 2),
            vec![
                0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.5, 0.5, 1.5, 0.5, 2.0, 2.0, 2.5, 2.5,
            ],
        )
        .map_err(|_| err("shape"))?;
        let y = array![0usize, 0, 0, 0, 0, 1, 1, 1];
        Ok((x, y))
    }

    fn cw_fit(cw: ClassWeight<f64>) -> Result<FittedSVC<f64, LinearKernel>, FerroError> {
        let (x, y) = class_weight_xy()?;
        SVC::new(LinearKernel)
            .with_c(1.0)
            .with_tol(1e-7)
            .with_max_iter(500_000)
            .with_class_weight(cw)
            .fit(&x, &y)
    }

    /// `class_weight` per-class C in the C-SVC SMO (REQ-8, #641).
    ///
    /// Oracle (live `SVC(kernel='linear', C=1.0, class_weight=...)` on the
    /// overlapping imbalanced binary set, R-CHAR-3):
    /// ```text
    /// X=[[0,0],[1,0],[0,1],[1,1],[0.5,0.5],[1.5,0.5],[2,2],[2.5,2.5]] y=[0,0,0,0,0,1,1,1]
    /// None     -> dual_coef_ [[-0.5,-1.0,1.0,0.5]]      intercept_ [-2.0]    support_ [1,3,5,6]
    /// balanced -> dual_coef_ [[-0.8,-0.8,1.3333,0.2667]] intercept_ [-1.6667] support_ [1,3,5,6]
    /// {0:1,1:5}-> support_ [1,3,4,5] intercept_ [-2.0]
    /// ```
    #[test]
    fn test_svc_class_weight_smoke() -> TestResult {
        // class_weight=None.
        let none = cw_fit(ClassWeight::None)?;
        assert_eq!(none.support().to_vec(), vec![1, 3, 5, 6]);
        let dc = none.dual_coef();
        for (c, &v) in [-0.5, -1.0, 1.0, 0.5].iter().enumerate() {
            assert!(
                (dc[[0, c]] - v).abs() < 1e-2,
                "None dual_coef_[0,{c}] = {} vs {v}",
                dc[[0, c]]
            );
        }
        let none_int = none.intercept()[0];
        assert!(
            (none_int - (-2.0)).abs() < 1e-2,
            "None intercept_ {none_int}"
        );

        // class_weight='balanced' (weights [0.8, 1.3333]).
        let bal = cw_fit(ClassWeight::Balanced)?;
        assert_eq!(bal.support().to_vec(), vec![1, 3, 5, 6]);
        let dcb = bal.dual_coef();
        for (c, &v) in [-0.8, -0.8, 1.3333, 0.2667].iter().enumerate() {
            assert!(
                (dcb[[0, c]] - v).abs() < 1e-2,
                "balanced dual_coef_[0,{c}] = {} vs {v}",
                dcb[[0, c]]
            );
        }
        let bal_int = bal.intercept()[0];
        assert!(
            (bal_int - (-1.6667)).abs() < 1e-2,
            "balanced intercept_ {bal_int}"
        );

        // class_weight={0:1, 1:5}.
        let exp = cw_fit(ClassWeight::Explicit(vec![(0, 1.0), (1, 5.0)]))?;
        assert_eq!(exp.support().to_vec(), vec![1, 3, 4, 5]);
        let exp_int = exp.intercept()[0];
        assert!(
            (exp_int - (-2.0)).abs() < 1e-2,
            "explicit intercept_ {exp_int}"
        );

        // None vs balanced MUST give different intercepts — fails if
        // class_weight were ignored (R-CHAR-1).
        assert!(
            (none_int - bal_int).abs() > 1e-2,
            "None intercept {none_int} must differ from balanced {bal_int}"
        );
        Ok(())
    }

    // -----------------------------------------------------------------------
    // REQ-9 smoke tests: probability / predict_proba (Platt scaling).
    //
    // These pin the DETERMINISTIC contract + STRUCTURAL invariants only, NOT
    // exact probA/probB or predict_proba values. sklearn's predict_proba is
    // RNG-CV-dependent (probA_ = -0.7749 at random_state=0 vs -1.0541 at
    // random_state=1 on the binary set), so exact values are NOT a stable
    // oracle (R-CHAR-3: the asserted invariants are sklearn's DOCUMENTED
    // contract — `_base.py:829-864` "columns correspond to classes_ in sorted
    // order", `predict_proba` rows are a probability distribution — never
    // copied from the ferrolearn side).
    // -----------------------------------------------------------------------

    #[test]
    fn test_svc_predict_proba_raises_when_probability_false() -> TestResult {
        // probability=false (default): predict_proba/predict_log_proba error,
        // mirroring sklearn's raise (`_base.py:820-827`/`856-860`).
        let m = binary_fit()?; // default probability=false
        let x = Array2::from_shape_vec((1, 2), vec![3.0, 3.0]).map_err(|_| err("shape"))?;
        assert!(m.predict_proba(&x).is_err());
        assert!(m.predict_log_proba(&x).is_err());
        Ok(())
    }

    fn proba_binary_fit() -> Result<FittedSVC<f64, LinearKernel>, FerroError> {
        let x = Array2::from_shape_vec(
            (10, 2),
            vec![
                1.0, 1.0, 2.0, 1.0, 1.0, 2.0, 2.0, 2.0, 1.5, 1.5, 5.0, 5.0, 6.0, 5.0, 5.0, 6.0,
                6.0, 6.0, 5.5, 5.5,
            ],
        )
        .map_err(|_| err("shape"))?;
        let y = array![0usize, 0, 0, 0, 0, 1, 1, 1, 1, 1];
        SVC::new(LinearKernel)
            .with_c(1.0)
            .with_tol(1e-6)
            .with_max_iter(200_000)
            .with_probability(true)
            .fit(&x, &y)
    }

    #[test]
    fn test_svc_predict_proba_binary_rows_sum_to_one() -> TestResult {
        let m = proba_binary_fit()?;
        let x = Array2::from_shape_vec((4, 2), vec![1.0, 1.0, 1.5, 1.5, 5.0, 5.0, 5.5, 5.5])
            .map_err(|_| err("shape"))?;
        let p = m.predict_proba(&x)?;
        assert_eq!(p.dim(), (4, 2));
        for s in 0..4 {
            let row_sum = p[[s, 0]] + p[[s, 1]];
            assert!((row_sum - 1.0).abs() < 1e-9, "row {s} sums to {row_sum}");
            for c in 0..2 {
                assert!(
                    p[[s, c]] >= 0.0 && p[[s, c]] <= 1.0,
                    "p[{s},{c}] = {} out of [0,1]",
                    p[[s, c]]
                );
            }
        }
        Ok(())
    }

    #[test]
    fn test_svc_predict_proba_binary_monotone_in_decision() -> TestResult {
        // STRUCTURAL invariant: P(classes_[1]) is monotone non-decreasing in
        // the (binary) decision_function value (higher decision -> higher
        // P(class_1)), per the sigmoid `1/(1+exp(A f + B))` contract.
        let m = proba_binary_fit()?;
        // A grid of query points sweeping from the class-0 to the class-1 side.
        let x = Array2::from_shape_vec(
            (5, 2),
            vec![1.0, 1.0, 2.5, 2.5, 3.5, 3.5, 4.5, 4.5, 6.0, 6.0],
        )
        .map_err(|_| err("shape"))?;
        let p = m.predict_proba(&x)?;
        let df = m.decision_function(&x)?;
        let bin = df.as_binary().ok_or_else(|| err("binary"))?;

        // Sort sample indices by decision value, then P(class_1) must be
        // non-decreasing along that order.
        let mut order: Vec<usize> = (0..5).collect();
        order.sort_by(|&a, &b| {
            bin[a]
                .partial_cmp(&bin[b])
                .unwrap_or(std::cmp::Ordering::Equal)
        });
        let mut prev = f64::NEG_INFINITY;
        for &s in &order {
            let p1 = p[[s, 1]];
            assert!(
                p1 >= prev - 1e-9,
                "P(class_1) not monotone in decision: sample {s} df={} p1={p1} prev={prev}",
                bin[s]
            );
            prev = p1;
        }
        Ok(())
    }

    #[test]
    fn test_svc_predict_log_proba_equals_log_of_proba() -> TestResult {
        let m = proba_binary_fit()?;
        let x = Array2::from_shape_vec((3, 2), vec![1.0, 1.0, 3.5, 3.5, 6.0, 6.0])
            .map_err(|_| err("shape"))?;
        let p = m.predict_proba(&x)?;
        let lp = m.predict_log_proba(&x)?;
        assert_eq!(lp.dim(), p.dim());
        for s in 0..3 {
            for c in 0..2 {
                assert_relative_eq!(lp[[s, c]], p[[s, c]].ln(), epsilon = 1e-12);
            }
        }
        Ok(())
    }

    fn proba_multiclass_fit() -> Result<FittedSVC<f64, LinearKernel>, FerroError> {
        let x = Array2::from_shape_vec(
            (9, 2),
            vec![
                0.0, 0.0, 0.5, 0.0, 0.0, 0.5, 5.0, 0.0, 5.5, 0.0, 5.0, 0.5, 0.0, 5.0, 0.5, 5.0,
                0.0, 5.5,
            ],
        )
        .map_err(|_| err("shape"))?;
        let y = array![0usize, 0, 0, 1, 1, 1, 2, 2, 2];
        SVC::new(LinearKernel)
            .with_c(1.0)
            .with_tol(1e-6)
            .with_max_iter(200_000)
            .with_probability(true)
            .fit(&x, &y)
    }

    #[test]
    fn test_svc_predict_proba_multiclass_rows_sum_to_one() -> TestResult {
        // 3-class: predict_proba is (n, 3), each row a probability
        // distribution (Wu-Lin-Weng coupling, `svm.cpp:2941`).
        let m = proba_multiclass_fit()?;
        let x = Array2::from_shape_vec((3, 2), vec![0.25, 0.25, 5.0, 0.25, 0.25, 5.0])
            .map_err(|_| err("shape"))?;
        let p = m.predict_proba(&x)?;
        assert_eq!(p.dim(), (3, 3));
        for s in 0..3 {
            let row_sum: f64 = (0..3).map(|c| p[[s, c]]).sum();
            assert!((row_sum - 1.0).abs() < 1e-9, "row {s} sums to {row_sum}");
            for c in 0..3 {
                assert!(
                    p[[s, c]] >= 0.0 && p[[s, c]] <= 1.0,
                    "p[{s},{c}] = {} out of [0,1]",
                    p[[s, c]]
                );
            }
        }
        Ok(())
    }

    #[test]
    fn test_sigmoid_predict_overflow_safe() {
        // sigmoid_predict matches `1/(1+exp(A f + B))` and is overflow-safe at
        // extreme decision values (`svm.cpp:2032-2040`).
        let a = -1.0f64;
        let b = 0.0;
        // f large positive -> fApB = -f large negative -> p -> 1.
        let p_pos = sigmoid_predict(1000.0, a, b);
        assert!(p_pos.is_finite() && (p_pos - 1.0).abs() < 1e-6);
        // f large negative -> p -> 0.
        let p_neg = sigmoid_predict(-1000.0, a, b);
        assert!(p_neg.is_finite() && p_neg.abs() < 1e-6);
        // f = 0 -> 1/(1+exp(0)) = 0.5.
        assert_relative_eq!(sigmoid_predict(0.0, a, b), 0.5, epsilon = 1e-12);
    }

    #[test]
    fn test_multiclass_probability_binary_reduces_to_pairwise() {
        // For k=2 the Wu-Lin-Weng coupling reduces to [r01, 1-r01].
        let mut r = Array2::<f64>::zeros((2, 2));
        r[[0, 1]] = 0.7;
        r[[1, 0]] = 0.3;
        let p = multiclass_probability(2, &r);
        assert_relative_eq!(p[0], 0.7, epsilon = 1e-6);
        assert_relative_eq!(p[1], 0.3, epsilon = 1e-6);
    }

    /// `compute_class_weight` matches `sklearn.utils.compute_class_weight`
    /// (`_classes.py:122-124` balanced formula) on the imbalanced set.
    #[test]
    fn test_compute_class_weight_balanced() {
        // 8 samples, 2 classes; class0 count=5, class1 count=3.
        // balanced[c] = 8 / (2 * count_c): [8/10, 8/6] = [0.8, 1.3333].
        let classes = [0usize, 1];
        let y = [0usize, 0, 0, 0, 0, 1, 1, 1];
        let w = compute_class_weight::<f64>(&ClassWeight::Balanced, &classes, &y);
        assert_relative_eq!(w[0], 0.8, epsilon = 1e-9);
        assert_relative_eq!(w[1], 8.0 / 6.0, epsilon = 1e-9);
        // None -> all 1.0.
        let wn = compute_class_weight::<f64>(&ClassWeight::None, &classes, &y);
        assert_eq!(wn, vec![1.0, 1.0]);
        // Explicit map, unlisted defaults to 1.0.
        let we = compute_class_weight::<f64>(&ClassWeight::Explicit(vec![(1, 5.0)]), &classes, &y);
        assert_eq!(we, vec![1.0, 5.0]);
    }
}