ferrolearn_linear/

sgd.rs

//! Stochastic Gradient Descent (SGD) linear models.
//!
//! This module provides [`SGDClassifier`] and [`SGDRegressor`], two linear
//! models trained using stochastic gradient descent. Both support online /
//! streaming learning via the [`PartialFit`] trait and a range of configurable
//! loss functions and learning-rate schedules.
//!
//! # Classifier
//!
//! ```
//! use ferrolearn_linear::sgd::{SGDClassifier, ClassifierLoss};
//! use ferrolearn_core::{Fit, Predict};
//! use ndarray::{array, Array2};
//!
//! let x = Array2::from_shape_vec((6, 2), vec![
//!     1.0, 2.0, 2.0, 3.0, 3.0, 1.0,
//!     8.0, 7.0, 9.0, 8.0, 7.0, 9.0,
//! ]).unwrap();
//! let y = array![0, 0, 0, 1, 1, 1];
//!
//! let model = SGDClassifier::<f64>::new();
//! let fitted = model.fit(&x, &y).unwrap();
//! let preds = fitted.predict(&x).unwrap();
//! assert_eq!(preds.len(), 6);
//! ```
//!
//! # Regressor
//!
//! ```
//! use ferrolearn_linear::sgd::{SGDRegressor, RegressorLoss};
//! use ferrolearn_core::{Fit, Predict};
//! use ndarray::{array, Array2};
//!
//! 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 = SGDRegressor::<f64>::new();
//! let fitted = model.fit(&x, &y).unwrap();
//! let preds = fitted.predict(&x).unwrap();
//! assert_eq!(preds.len(), 4);
//! ```
//!
//! ## REQ status (per `.design/linear/sgd.md`, mirrors `sklearn/linear_model/_stochastic_gradient.py` + `_sgd_fast.pyx.tp` @ 1.5.2, commit 156ef14)
//!
//! Parity is framed on the deterministic schedule/loss/penalty math + defaults;
//! random-shuffle full-fit weight parity is out of scope (cross-PRNG boundary,
//! `_sgd_fast.pyx.tp:579-580` vs `StdRng`). Two states only per R-DEFER-2.
//!
//! | REQ | Status | Evidence |
//! |---|---|---|
//! | REQ-1 (classifier losses hinge/log/modified_huber/squared_error incl. Hinge boundary) | SHIPPED | `impl Loss for Hinge/LogLoss/ModifiedHuber/SquaredError`. Hinge `gradient` now uses the NON-strict boundary `margin <= 1` matching `_sgd_fast.pyx.tp:224` (`if z <= threshold: return -y`). Consumer: `fn dispatch_train_binary` -> `Fit for SGDClassifier` -> `impl PipelineEstimator for SGDClassifier`. Tests: `test_hinge_loss_*`, divergence `sgd_hinge_gradient_boundary`. Closed #539. |
//! | REQ-2 (squared_hinge, perceptron) | SHIPPED | `pub struct SquaredHinge` (`loss = (1-py)^2 if >0 else 0`, `gradient = -2y(1-py)`, `_sgd_fast.pyx.tp:248-258` with `threshold=1.0`, `_stochastic_gradient.py:511`) + `pub struct Perceptron` (`Hinge(threshold=0.0)`: `loss = max(0,-py)`, `gradient = -y if py<=0 else 0`, `_sgd_fast.pyx.tp:216-226`, `_stochastic_gradient.py:512`); `enum ClassifierLoss::{SquaredHinge,Perceptron}` wired in `fn dispatch_train_binary`. Consumer: `fn dispatch_train_binary` -> `fn fit_ova` -> `Fit for SGDClassifier` -> `impl PipelineEstimator for SGDClassifier`. Tests: divergence `sgd_squared_hinge_loss` (live oracle coef `[0.0569485774276016, 0.09335170687740356]` intercept `-0.20237316143907`), `sgd_perceptron_loss` (live oracle coef `[0.009957048471181063, 0.009961042575429069]` intercept `-0.04`). Closes #523. |
//! | REQ-3 (regressor losses incl. squared_epsilon_insensitive) | SHIPPED | `pub struct SquaredEpsilonInsensitive<F> { epsilon }` (`loss = max(0,|y-p|-eps)^2`, `gradient = -2(z-eps) if z>eps; 2(-z-eps) if z<-eps; else 0` for `z=y-p`, `_sgd_fast.pyx.tp:375-387`, `_stochastic_gradient.py:1405` default `epsilon=0.1`) + `enum RegressorLoss::SquaredEpsilonInsensitive(F)` wired in `fn dispatch_train_regressor`. Consumer: `fn dispatch_train_regressor` -> `Fit for SGDRegressor` -> `impl PipelineEstimator for SGDRegressor`. Test: divergence `sgd_squared_epsilon_insensitive_loss` (live oracle single-sample coef `[0.9558857922397863, -0.47794289611989316]` intercept `0.478752180393125`, multi-sample shuffle=false coef `[0.5631419328099845, 0.41545070758814734]` intercept `0.16944283314514064`). Closes #524. |
//! | REQ-4 (L2 penalty = clamped wscale shrink) | SHIPPED | `fn train_binary_sgd`/`train_regressor_sgd` apply `shrink = max(0, 1 - eta*alpha)` then `w = w*shrink - eta*grad*x`, mirroring `w.scale(max(0, 1-eta*alpha))` (`_sgd_fast.pyx.tp:632-635`); intercept unregularized. Consumer: `Fit for SGDRegressor`/`SGDClassifier` -> `PipelineEstimator`. Test: divergence `sgd_l2_wscale_clamp`. Closed #525. |
//! | REQ-5 (l1/elasticnet + l1_ratio) | SHIPPED | `enum Penalty {L2,L1,ElasticNet}` + `pub penalty`/`pub l1_ratio` fields on `SGDClassifier`/`SGDRegressor` with `fn with_penalty`/`fn with_l1_ratio` builders (defaults `L2`/`0.15`, `_stochastic_gradient.py:1231-1256`). `fn train_binary_sgd`/`train_regressor_sgd` derive `eff` via `fn effective_l1_ratio` (`L2->0`, `L1->1`, `ElasticNet->l1_ratio`, `_sgd_fast.pyx.tp:558-561`), apply the L2 shrink `max(0, 1-(1-eff)*eta*alpha)` BEFORE the gradient add (`:632-635`), then the Tsuruoka cumulative-penalty L1 truncation with fit-persistent scalar `u` and per-feature `q` AFTER (`:656-658,750-778`, `wscale=1`). Consumer: `Fit for SGDRegressor`/`SGDClassifier` -> `PipelineEstimator`. Tests: divergence `sgd_l1_truncated_gradient` (live oracle coef [0.9204,-0.4452]), `sgd_elasticnet_l1_ratio` (l1_ratio=0.3, coef [0.92340705,-0.45723495]). Closed #526. NOTE (partial_fit+l1): `u`/`q` are scoped per `train_*_sgd` call, so they persist across the epochs of a single `fit` (the parity-critical path) but reset per `partial_fit` call. This MATCHES sklearn, which re-allocates `q=np.zeros(...)`/`u=0.0` at the top of every `_plain_sgd` call (`_sgd_fast.pyx.tp:551-556`) and only carries `t_` across `partial_fit` (`_stochastic_gradient.py` re-invokes `_plain_sgd` per call). The full `fit` path is exact. |
//! | REQ-6 (constant + invscaling schedules) | SHIPPED | `fn compute_lr`: `Constant => eta0`, `InvScaling => eta0 / t^power_t` (`_sgd_fast.pyx.tp:479,593-594`). Consumer: per-step in `fn train_binary_sgd`/`train_regressor_sgd`. Tests: `test_constant_lr`, `test_invscaling_lr`. |
//! | REQ-7 (optimal schedule t0 offset) | SHIPPED | `fn compute_lr` Optimal arm now `1/(alpha*(optimal_init + t - 1))` with `optimal_init` from `fn optimal_init` (`typw=sqrt(1/sqrt(alpha))`, `e0=typw/max(1,|gradient(1,-typw)|)`, `optimal_init=1/(e0*alpha)`), mirroring `_sgd_fast.pyx.tp:565-570,592`. Computed once per fit before the epoch loop. Consumer: `fn train_*_sgd`. Tests: `test_optimal_lr`, `test_optimal_init_matches_oracle`, divergence `sgd_optimal_schedule_t0_offset`. Closed #527. |
//! | REQ-8 (adaptive /5 + n_iter_no_change trigger) | SHIPPED | `fn convergence_tail` (shared by `fn train_binary_sgd`/`train_regressor_sgd`) divides `current_eta` by 5 (not 2) when `no_improve_count >= n_iter_no_change` AND the schedule is `Adaptive` AND `eta > 1e-6`, resetting the count — the SAME `best_loss`/`sumloss > best_loss - tol*n` machinery as convergence (`_sgd_fast.pyx.tp:697-707`). The old `>= prev_loss` 5-epoch `/2` trigger is deleted. Consumer: `Fit for SGDRegressor`/`SGDClassifier` -> `PipelineEstimator`. Test: divergence `sgd_adaptive_schedule_divisor` (live oracle coef `[0.8065190275590332, 0.15336844797680402]` intercept `0.12731338963662575`, n_iter_ 80). Closes #528. |
//! | REQ-9 (default params per estimator) | SHIPPED (classifier defaults) | `SGDClassifier::new` now sets `learning_rate=Optimal, eta0=0.0, power_t=0.5` (`_stochastic_gradient.py:1242-1244`); `fn schedule_requires_eta0` gates the `eta0>0` validation to constant/invscaling/adaptive (`_stochastic_gradient.py:149-153`). Consumer: `Fit for SGDClassifier`. Tests: divergence `sgd_classifier_default_learning_rate`, `test_sgd_classifier_default`, `test_sgd_classifier_optimal_eta0_zero_ok`. Closed #529. `epsilon` is validated to `[0, inf)` in `fn validate_reg_params` (`RegressorLoss::{Huber,EpsilonInsensitive,SquaredEpsilonInsensitive}(e)` reject `e < 0` with `FerroError::InvalidParameter { name: "epsilon" }`, mirroring `_stochastic_gradient.py:2024` `"epsilon": [Interval(Real, 0, None, closed="left")]`; test `sgd_epsilon_negative_rejected`, closes #544). Remaining missing fields (`fit_intercept`, `early_stopping`, `validation_fraction`, `average`, `warm_start`, `class_weight`, `C`) tracked under their own blockers. (`penalty`/`l1_ratio` shipped under REQ-5; `shuffle` under REQ-12; `n_iter_no_change` folded into REQ-10.) |
//! | REQ-10 (convergence best_loss/n_iter_no_change/sumloss) | SHIPPED | `fn convergence_tail` (shared by `fn train_binary_sgd`/`train_regressor_sgd`) tracks `best_loss` (running min, init `+inf`) and increments `no_improve_count` when `tol_active && sumloss > best_loss - tol*n_samples`, resetting otherwise; breaks once `no_improve_count >= hyper.n_iter_no_change` (non-adaptive), exactly mirroring `_sgd_fast.pyx.tp:688-707`. `sumloss` is now the SUM of per-sample losses over the epoch (the `/= n_samples` mean division is removed, `_sgd_fast.pyx.tp:597`); `tol_active = hyper.tol > -inf` encodes sklearn's `tol=None -> -INFINITY` disable (`:690`). The per-sample gradient is clipped to `[-1e12, 1e12]` via `fn max_dloss` before the update (`_sgd_fast.pyx.tp:546,613-620`). `n_iter_no_change` is now a settable `pub` field (default 5) with `fn with_n_iter_no_change` on both estimators, threaded through `SGDHyper`/`clf_hyper`/`reg_hyper` (`_stochastic_gradient.py` `n_iter_no_change=5`). Consumer: `Fit for SGDRegressor`/`SGDClassifier` -> `PipelineEstimator`. Test: divergence `sgd_convergence_n_iter_no_change` (live oracle coef `[0.8037686404055491, 0.16059017315681692]` intercept `0.12903834217696583`, n_iter_ 49). Closes #530. |
//! | REQ-11 (fit_intercept) | SHIPPED | `pub fit_intercept: bool` field on `SGDClassifier`/`SGDRegressor` + `fn with_fit_intercept` builders (default `true`, sklearn `_stochastic_gradient.py` `fit_intercept=True`, constraint `["boolean"]` at `:86`), threaded through `SGDHyper.fit_intercept` + `fn clf_hyper`/`reg_hyper`. `fn train_binary_sgd`/`train_regressor_sgd` gate the intercept update: `if hyper.fit_intercept { *intercept = *intercept - eta * grad; }`, mirroring `if fit_intercept == 1: intercept_update = update; ... intercept += intercept_update * intercept_decay` (`_sgd_fast.pyx.tp:639-644`, `intercept_decay=1` on the standard path). When `false` the intercept is never modified and stays at its init value `0` (`b = F::zero()` before training in `fn fit_ova`/regressor `Fit::fit`), so `coef_` matches sklearn and `intercept_` is exactly `0`. Consumer: `Fit for SGDRegressor`/`SGDClassifier` -> `PipelineEstimator`. Test: divergence `sgd_fit_intercept_false` (live oracle coef `[0.5326796739094939, 0.44573604649819804]`, intercept exactly `0.0`). Closes #531. |
//! | REQ-12 (shuffle flag) | SHIPPED | `pub shuffle: bool` field on `SGDClassifier`/`SGDRegressor` + `fn with_shuffle` builders (default `true`, `_stochastic_gradient.py:107` `shuffle=True`, constraint `["boolean"]` at `:89`), threaded through `SGDHyper.shuffle` + `fn clf_hyper`/`reg_hyper`. `fn train_binary_sgd`/`train_regressor_sgd` gate the per-epoch shuffle: `if hyper.shuffle { indices.shuffle(&mut rng); }`, mirroring `if shuffle: dataset.shuffle(seed)` (`_sgd_fast.pyx.tp:579-580`); when off, `indices` stays `0..n-1` each epoch matching sklearn's no-shuffle index order (`:581` `for i in range(n_samples)`). Consumer: `Fit for SGDRegressor`/`SGDClassifier` -> `PipelineEstimator`. Tests: divergence `sgd_shuffle_false_multisample_kernel_parity` (4-sample/2-feature/5-epoch L2 oracle coef `[0.5103165909636498, 0.42319810364130317]` intercept `0.16255331549195393`; elasticnet l1_ratio=0.3 oracle coef `[0.5102136050112174, 0.4230749783888256]` intercept `0.16265294456399926`). Closes #532. This `shuffle=false` parity ALSO validates REQ-4/REQ-5/REQ-6 (L2 shrink + elasticnet truncated gradient + constant schedule) over MULTIPLE samples and epochs against the live oracle — previously only single-sample. |
//! | REQ-13 (early_stopping + validation_fraction + n_iter_no_change-on-val-score) | SHIPPED (for the verifiable logic; the validation-split SELECTION is numpy-RNG-coupled so full fitted-coef parity is NOT oracle-verifiable — same barrier as `shuffle`) | `pub early_stopping: bool` (default `false`) + `pub validation_fraction: F` (default `0.1` via `cst`) fields on `SGDClassifier`/`SGDRegressor` + `fn with_early_stopping`/`with_validation_fraction` builders (`_stochastic_gradient.py:114-115`, constraints `["boolean"]`/`Interval(Real, 0, 1, closed="neither")` at `:524-525`), threaded through `SGDHyper.early_stopping`/`validation_fraction` (`fn clf_hyper`/`reg_hyper`; one-class hardwires `false`/`0.1`). `fn validate_validation_fraction` rejects `validation_fraction` outside the OPEN interval `(0, 1)` (`FerroError::InvalidParameter`, called from `fn validate_clf_params`/`validate_reg_params`). When `early_stopping`, the `Fit` path (`fn fit_ova` for the classifier, `SGDRegressor::fit_with_sample_weight` for the regressor) splits the data via `fn make_validation_split` BEFORE the kernel: a seeded (`StdRng::seed_from_u64(random_state.unwrap_or(0))`) hold-out of `fn validation_count` (`ceil(validation_fraction*n)`, clamped to `[1,n-1]`) samples — STRATIFIED per class for the classifier (mirrors `StratifiedShuffleSplit`, `:280-287`; the multiclass split is computed ONCE on the full `y` and SHARED across OvA subproblems, `:796`), plain `ShuffleSplit` for the regressor — returning `Err` on an empty train/val subset (`:295-307`). The kernel trains on the TRAIN subset only (sklearn does NOT refit on full data) and, at epoch-end, when early stopping, scores the CURRENT (weights,intercept) on the held-out val set via `fn convergence_tail_score` with `best_score` init `-inf`: `if tol_active && score < best_score + tol { no_improve++ } else { 0 }; if score > best_score { best_score = score }` then the SHARED `no_improve >= n_iter_no_change` adaptive-÷5/break tail (`_sgd_fast.pyx.tp:678-707`). The val score is `fn r2_score` for the regressor (`1 - SS_res/SS_tot`, `SS_tot==0` -> `1.0`/`0.0` edge) and `fn binary_accuracy` for each classifier subproblem (relabeled `{-1,+1}` target, `decision>=0 -> +1` tie convention), mirroring `_ValidationScoreCallback.__call__` = `est.score(X_val,y_val)` (R²/accuracy, `_stochastic_gradient.py:79`, `fit_binary`'s `classes=[-1,1]`+`y_i` callback at `:451-454`). `early_stopping=false` leaves the training-loss `fn convergence_tail` path byte-identical (the 25 prior divergence tests stay green). Consumer: `Fit for SGDRegressor`/`SGDClassifier` -> `impl PipelineEstimator`. VERIFIED DETERMINISTICALLY: `fn r2_score`/`binary_accuracy` against the live `sklearn.metrics.r2_score`/`accuracy_score` oracle in `mod tests` (`test_validation_r2_matches_sklearn` -> `0.8887362637362637`, `test_validation_r2_constant_y_edge_cases` -> `1.0`/`0.0`, `test_validation_binary_accuracy_matches_sklearn` -> `0.75`); the `(0,1)` constraint (divergence `sgd_validation_fraction_invalid`); the behavioral early-stop (divergence `sgd_early_stopping_stops_early` — `early_stopping=true` yields a finite model DIFFERENT from `early_stopping=false`; `sgd_early_stopping_classifier_valid`). NOT VERIFIABLE (honest): the validation-subset SELECTION uses numpy Mersenne-Twister (`ShuffleSplit`/`StratifiedShuffleSplit`) whereas ferrolearn uses `StdRng`, so the exact held-out indices — and hence the full fitted `coef_` under early stopping — are NOT cross-impl reproducible (`_stochastic_gradient.py:284-287`, the SAME PRNG barrier as `shuffle`, `_sgd_fast.pyx.tp:579-580`). `_sgd_fast.pyx.tp:678-689` / `_stochastic_gradient.py:63-79,257-310,524-525`. Closes #533. NOTE: early stopping on `partial_fit` is OFF (sklearn raises `early_stopping should be False with partial_fit`, `:147-148`); the kernel receives `val_set=None` there. |
//! | REQ-14 (average / ASGD) | SHIPPED | `pub average: usize` field on `SGDClassifier`/`SGDRegressor` (default `0` = OFF) + `fn with_average` builders (sklearn `average=True`≡`1`, `average=N`≡`N`, `average=False`≡`0`; `_stochastic_gradient.py:1256,2068`), threaded through `SGDHyper.average` + `fn clf_hyper`/`reg_hyper` (one-class path hardwires `0`). `fn train_binary_sgd`/`train_regressor_sgd` allocate `average_coef`/`average_intercept` before the epoch loop and, AFTER the weight/intercept update + L1 truncation, when `hyper.average > 0 && t >= hyper.average`, accumulate the DIRECT running mean `avg += (current - avg) / (t - average + 1)` — the plain-array equivalent of sklearn's lazy `w.add_average(..., t - average + 1)` / `average_intercept += (intercept - average_intercept) / (t - average + 1)` (`_sgd_fast.pyx.tp:646-654`); the accumulator is passive (does NOT alter the live trajectory). FINALIZE: at fit-end, `if hyper.average > 0 && hyper.average <= t { weights = average_coef; intercept = average_intercept; }`, mirroring `if self.average > 0: if self.average <= self.t_ - 1: coef_ = average_coef` (`_stochastic_gradient.py:834-836`); `t` (= `n_iter_ * n_samples`, `initial_t=0`) equals sklearn's `self.t_ - 1` (sklearn inits `self.t_ = 1`). `average=0` skips both blocks, leaving the trajectory byte-identical (the 22 prior divergence tests stay green). Consumer: `Fit for SGDRegressor`/`SGDClassifier` -> `PipelineEstimator`. Tests: divergence `sgd_average_from_start` (`SGDRegressor(average=True)` live oracle coef `[0.42614902504529534, 0.3665230497098742]` intercept `0.14648807826338486`), `sgd_average_threshold` (`SGDRegressor(average=20)`, begins mid-run, oracle coef `[0.5042444287230554, 0.41888001003992603]` intercept `0.16902090306985734`), `sgd_average_classifier` (`SGDClassifier(average=True)` oracle coef `[0.11902998815794437, 0.060826180676538694]` intercept `-0.10666666666666665`). `_sgd_fast.pyx.tp:646-654` / `_stochastic_gradient.py:834-836`. Closes #534. NOTE: averaging on the `partial_fit` path is OFF (`average` not yet carried into the `partial_fit_ova` hyper, which sets `max_iter=1`) — the full `fit` path (the parity-critical one) is exact; partial_fit ASGD state carry-over is a follow-up. |
//! | REQ-15 (class_weight + sample_weight) | SHIPPED | `pub enum ClassWeight<F> {None,Balanced,Explicit(Vec<(usize,F)>)}` + `pub class_weight` field on `SGDClassifier` (default `ClassWeight::None`) with `fn with_class_weight`; `fn compute_class_weight` returns the expanded per-class weights (`None->1.0`; `Balanced-> n_samples/(n_classes*count_c)`; `Explicit->1.0 default, override by label`) faithful to `sklearn.utils.compute_class_weight` (`sklearn/utils/class_weight.py:63-81`, `_stochastic_gradient.py:624`). `fn fit_with_sample_weight` on `SGDClassifier` AND `SGDRegressor` validates `sample_weight.len()==n_samples` (else `ShapeMismatch`); `Fit::fit` delegates with `ones(n)` (byte-identical default path — the 17 prior divergence tests stay green). `fn fit_ova` builds the per-subproblem per-sample weight `w_i = class_weight_for_sample(i) * sample_weight[i]` with the sklearn OvA mapping (binary `pos=expanded[1]`/`neg=expanded[0]`, `_stochastic_gradient.py:765-766`; multiclass class k `pos=expanded[k]`/`neg=1.0`, `:816`) and passes `&[F]` into `fn train_binary_sgd`. The kernel scales ONLY the gradient term `g = grad * sample_w[i]` (`update *= class_weight*sample_weight`, `_sgd_fast.pyx.tp:630`): the weight data term `w[j]*shrink - eta*g*x[j]` and the (gated) intercept gradient term `-eta*g` use `g`; the L2 shrink (`:632-635`), L1 truncation (`:656-658`), one-class `-2*eta*alpha` offset (`:642`) and the unweighted `sumloss` (`:597`) are UNSCALED. `fn train_regressor_sgd` mirrors the same scaling (`class_weight=1` for regression). Consumer: `Fit for SGDClassifier`/`SGDRegressor` -> `PipelineEstimator`; `fit_with_sample_weight` consumed by `Fit::fit`. Tests: divergence `sgd_class_weight_balanced` (oracle coef `[0.4806667587635881, 0.4620316761984426]` intercept `-1.2811684177087947`), `sgd_class_weight_explicit` (`{0:1.0,1:3.0}` coef `[0.5705300651778317, 0.5660417632427646]` intercept `-1.7542279278451731`), `sgd_sample_weight` (coef `[0.25648548424261425, 0.7995046753090618]` intercept `-1.221373410658307`), `sgd_class_weight_balanced_multiclass` (class-0 coef `[-0.586000112348521, -0.369263665877338]` + argmax preds), `sgd_regressor_sample_weight` (coef `[0.9425558668838198, 1.3974216923953962]` intercept `0.7259434415390171`). `_sgd_fast.pyx.tp:599-602,630` / `_stochastic_gradient.py:624,765-766,816`. Closes #535. NOTE: `class_weight`/`sample_weight` on the `partial_fit` path are uniform `1.0` (no `class_weight`/`sample_weight` arg on `PartialFit` yet) — tracked under the partial_fit surface, not this REQ. |
//! | REQ-16 (partial_fit semantics) | SHIPPED | `fn partial_fit (PartialFit for SGDClassifier/FittedSGDClassifier/SGDRegressor/FittedSGDRegressor)` sets `max_iter=1` and carries `self.t` (`_stochastic_gradient.py:581-674`). Consumer: `PartialFit` trait (`ferrolearn-core`). Tests: `test_sgd_*_partial_fit*`. |
//! | REQ-17 (multiclass one-vs-all) | SHIPPED | `fn fit_ova` (one binary per class) + `fn predict` argmax (`_stochastic_gradient.py:788-844`). Consumer: `Fit for SGDClassifier` -> `PipelineEstimator`. Test: `test_sgd_classifier_multiclass`. |
//! | REQ-18 (SGDOneClassSVM) | SHIPPED | `pub struct SGDOneClassSVM<F>` (`nu`/`fit_intercept`/`max_iter`/`tol`/`shuffle`/`learning_rate`/`eta0`/`power_t`/`random_state`/`n_iter_no_change` + `new`/`#[must_use]` builders, defaults `_stochastic_gradient.py:2245-2281`) with `fn fit_one_class` + `impl Fit<Array2<F>, ()> for SGDOneClassSVM` (X-only fit, `y` ignored, `_stochastic_gradient.py:2554`): builds `y = ones(n)`, `alpha = nu/2` (`:2588`), `penalty = L2`, `l1_ratio = 0`, `one_class = true` (`:2262-2289,2312`), inits the SGD intercept `b = 1` (offset init 0 -> `1 - 0`, `:2238,2325`), calls the reused `fn train_binary_sgd` Hinge kernel, then stores `coef_ = w`, `offset_ = 1 - b` (`:2377`). The one-class intercept term lives in `fn train_binary_sgd`: when `hyper.one_class` the gated intercept update gains `- 2*eta*alpha` (`intercept_update = -eta*grad - 2*eta*alpha`), mirroring `_sgd_fast.pyx.tp:641-642` (`if one_class: intercept_update -= 2.*eta*alpha`); `pub one_class: bool` was added to `SGDHyper` (default `false` via `fn clf_hyper`/`reg_hyper`, leaving the clf/reg intercept update byte-identical — the existing 15 divergence tests stay green). `pub struct FittedSGDOneClassSVM<F>` exposes `coef()`/`offset()`/`decision_function()` (`X·coef_ - offset_`, `:2622`)/`score_samples()` (`+ offset_ = X·coef_`, `:2639`) and `impl Predict<Array2<F>>` returning `Array1<isize>` of `+1`/`-1` (`(decision >= 0) ? +1 : -1`, `:2655-2657`). Consumer: `pub use sgd::{SGDOneClassSVM, FittedSGDOneClassSVM}` from `ferrolearn-linear/src/lib.rs` (the grandfathered public-API boundary, matching `SGDClassifier`/`SGDRegressor`). Tests: divergence `sgd_one_class_svm_decision` (live oracle nu=0.5/eta0=0.05/constant/max_iter=10/shuffle=false: coef `[0.009883660184666337, 0.009883660184666337]`, offset `1.1102230246251565e-16`, 1e-7) and `sgd_one_class_svm_predict` (nu=0.8/eta0=0.1/max_iter=15: coef `[0.20020636453962284, 0.12292535592963398]`, offset `0.10000000000000009`, predict `[1,-1,1,-1]`). `_stochastic_gradient.py:2084-2668` / `_sgd_fast.pyx.tp:639-644`. Closes #536. |
//! | REQ-19 (anti-pattern cleanup) | SHIPPED | `fn compute_lr`'s `_Phantom` arm returns `eta0` (the `unreachable!()` macro was removed earlier), and every production `F::from(<f64 literal>).unwrap()` / `F::from(<literal>).unwrap_or_else(|| ...)` constant-construction site is now `fn cst<F: Float>(x: f64) -> F { F::from(x).unwrap_or_else(F::zero) }` (a private module-level infallible-for-f32/f64 constant helper, defined after the imports). 23 call sites replaced: LogLoss `18.0`/`-18.0`/`1e18` (`_sgd_fast.pyx.tp:267-283`), SquaredError/Huber `0.5` (`:291-295,315-331`), ModifiedHuber `4.0`/`-2.0` (`:178-194`), SquaredHinge/SquaredEpsilonInsensitive/intercept/one-class `2.0` (`:254-258,379-387,641-642,2588`), and the `SGDClassifier`/`SGDRegressor`/`SGDOneClassSVM` `::new` defaults (`0.0`/`0.0001`/`0.15`/`1e-3`/`0.5`/`0.25`/`0.01`/`0.01`/`0.5`, `_stochastic_gradient.py:1242-1256,2042-2068,2245-2281`). No numeric literal changed -> byte-identical for f32/f64; all 25 `divergence_sgd_fit` + full lib/doctest suites stay green. No production panicking constant-conversion remains outside `#[cfg(test)]` in `sgd.rs` (verified by grep). Per R-APG-1 / R-CODE-2. The runtime `F::from(<usize>)` conversions (`t`, `n_samples`, `num_iter`, `count`, `from_usize`) and the deliberately-non-zero-fallback constants (`max_dloss` `1e12`->`F::max_value`, `eta_floor` `1e-6`, `divisor` `5.0`) already used `unwrap_or_else` and were already gate-compliant. Closes #537. |
//! | REQ-20 (ferray substrate migration) | NOT-STARTED | blocker #538. Still `ndarray` + `StdRng` (R-SUBSTRATE-1). |
//! | REQ-21 (non-finite input rejected) | SHIPPED | All three SGD fit entries reject any NaN/+/-inf in their float inputs BEFORE the SGD kernel with `FerroError::InvalidParameter`, mirroring sklearn's `_validate_data(force_all_finite=True)` (`_stochastic_gradient.py:1476` clf/reg base, `:2392` one-class) + `_check_sample_weight` (`:1501`) → `ValueError("Input X contains NaN.")` / `"Input y contains NaN."` / `"... contains infinity ..."`. `SGDClassifier::fit_with_sample_weight` checks X + `sample_weight` (`y: Array1<usize>` finite by type); `SGDRegressor::fit_with_sample_weight` checks X + y (`Array1<F>`) + `sample_weight`; the SEPARATE `SGDOneClassSVM::fit_one_class` arm (X-only fit, no y/sample_weight) checks X. `Fit::fit` delegates to the `fit_with_sample_weight` entries with unit weights, so the guard covers the default path too. `.iter().any(|v| !v.is_finite())` catches NaN and Inf; finite paths byte-identical. Verified vs the live sklearn 1.5.2 oracle (R-CHAR-3): NaN/+inf/-inf in X for SGDClassifier/SGDRegressor, NaN/inf in y + sample_weight for SGDRegressor, NaN in sample_weight for SGDClassifier all raise `ValueError` (`tests/divergence_linear_nonfinite_batch4.rs::sgd_*`). Non-test consumer: the existing `Fit for SGDClassifier`/`SGDRegressor` + `pub use sgd::{...}` boundary. (#2263) |

use ferrolearn_core::error::FerroError;
use ferrolearn_core::introspection::HasCoefficients;
use ferrolearn_core::pipeline::{FittedPipelineEstimator, PipelineEstimator};
use ferrolearn_core::traits::{Fit, PartialFit, Predict};
use ndarray::{Array1, Array2, ScalarOperand};
use num_traits::{Float, FromPrimitive, ToPrimitive};
use rand::SeedableRng;
use rand::seq::SliceRandom;

/// Convert an `f64` literal constant to `F`. The conversion is infallible for
/// the supported real types (`f32`/`f64`); the `F::zero()` fallback is
/// unreachable for those and exists only to keep the call non-panicking
/// (no `.unwrap()` in production, per the anti-pattern gate R-APG-1 / R-CODE-2).
#[inline]
fn cst<F: Float>(x: f64) -> F {
    F::from(x).unwrap_or_else(F::zero)
}

// ---------------------------------------------------------------------------
// Loss functions
// ---------------------------------------------------------------------------

/// A loss function for SGD optimization.
///
/// Provides the loss value and its gradient with respect to the prediction.
pub trait Loss<F: Float>: Clone + Send + Sync {
    /// Compute the loss for a single sample.
    fn loss(&self, y_true: F, y_pred: F) -> F;

    /// Compute the gradient of the loss with respect to `y_pred`.
    fn gradient(&self, y_true: F, y_pred: F) -> F;
}

/// Hinge loss for linear SVM-style classification.
///
/// `L(y, p) = max(0, 1 - y * p)` where `y in {-1, +1}`.
#[derive(Debug, Clone, Copy)]
pub struct Hinge;

impl<F: Float> Loss<F> for Hinge {
    fn loss(&self, y_true: F, y_pred: F) -> F {
        let margin = y_true * y_pred;
        if margin < F::one() {
            F::one() - margin
        } else {
            F::zero()
        }
    }

    fn gradient(&self, y_true: F, y_pred: F) -> F {
        // sklearn `Hinge.dloss` uses a NON-strict boundary at the threshold
        // (`_sgd_fast.pyx.tp:224`: `if z <= self.threshold: return -y`), so at
        // the exact margin `z == 1` the gradient is `-y`, not `0`.
        let margin = y_true * y_pred;
        if margin <= F::one() {
            -y_true
        } else {
            F::zero()
        }
    }
}

/// Squared hinge loss for (quadratically penalized) linear SVM classification.
///
/// `L(y, p) = max(0, 1 - y * p)^2` where `y in {-1, +1}`. This is sklearn's
/// `SquaredHinge(threshold=1.0)` (`_sgd_fast.pyx.tp:232-258`); the
/// `squared_hinge` classifier loss maps to it (`_stochastic_gradient.py:511`).
#[derive(Debug, Clone, Copy)]
pub struct SquaredHinge;

impl<F: Float> Loss<F> for SquaredHinge {
    fn loss(&self, y_true: F, y_pred: F) -> F {
        // `_sgd_fast.pyx.tp:248-252`: `z = threshold - p*y; z*z if z > 0 else 0`
        // with `threshold = 1.0` (`_stochastic_gradient.py:511`).
        let z = F::one() - y_pred * y_true;
        if z > F::zero() { z * z } else { F::zero() }
    }

    fn gradient(&self, y_true: F, y_pred: F) -> F {
        // `_sgd_fast.pyx.tp:254-258`: `z = threshold - p*y; -2*y*z if z > 0 else 0`.
        let z = F::one() - y_pred * y_true;
        if z > F::zero() {
            -cst::<F>(2.0) * y_true * z
        } else {
            F::zero()
        }
    }
}

/// Perceptron loss for linear classification.
///
/// `L(y, p) = max(0, -y * p)` where `y in {-1, +1}`. This is sklearn's
/// `Hinge(threshold=0.0)` (`_sgd_fast.pyx.tp:200-226`); the `perceptron`
/// classifier loss maps to it (`_stochastic_gradient.py:512`). The existing
/// [`Hinge`] hardcodes `threshold = 1.0`, so this is a separate type.
#[derive(Debug, Clone, Copy)]
pub struct Perceptron;

impl<F: Float> Loss<F> for Perceptron {
    fn loss(&self, y_true: F, y_pred: F) -> F {
        // `_sgd_fast.pyx.tp:216-220`: `z = p*y; threshold - z if z <= threshold
        // else 0` with `threshold = 0.0` (`_stochastic_gradient.py:512`), i.e.
        // `max(0, -z)`.
        let z = y_pred * y_true;
        if z <= F::zero() { -z } else { F::zero() }
    }

    fn gradient(&self, y_true: F, y_pred: F) -> F {
        // `_sgd_fast.pyx.tp:222-226`: `z = p*y; -y if z <= threshold else 0`
        // with `threshold = 0.0`.
        let z = y_pred * y_true;
        if z <= F::zero() { -y_true } else { F::zero() }
    }
}

/// Log loss (logistic regression / cross-entropy).
///
/// `L(y, p) = log(1 + exp(-y * p))` where `y in {-1, +1}`.
#[derive(Debug, Clone, Copy)]
pub struct LogLoss;

impl<F: Float> Loss<F> for LogLoss {
    fn loss(&self, y_true: F, y_pred: F) -> F {
        let z = y_true * y_pred;
        if z > cst(18.0) {
            (-z).exp()
        } else if z < cst(-18.0) {
            -z
        } else {
            (F::one() + (-z).exp()).ln()
        }
    }

    fn gradient(&self, y_true: F, y_pred: F) -> F {
        let z = y_true * y_pred;
        let exp_nz = if z > cst(18.0) {
            (-z).exp()
        } else if z < cst(-18.0) {
            cst(1e18)
        } else {
            (-z).exp()
        };
        -y_true * exp_nz / (F::one() + exp_nz)
    }
}

/// Squared error loss for regression.
///
/// `L(y, p) = 0.5 * (y - p)^2`.
#[derive(Debug, Clone, Copy)]
pub struct SquaredError;

impl<F: Float> Loss<F> for SquaredError {
    fn loss(&self, y_true: F, y_pred: F) -> F {
        let diff = y_true - y_pred;
        cst::<F>(0.5) * diff * diff
    }

    fn gradient(&self, y_true: F, y_pred: F) -> F {
        y_pred - y_true
    }
}

/// Modified Huber loss for classification.
///
/// Smooth approximation to hinge with quadratic behaviour near the margin:
///
/// ```text
/// L(y, p) = max(0, 1 - y*p)^2   if y*p >= -1
///         = -4 * y * p            otherwise
/// ```
#[derive(Debug, Clone, Copy)]
pub struct ModifiedHuber;

impl<F: Float> Loss<F> for ModifiedHuber {
    fn loss(&self, y_true: F, y_pred: F) -> F {
        let z = y_true * y_pred;
        if z >= -F::one() {
            let margin = F::one() - z;
            if margin > F::zero() {
                margin * margin
            } else {
                F::zero()
            }
        } else {
            -cst::<F>(4.0) * z
        }
    }

    fn gradient(&self, y_true: F, y_pred: F) -> F {
        let z = y_true * y_pred;
        if z >= -F::one() {
            if z < F::one() {
                cst::<F>(-2.0) * y_true * (F::one() - z)
            } else {
                F::zero()
            }
        } else {
            -cst::<F>(4.0) * y_true
        }
    }
}

/// Huber loss for robust regression.
///
/// `L(y, p) = 0.5 * (y - p)^2` if `|y - p| <= epsilon`, else
/// `epsilon * (|y - p| - 0.5 * epsilon)`.
#[derive(Debug, Clone, Copy)]
pub struct Huber<F> {
    /// Threshold parameter for switching from quadratic to linear loss.
    pub epsilon: F,
}

impl<F: Float + Send + Sync> Loss<F> for Huber<F> {
    fn loss(&self, y_true: F, y_pred: F) -> F {
        let diff = y_true - y_pred;
        let abs_diff = diff.abs();
        if abs_diff <= self.epsilon {
            cst::<F>(0.5) * diff * diff
        } else {
            self.epsilon * (abs_diff - cst::<F>(0.5) * self.epsilon)
        }
    }

    fn gradient(&self, y_true: F, y_pred: F) -> F {
        let diff = y_pred - y_true;
        let abs_diff = diff.abs();
        if abs_diff <= self.epsilon {
            diff
        } else if diff > F::zero() {
            self.epsilon
        } else {
            -self.epsilon
        }
    }
}

/// Epsilon-insensitive loss for support vector regression.
///
/// `L(y, p) = max(0, |y - p| - epsilon)`.
#[derive(Debug, Clone, Copy)]
pub struct EpsilonInsensitive<F> {
    /// Insensitivity margin.
    pub epsilon: F,
}

impl<F: Float + Send + Sync> Loss<F> for EpsilonInsensitive<F> {
    fn loss(&self, y_true: F, y_pred: F) -> F {
        let diff = (y_true - y_pred).abs();
        if diff > self.epsilon {
            diff - self.epsilon
        } else {
            F::zero()
        }
    }

    fn gradient(&self, y_true: F, y_pred: F) -> F {
        let diff = y_pred - y_true;
        if diff > self.epsilon {
            F::one()
        } else if diff < -self.epsilon {
            -F::one()
        } else {
            F::zero()
        }
    }
}

/// Squared epsilon-insensitive loss for support vector regression.
///
/// `L(y, p) = max(0, |y - p| - epsilon)^2`. This is sklearn's
/// `SquaredEpsilonInsensitive` (`_sgd_fast.pyx.tp:364-388`); the
/// `squared_epsilon_insensitive` regressor loss maps to it
/// (`_stochastic_gradient.py:1405`, default `epsilon = DEFAULT_EPSILON = 0.1`).
#[derive(Debug, Clone, Copy)]
pub struct SquaredEpsilonInsensitive<F> {
    /// Insensitivity margin.
    pub epsilon: F,
}

impl<F: Float + Send + Sync> Loss<F> for SquaredEpsilonInsensitive<F> {
    fn loss(&self, y_true: F, y_pred: F) -> F {
        // `_sgd_fast.pyx.tp:375-377`: `ret = |y - p| - epsilon;
        // ret*ret if ret > 0 else 0`.
        let ret = (y_true - y_pred).abs() - self.epsilon;
        if ret > F::zero() {
            ret * ret
        } else {
            F::zero()
        }
    }

    fn gradient(&self, y_true: F, y_pred: F) -> F {
        // `_sgd_fast.pyx.tp:379-387`: `z = y - p;
        // -2*(z-epsilon) if z > epsilon; 2*(-z-epsilon) if z < -epsilon; else 0`.
        let two = cst::<F>(2.0);
        let z = y_true - y_pred;
        if z > self.epsilon {
            -two * (z - self.epsilon)
        } else if z < -self.epsilon {
            two * (-z - self.epsilon)
        } else {
            F::zero()
        }
    }
}

// ---------------------------------------------------------------------------
// Learning rate schedules
// ---------------------------------------------------------------------------

/// Learning rate schedule for SGD.
#[derive(Debug, Clone, Copy)]
pub enum LearningRateSchedule<F> {
    /// Fixed learning rate `eta0` throughout training.
    Constant,
    /// Optimal schedule: `eta = 1 / (alpha * t)`.
    Optimal,
    /// Inverse scaling: `eta = eta0 / t^power_t`.
    InvScaling,
    /// Adaptive: starts at `eta0`, halved when loss fails to decrease for
    /// 5 consecutive epochs. Stops when `eta < 1e-6`.
    Adaptive,
    #[doc(hidden)]
    _Phantom(std::marker::PhantomData<F>),
}

/// Compute the learning rate for a given step.
///
/// `optimal_init` is the `t0` offset of the `optimal` schedule, derived once
/// per fit from `alpha` and the loss's `dloss(1, -typw)` bound (see
/// [`optimal_init`]). It is ignored by the other schedules.
fn compute_lr<F: Float>(
    schedule: &LearningRateSchedule<F>,
    eta0: F,
    alpha: F,
    power_t: F,
    optimal_init: F,
    t: usize,
) -> F {
    let t_f = F::from(t.max(1)).unwrap_or_else(F::one);
    match schedule {
        LearningRateSchedule::Constant => eta0,
        // sklearn `_sgd_fast.pyx.tp:592`: `eta = 1/(alpha*(optimal_init+t-1))`,
        // so the first sample (t=1) sees `eta = 1/(alpha*optimal_init) = e0`.
        LearningRateSchedule::Optimal => F::one() / (alpha * (optimal_init + t_f - F::one())),
        LearningRateSchedule::InvScaling => eta0 / t_f.powf(power_t),
        LearningRateSchedule::Adaptive => eta0,
        // `_Phantom` is an uninhabited marker arm; fall back to `eta0` rather
        // than aborting (R-APG-1 forbids the unreach macro in production).
        LearningRateSchedule::_Phantom(_) => eta0,
    }
}

/// Compute the `optimal` schedule's `t0` offset `optimal_init`.
///
/// Mirrors `_sgd_fast.pyx.tp:565-570`:
/// `typw = sqrt(1/sqrt(alpha))`,
/// `initial_eta0 = typw / max(1, dloss(1, -typw))`,
/// `optimal_init = 1/(initial_eta0 * alpha)`.
///
/// sklearn calls `loss.dloss(1.0, -typw)` where the cython signature is
/// `dloss(self, y, p)` (y first, p second), so `y = 1.0`, `p = -typw`. The
/// ferrolearn signature is `gradient(y_true, y_pred)`, mapping to
/// `gradient(1.0, -typw)`; its absolute value matches `max(1.0, dloss(...))`.
/// Returns `1.0` when `alpha == 0` (the schedule is unused / guarded upstream).
fn optimal_init<F, L>(loss_fn: &L, alpha: F) -> F
where
    F: Float,
    L: Loss<F>,
{
    if alpha <= F::zero() {
        return F::one();
    }
    let typw = (F::one() / alpha.sqrt()).sqrt();
    let dloss = loss_fn.gradient(F::one(), -typw).abs();
    let initial_eta0 = typw / dloss.max(F::one());
    F::one() / (initial_eta0 * alpha)
}

// ---------------------------------------------------------------------------
// Class weights
// ---------------------------------------------------------------------------

/// Per-class weighting strategy for [`SGDClassifier`].
///
/// Mirrors sklearn's `class_weight` parameter
/// (`_stochastic_gradient.py` constraint `[dict, "balanced", None]`); the
/// expanded per-class weights are computed by [`compute_class_weight`] following
/// `sklearn.utils.compute_class_weight` semantics and fed into the per-sample
/// `update *= class_weight * sample_weight` scaling (`_sgd_fast.pyx.tp:630`).
#[derive(Debug, Clone)]
pub enum ClassWeight<F> {
    /// Uniform weights (all classes weighted `1.0`). The default.
    None,
    /// Balanced weights `n_samples / (n_classes * count_c)` per class `c`,
    /// matching `sklearn.utils.compute_class_weight("balanced", ...)`
    /// (`class_weight.py:73`).
    Balanced,
    /// Explicit class-label -> weight map. Classes absent from the map default
    /// to `1.0`, matching the dict branch of `compute_class_weight`
    /// (`class_weight.py:77-81`).
    Explicit(Vec<(usize, F)>),
}

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

// ---------------------------------------------------------------------------
// Penalty (regularization term)
// ---------------------------------------------------------------------------

/// Regularization penalty for SGD.
///
/// Mirrors sklearn's `penalty` parameter
/// (`_stochastic_gradient.py:997-1012`). The effective `l1_ratio` passed to the
/// kernel is derived per `_sgd_fast.pyx.tp:558-561`: `l2 -> 0.0`, `l1 -> 1.0`,
/// `elasticnet -> user l1_ratio`.
#[derive(Debug, Clone, Copy)]
pub enum Penalty {
    /// L2 (ridge) penalty — multiplicative `wscale` shrink only (the default).
    L2,
    /// L1 (lasso) penalty — Tsuruoka cumulative-penalty truncated gradient.
    L1,
    /// Elastic-net — convex mix of L2 and L1 controlled by `l1_ratio`.
    ElasticNet,
}

/// Compute the effective `l1_ratio` for the truncated-gradient kernel.
///
/// Mirrors `_sgd_fast.pyx.tp:558-561`: `L2 -> 0.0`, `L1 -> 1.0`,
/// `ElasticNet -> user l1_ratio`.
fn effective_l1_ratio<F: Float>(penalty: Penalty, l1_ratio: F) -> F {
    match penalty {
        Penalty::L2 => F::zero(),
        Penalty::L1 => F::one(),
        Penalty::ElasticNet => l1_ratio,
    }
}

/// The `MAX_DLOSS` gradient clip bound (`_sgd_fast.pyx.tp:546`, `1e12`).
///
/// sklearn clips `dloss` to `[-MAX_DLOSS, MAX_DLOSS]` before forming the update
/// `update = -eta * dloss` (`_sgd_fast.pyx.tp:613-620`) to avoid numerical
/// instabilities. Falls back to `F::max_value()` (an even looser, never-active
/// bound) if `1e12` is not representable in `F` — so the clamp is always a safe
/// no-op widening rather than a panic.
#[inline]
fn max_dloss<F: Float>() -> F {
    F::from(1e12_f64).unwrap_or_else(F::max_value)
}

/// Coefficient of determination `R^2 = 1 - SS_res / SS_tot` of a linear model
/// `(weights, intercept)` on `(x_val, y_val)`.
///
/// Mirrors `sklearn.metrics.r2_score` for the dense single-output, uniformly
/// weighted case (`RegressorMixin.score` -> `r2_score`, the regressor
/// `_ValidationScoreCallback.__call__`, `_stochastic_gradient.py:79`):
/// `SS_res = sum((y - y_pred)^2)`, `SS_tot = sum((y - mean(y))^2)`,
/// `R^2 = 1 - SS_res/SS_tot`. The degenerate `SS_tot == 0` (constant `y_val`)
/// case follows sklearn's `_metrics/_regression.py` convention: a perfect
/// `SS_res == 0` scores `1.0`, otherwise `0.0`
/// (`r2_score` `nonzero_denominator`/`nonzero_numerator` branch). Returns `0.0`
/// for an empty validation set (no information).
#[must_use]
fn r2_score<F: Float>(
    weights: &Array1<F>,
    intercept: F,
    x_val: &Array2<F>,
    y_val: &Array1<F>,
) -> F {
    let n = y_val.len();
    if n == 0 {
        return F::zero();
    }
    let n_f = F::from(n).unwrap_or_else(F::one);
    let mean = y_val.iter().fold(F::zero(), |acc, &v| acc + v) / n_f;
    let mut ss_res = F::zero();
    let mut ss_tot = F::zero();
    let n_features = weights.len();
    for i in 0..n {
        let xi = x_val.row(i);
        let mut pred = intercept;
        for j in 0..n_features {
            pred = pred + weights[j] * xi[j];
        }
        let res = y_val[i] - pred;
        ss_res = ss_res + res * res;
        let dev = y_val[i] - mean;
        ss_tot = ss_tot + dev * dev;
    }
    if ss_tot > F::zero() {
        F::one() - ss_res / ss_tot
    } else if ss_res > F::zero() {
        // constant `y_val` but imperfect prediction -> R^2 = 0 (sklearn).
        F::zero()
    } else {
        // constant `y_val` and perfect prediction -> R^2 = 1 (sklearn).
        F::one()
    }
}

/// Binary classification accuracy of a linear decision on `(x_val, y_val)` where
/// `y_val` is the relabeled `{-1, +1}` target.
///
/// Mirrors the classifier `_ValidationScoreCallback.__call__`
/// (`_stochastic_gradient.py:79`, `est.score(X_val, y_val)` ->
/// `ClassifierMixin.score` -> `accuracy_score`) for one One-vs-All binary
/// subproblem: the callback is built with `classes = np.array([-1, 1])` and the
/// relabeled binary target `y_i` (`fit_binary`, `:451-454`), so the score is the
/// fraction of validation samples whose binary decision `sign(w·x + b)` matches
/// the relabeled label. The decision uses the `>= 0 -> +1` tie convention
/// matching [`FittedSGDClassifier::predict`] (`scores[i] >= 0 -> classes[1]`).
/// Returns `0.0` for an empty validation set.
#[must_use]
fn binary_accuracy<F: Float>(
    weights: &Array1<F>,
    intercept: F,
    x_val: &Array2<F>,
    y_val: &Array1<F>,
) -> F {
    let n = y_val.len();
    if n == 0 {
        return F::zero();
    }
    let n_features = weights.len();
    let mut correct = 0usize;
    for i in 0..n {
        let xi = x_val.row(i);
        let mut decision = intercept;
        for j in 0..n_features {
            decision = decision + weights[j] * xi[j];
        }
        let pred = if decision >= F::zero() {
            F::one()
        } else {
            -F::one()
        };
        // `y_val[i]` is the relabeled `{-1, +1}` target; a positive product
        // means `pred` and the label agree.
        if pred * y_val[i] > F::zero() {
            correct += 1;
        }
    }
    F::from(correct).unwrap_or_else(F::zero) / F::from(n).unwrap_or_else(F::one)
}

/// Shared SGD epoch-end convergence / adaptive-eta tail.
///
/// Mirrors `_sgd_fast.pyx.tp:688-707` exactly. `epoch_sumloss` is the SUM of
/// per-sample losses over the epoch (`:597`, NOT the printed mean `sumloss /
/// train_count`). The criterion increments `no_improve_count` whenever
/// `sumloss > best_loss - tol * train_count` (an epoch that fails to beat the
/// running minimum by at least `tol*n`) and resets it otherwise; `best_loss`
/// tracks the running minimum. Once `no_improve_count >= n_iter_no_change`,
/// under the adaptive schedule (`eta > 1e-6`) `eta` is divided by 5 and the
/// count reset (`:699-701`); otherwise the caller breaks (convergence, `:702`).
///
/// Returns `true` iff the epoch loop should `break` (convergence). The two
/// branches are mutually exclusive, exactly as upstream: adaptive decays eta
/// and keeps running, non-adaptive (or eta already `<= 1e-6`) stops.
#[allow(clippy::too_many_arguments, reason = "mirrors the upstream epoch tail")]
#[inline]
fn convergence_tail<F: Float>(
    epoch_sumloss: F,
    best_loss: &mut F,
    no_improve_count: &mut usize,
    current_eta: &mut F,
    tol_active: bool,
    tol: F,
    n_samples: usize,
    n_iter_no_change: usize,
    adaptive: bool,
) -> bool {
    // `_sgd_fast.pyx.tp:690-693`: training-loss branch (early_stopping=False).
    let n = F::from(n_samples).unwrap_or_else(F::zero);
    if tol_active && epoch_sumloss > *best_loss - tol * n {
        *no_improve_count += 1;
    } else {
        *no_improve_count = 0;
    }
    // `:694-695`: track the running minimum.
    if epoch_sumloss < *best_loss {
        *best_loss = epoch_sumloss;
    }
    // `:698-707`: convergence break OR adaptive eta/=5.
    if *no_improve_count >= n_iter_no_change {
        // `:699`: `if learning_rate == ADAPTIVE and eta > 1e-6`.
        let eta_floor = F::from(1e-6_f64).unwrap_or_else(F::zero);
        let divisor = F::from(5.0_f64).unwrap_or_else(F::one);
        if adaptive && *current_eta > eta_floor {
            *current_eta = *current_eta / divisor;
            *no_improve_count = 0;
            false
        } else {
            true
        }
    } else {
        false
    }
}

/// Score-based SGD epoch-end convergence / adaptive-eta tail (early stopping).
///
/// Mirrors the `if early_stopping:` branch of `_sgd_fast.pyx.tp:678-707`. The
/// validation `score` (R^2 for the regressor, binary accuracy for the
/// classifier subproblem) replaces the training loss, and — because a HIGHER
/// score is better — `best_score` is initialized to `-inf` and the criterion
/// flips sense relative to [`convergence_tail`]: `no_improvement_count`
/// increments when `score < best_score + tol` (the epoch fails to beat the best
/// score so far by at least `tol`, `:682-685`) and resets otherwise; `best_score`
/// tracks the running MAXIMUM (`:686-687`). The shared
/// `no_improvement_count >= n_iter_no_change` tail (`:698-707`, adaptive eta/=5
/// or break) is identical to the loss path.
///
/// Returns `true` iff the epoch loop should `break` (convergence).
#[allow(clippy::too_many_arguments, reason = "mirrors the upstream epoch tail")]
#[inline]
fn convergence_tail_score<F: Float>(
    score: F,
    best_score: &mut F,
    no_improve_count: &mut usize,
    current_eta: &mut F,
    tol_active: bool,
    tol: F,
    n_iter_no_change: usize,
    adaptive: bool,
) -> bool {
    // `_sgd_fast.pyx.tp:682-685`: validation-score branch (early_stopping=True).
    if tol_active && score < *best_score + tol {
        *no_improve_count += 1;
    } else {
        *no_improve_count = 0;
    }
    // `:686-687`: track the running maximum.
    if score > *best_score {
        *best_score = score;
    }
    // `:698-707`: convergence break OR adaptive eta/=5 (shared with the loss
    // path — byte-identical logic).
    if *no_improve_count >= n_iter_no_change {
        let eta_floor = F::from(1e-6_f64).unwrap_or_else(F::zero);
        let divisor = F::from(5.0_f64).unwrap_or_else(F::one);
        if adaptive && *current_eta > eta_floor {
            *current_eta = *current_eta / divisor;
            *no_improve_count = 0;
            false
        } else {
            true
        }
    } else {
        false
    }
}

// ---------------------------------------------------------------------------
// Classifier loss enum
// ---------------------------------------------------------------------------

/// Available loss functions for [`SGDClassifier`].
#[derive(Debug, Clone, Copy)]
pub enum ClassifierLoss {
    /// Hinge loss (linear SVM).
    Hinge,
    /// Squared hinge loss (quadratically penalized SVM,
    /// `_stochastic_gradient.py:511`).
    SquaredHinge,
    /// Perceptron loss (`Hinge(threshold=0.0)`,
    /// `_stochastic_gradient.py:512`).
    Perceptron,
    /// Log loss (logistic regression).
    Log,
    /// Squared error loss.
    SquaredError,
    /// Modified Huber loss.
    ModifiedHuber,
}

/// Available loss functions for [`SGDRegressor`].
#[derive(Debug, Clone, Copy)]
pub enum RegressorLoss<F> {
    /// Squared error loss (default).
    SquaredError,
    /// Huber loss with the given epsilon.
    Huber(F),
    /// Epsilon-insensitive loss with the given epsilon.
    EpsilonInsensitive(F),
    /// Squared epsilon-insensitive loss with the given epsilon
    /// (`_stochastic_gradient.py:1405`, `_sgd_fast.pyx.tp:364-388`).
    SquaredEpsilonInsensitive(F),
}

// ---------------------------------------------------------------------------
// SGDClassifier
// ---------------------------------------------------------------------------

/// Stochastic Gradient Descent classifier.
///
/// Supports binary classification via a decision boundary and multiclass
/// classification via one-vs-all decomposition.
///
/// # Type Parameters
///
/// - `F`: The floating-point type (`f32` or `f64`).
///
/// # Examples
///
/// ```
/// use ferrolearn_linear::sgd::SGDClassifier;
/// use ferrolearn_core::{Fit, Predict};
/// use ndarray::{array, Array2};
///
/// let x = Array2::from_shape_vec((6, 2), vec![
///     1.0, 2.0, 2.0, 3.0, 3.0, 1.0,
///     8.0, 7.0, 9.0, 8.0, 7.0, 9.0,
/// ]).unwrap();
/// let y = array![0, 0, 0, 1, 1, 1];
///
/// let clf = SGDClassifier::<f64>::new();
/// let fitted = clf.fit(&x, &y).unwrap();
/// let preds = fitted.predict(&x).unwrap();
/// ```
#[derive(Debug, Clone)]
pub struct SGDClassifier<F> {
    /// The loss function to use.
    pub loss: ClassifierLoss,
    /// The learning rate schedule.
    pub learning_rate: LearningRateSchedule<F>,
    /// Initial learning rate.
    pub eta0: F,
    /// Regularization strength (`alpha`).
    pub alpha: F,
    /// Regularization penalty (`l2`/`l1`/`elasticnet`). Defaults to `L2`.
    pub penalty: Penalty,
    /// Elastic-net mixing parameter; only used when `penalty == ElasticNet`.
    /// Defaults to `0.15` (sklearn default).
    pub l1_ratio: F,
    /// Maximum number of passes over the training data.
    pub max_iter: usize,
    /// Convergence tolerance. Training stops when the loss improvement
    /// is below this threshold.
    pub tol: F,
    /// Optional random seed for sample shuffling.
    pub random_state: Option<u64>,
    /// Power parameter for inverse scaling schedule.
    pub power_t: F,
    /// Whether to shuffle the training data after each epoch. Defaults to
    /// `true` (sklearn `SGDClassifier(shuffle=True)`,
    /// `_stochastic_gradient.py:107`).
    pub shuffle: bool,
    /// Number of consecutive epochs with no loss improvement (beyond `tol`)
    /// before convergence triggers, or — under the `adaptive` schedule — before
    /// `eta` is divided by 5. Defaults to `5` (sklearn
    /// `_stochastic_gradient.py` `n_iter_no_change=5`, `_sgd_fast.pyx.tp:698`).
    pub n_iter_no_change: usize,
    /// Whether to fit (update) the intercept. Defaults to `true` (sklearn
    /// `SGDClassifier(fit_intercept=True)`, `_stochastic_gradient.py:104`,
    /// constraint `["boolean"]` at `:86`). When `false` the intercept is never
    /// updated and stays at its init value `0` (`_sgd_fast.pyx.tp:639-644`: the
    /// intercept update is gated on `if fit_intercept == 1`).
    pub fit_intercept: bool,
    /// Per-class weighting strategy. Defaults to [`ClassWeight::None`] (uniform).
    /// Mirrors sklearn's `class_weight` parameter (default `None`); the expanded
    /// weights scale the per-sample gradient term via
    /// `update *= class_weight * sample_weight` (`_sgd_fast.pyx.tp:599-602,630`).
    pub class_weight: ClassWeight<F>,
    /// Averaged-SGD (ASGD) threshold. `0` disables averaging (the default,
    /// matching sklearn `average=False`); `1` averages from the first step
    /// (sklearn `average=True`); `N > 1` begins averaging once the global step
    /// counter `t >= N` (sklearn `average=N`). The averaged weights/intercept
    /// replace the plain ones at fit-end when `average <= self.t_ - 1`
    /// (`_sgd_fast.pyx.tp:646-654`, `_stochastic_gradient.py:834-836`).
    pub average: usize,
    /// Whether to stop training early based on a held-out validation score.
    /// Defaults to `false` (sklearn `SGDClassifier(early_stopping=False)`,
    /// `_stochastic_gradient.py:114`, constraint `["boolean"]` at `:524`). When
    /// `true`, [`validation_fraction`](Self::validation_fraction) of the training
    /// data is held out (stratified) as a validation set and the epoch-end
    /// convergence rule uses the validation accuracy of each One-vs-All binary
    /// subproblem instead of the training loss (`_sgd_fast.pyx.tp:678-687`).
    pub early_stopping: bool,
    /// Fraction of the training data held out as the validation set when
    /// [`early_stopping`](Self::early_stopping) is `true`. Defaults to `0.1`
    /// (sklearn `validation_fraction=0.1`, `_stochastic_gradient.py:115`). Must
    /// lie in the open interval `(0, 1)`
    /// (constraint `Interval(Real, 0, 1, closed="neither")` at `:525`).
    pub validation_fraction: F,
}

impl<F: Float> SGDClassifier<F> {
    /// Create a new `SGDClassifier` with default settings.
    ///
    /// Defaults match scikit-learn's `SGDClassifier.__init__`
    /// (`_stochastic_gradient.py:1231-1256`): `loss = Hinge`,
    /// `learning_rate = Optimal`, `eta0 = 0.0`, `alpha = 0.0001`,
    /// `penalty = L2`, `l1_ratio = 0.15`, `max_iter = 1000`, `tol = 1e-3`,
    /// `power_t = 0.5`.
    #[must_use]
    pub fn new() -> Self {
        Self {
            loss: ClassifierLoss::Hinge,
            learning_rate: LearningRateSchedule::Optimal,
            eta0: cst(0.0),
            alpha: cst(0.0001),
            penalty: Penalty::L2,
            l1_ratio: cst(0.15),
            max_iter: 1000,
            tol: cst(1e-3),
            random_state: None,
            power_t: cst(0.5),
            shuffle: true,
            n_iter_no_change: 5,
            fit_intercept: true,
            class_weight: ClassWeight::None,
            average: 0,
            early_stopping: false,
            validation_fraction: cst(0.1),
        }
    }

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

    /// Set the learning rate schedule.
    #[must_use]
    pub fn with_learning_rate(mut self, lr: LearningRateSchedule<F>) -> Self {
        self.learning_rate = lr;
        self
    }

    /// Set the initial learning rate.
    #[must_use]
    pub fn with_eta0(mut self, eta0: F) -> Self {
        self.eta0 = eta0;
        self
    }

    /// Set the regularization strength.
    #[must_use]
    pub fn with_alpha(mut self, alpha: F) -> Self {
        self.alpha = alpha;
        self
    }

    /// Set the regularization penalty (`l2`/`l1`/`elasticnet`).
    #[must_use]
    pub fn with_penalty(mut self, penalty: Penalty) -> Self {
        self.penalty = penalty;
        self
    }

    /// Set the elastic-net mixing parameter (`l1_ratio`, used only when
    /// `penalty == ElasticNet`).
    #[must_use]
    pub fn with_l1_ratio(mut self, l1_ratio: F) -> Self {
        self.l1_ratio = l1_ratio;
        self
    }

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

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

    /// Set the random seed for reproducibility.
    #[must_use]
    pub fn with_random_state(mut self, seed: u64) -> Self {
        self.random_state = Some(seed);
        self
    }

    /// Set the power parameter for inverse scaling.
    #[must_use]
    pub fn with_power_t(mut self, power_t: F) -> Self {
        self.power_t = power_t;
        self
    }

    /// Set whether the training data is shuffled after each epoch.
    ///
    /// Mirrors sklearn's `shuffle` parameter (default `True`,
    /// `_stochastic_gradient.py:107`). With `false` the samples are visited in
    /// index order `0..n-1` every epoch (`_sgd_fast.pyx.tp:579-581`), making the
    /// fit fully deterministic and cross-impl comparable to sklearn.
    #[must_use]
    pub fn with_shuffle(mut self, shuffle: bool) -> Self {
        self.shuffle = shuffle;
        self
    }

    /// Set the number of consecutive non-improving epochs before convergence
    /// (or, under the `adaptive` schedule, before `eta` is divided by 5).
    ///
    /// Mirrors sklearn's `n_iter_no_change` parameter (default `5`,
    /// `_stochastic_gradient.py`); the epoch-end stop rule at
    /// `_sgd_fast.pyx.tp:698` triggers once `no_improvement_count` reaches it.
    #[must_use]
    pub fn with_n_iter_no_change(mut self, n_iter_no_change: usize) -> Self {
        self.n_iter_no_change = n_iter_no_change;
        self
    }

    /// Set whether the intercept (bias) term is fit.
    ///
    /// Mirrors sklearn's `fit_intercept` parameter (default `True`,
    /// `_stochastic_gradient.py:104`, constraint `["boolean"]` at `:86`). With
    /// `false` the intercept update is skipped every step
    /// (`_sgd_fast.pyx.tp:639-644`: `if fit_intercept == 1: ... intercept += ...`)
    /// and the intercept stays at its init value `0`.
    #[must_use]
    pub fn with_fit_intercept(mut self, fit_intercept: bool) -> Self {
        self.fit_intercept = fit_intercept;
        self
    }

    /// Set the per-class weighting strategy.
    ///
    /// Mirrors sklearn's `class_weight` parameter (default `None`). The expanded
    /// per-class weights (via [`compute_class_weight`]) scale only the
    /// gradient-derived part of each per-sample update
    /// (`_sgd_fast.pyx.tp:599-602,630`); the L2 shrink, L1 truncation and the
    /// one-class offset are left unscaled.
    #[must_use]
    pub fn with_class_weight(mut self, class_weight: ClassWeight<F>) -> Self {
        self.class_weight = class_weight;
        self
    }

    /// Set the averaged-SGD (ASGD) threshold.
    ///
    /// Mirrors sklearn's `average` parameter (default `False`,
    /// `_stochastic_gradient.py:1256`). `0` disables averaging (the plain SGD
    /// trajectory, byte-identical to the unaveraged kernel); `with_average(1)`
    /// is sklearn `average=True` (average from the first step); `with_average(N)`
    /// is sklearn `average=N` (begin averaging once the global step counter
    /// `t >= N`). The running mean of the post-update weights/intercept replaces
    /// the plain `coef_`/`intercept_` at fit-end when `average <= self.t_ - 1`
    /// (`_sgd_fast.pyx.tp:646-654`, `_stochastic_gradient.py:834-836`).
    #[must_use]
    pub fn with_average(mut self, average: usize) -> Self {
        self.average = average;
        self
    }

    /// Enable or disable early stopping on a held-out validation score.
    ///
    /// Mirrors sklearn's `early_stopping` parameter (default `False`,
    /// `_stochastic_gradient.py:114`, constraint `["boolean"]` at `:524`). When
    /// enabled, [`with_validation_fraction`](Self::with_validation_fraction) of
    /// the data is held out (stratified per class) and each One-vs-All binary
    /// subproblem's epoch-end convergence is driven by its validation accuracy
    /// rather than the training loss (`_sgd_fast.pyx.tp:678-687`).
    #[must_use]
    pub fn with_early_stopping(mut self, early_stopping: bool) -> Self {
        self.early_stopping = early_stopping;
        self
    }

    /// Set the fraction of the training data held out for early-stopping
    /// validation.
    ///
    /// Mirrors sklearn's `validation_fraction` parameter (default `0.1`,
    /// `_stochastic_gradient.py:115`, constraint
    /// `Interval(Real, 0, 1, closed="neither")` at `:525`). Only used when
    /// [`early_stopping`](Self::early_stopping) is `true`; validated to the open
    /// interval `(0, 1)` at fit time.
    #[must_use]
    pub fn with_validation_fraction(mut self, validation_fraction: F) -> Self {
        self.validation_fraction = validation_fraction;
        self
    }
}

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

/// Extract hyperparameter bundle from an `SGDClassifier`.
fn clf_hyper<F: Float>(clf: &SGDClassifier<F>) -> SGDHyper<F> {
    SGDHyper {
        learning_rate: clf.learning_rate,
        eta0: clf.eta0,
        alpha: clf.alpha,
        max_iter: clf.max_iter,
        tol: clf.tol,
        random_state: clf.random_state,
        power_t: clf.power_t,
        penalty: clf.penalty,
        l1_ratio: clf.l1_ratio,
        shuffle: clf.shuffle,
        n_iter_no_change: clf.n_iter_no_change,
        fit_intercept: clf.fit_intercept,
        one_class: false,
        average: clf.average,
        early_stopping: clf.early_stopping,
        validation_fraction: clf.validation_fraction,
    }
}

/// Internal hyperparameter bundle shared between Fit and PartialFit paths.
#[derive(Debug, Clone)]
struct SGDHyper<F> {
    learning_rate: LearningRateSchedule<F>,
    eta0: F,
    alpha: F,
    max_iter: usize,
    tol: F,
    random_state: Option<u64>,
    power_t: F,
    /// Regularization penalty (l2/l1/elasticnet).
    penalty: Penalty,
    /// Elastic-net mixing parameter (only meaningful for `ElasticNet`).
    l1_ratio: F,
    /// Whether to shuffle the sample order each epoch (`_sgd_fast.pyx.tp:579`).
    shuffle: bool,
    /// Number of consecutive non-improving epochs before convergence /
    /// adaptive-eta decay triggers (`_stochastic_gradient.py` default 5,
    /// `_sgd_fast.pyx.tp:698`).
    n_iter_no_change: usize,
    /// Whether to fit (update) the intercept each step. When `false` the
    /// intercept update is skipped (`_sgd_fast.pyx.tp:639-644`).
    fit_intercept: bool,
    /// Whether this is a one-class SVM fit. When `true` the (gated) intercept
    /// update gains the extra `- 2*eta*alpha` term, mirroring the `if one_class`
    /// branch in `_sgd_fast.pyx.tp:641-642`
    /// (`intercept_update -= 2. * eta * alpha`). `false` for the standard
    /// classifier/regressor paths, leaving their intercept update byte-identical.
    one_class: bool,
    /// Averaged-SGD (ASGD) threshold. `0` disables averaging (the default,
    /// byte-identical to the plain SGD trajectory). `N > 0` begins accumulating
    /// the running mean of the post-update weights/intercept once the global step
    /// counter `t >= N`, mirroring `if 0 < average <= t` (`_sgd_fast.pyx.tp:646`).
    /// sklearn `average=True` maps to `N = 1`; `average=N` maps to `N`.
    average: usize,
    /// Whether to use early stopping on a held-out validation score. When
    /// `true` the epoch-end convergence rule scores the current weights on the
    /// validation set (R^2 / accuracy) instead of the training loss
    /// (`_sgd_fast.pyx.tp:678-687`, `_stochastic_gradient.py:114`). The
    /// validation set itself is split off in the `Fit` path BEFORE the kernel
    /// and passed in separately; this flag only selects the score-based
    /// epoch-end branch. `false` (the default) leaves the training-loss
    /// convergence path byte-identical.
    early_stopping: bool,
    /// Fraction of the training data held out as the validation set when
    /// `early_stopping` is `true` (`_stochastic_gradient.py:115`, default `0.1`,
    /// constraint `Interval(Real, 0, 1, closed="neither")` at `:525`). Carried
    /// for validation/documentation; the actual split happens in the `Fit` path.
    validation_fraction: F,
}

/// Train a single binary classifier via SGD, updating `weights` and
/// `intercept` in place. `y_binary` must be in `{-1, +1}`.
///
/// `sample_w[i]` is the per-sample weight `class_weight_i * sample_weight_i`
/// for sample `i` (`_sgd_fast.pyx.tp:599-602,630`). It scales ONLY the
/// gradient-derived part of the update (`update *= class_weight * sample_weight`
/// at `:630`); the L2 shrink (`:632-635`), the L1 truncation (`:656-658`) and
/// the one-class `-2*eta*alpha` offset (`:642`) are left unscaled. An all-ones
/// `sample_w` (the default `fit` path) reproduces the byte-identical unweighted
/// behaviour. `sample_w.len()` must equal `x.nrows()`.
///
/// Returns the cumulative loss and the step counter after training.
#[allow(
    clippy::too_many_arguments,
    reason = "threads the per-sample weight vector"
)]
fn train_binary_sgd<F, L>(
    x: &Array2<F>,
    y_binary: &Array1<F>,
    weights: &mut Array1<F>,
    intercept: &mut F,
    loss_fn: &L,
    hyper: &SGDHyper<F>,
    initial_t: usize,
    sample_w: &[F],
    val_set: Option<(&Array2<F>, &Array1<F>)>,
) -> (F, usize)
where
    F: Float + ScalarOperand + Send + Sync + 'static,
    L: Loss<F>,
{
    let n_samples = x.nrows();
    let n_features = x.ncols();
    let mut t = initial_t;
    // Epoch-end convergence state, mirroring `_sgd_fast.pyx.tp:525,532-534`:
    // `best_loss = INFINITY`, `no_improvement_count = 0`. `current_eta` carries
    // the adaptive-schedule eta (`eta = eta / 5` decay, `:700`). When early
    // stopping is active `best_score = -INFINITY` instead (higher score is
    // better — `_sgd_fast.pyx.tp:533`).
    let mut best_loss = F::infinity();
    let mut best_score = F::neg_infinity();
    // Early stopping uses the score branch only when a validation set was split
    // off in the `Fit` path (`val_set.is_some()`). The relabeled `{-1,+1}`
    // binary validation target is scored with accuracy (`binary_accuracy`,
    // `_stochastic_gradient.py:451-454,79`).
    let early_stopping = hyper.early_stopping && val_set.is_some();
    let mut current_eta = hyper.eta0;
    let mut no_improve_count: usize = 0;
    // `tol = None` upstream becomes `-INFINITY`, disabling the stop rule
    // (`tol > -INFINITY` is false, `:690`). ferrolearn encodes that as a
    // finite/-inf `tol`; the criterion is active iff `tol > -inf`.
    let tol_active = hyper.tol > F::neg_infinity();
    let max_dloss = max_dloss::<F>();
    let mut indices: Vec<usize> = (0..n_samples).collect();
    // `optimal_init` (the `optimal` schedule's t0 offset) depends on the loss
    // and alpha, so it is computed once per fit, before the epoch loop
    // (`_sgd_fast.pyx.tp:565-570`).
    let opt_init = optimal_init(loss_fn, hyper.alpha);

    // Effective l1_ratio from the penalty (`_sgd_fast.pyx.tp:558-561`):
    // `L2 -> 0.0`, `L1 -> 1.0`, `ElasticNet -> user l1_ratio`.
    let eff = effective_l1_ratio(hyper.penalty, hyper.l1_ratio);
    let apply_l1 = matches!(hyper.penalty, Penalty::L1 | Penalty::ElasticNet);
    // Tsuruoka cumulative-penalty state. `u` (scalar) accumulates the total L1
    // penalty applied so far; `q` (per-feature) records how much penalty has
    // actually been applied to each weight. Both persist for the WHOLE fit —
    // allocated once before the epoch loop, mirroring `q = np.zeros(...)` and
    // `u = 0.0` allocated once per `_plain_sgd` call (`_sgd_fast.pyx.tp:551-556`).
    let mut u = F::zero();
    let mut q: Array1<F> = Array1::zeros(n_features);

    // Averaged-SGD (ASGD) accumulators (`_sgd_fast.pyx.tp:646-654`). When
    // `hyper.average > 0`, once the global step `t >= average` we maintain the
    // running mean of the POST-update weights/intercept. This is the DIRECT
    // running-mean form of sklearn's lazy `w.add_average` (a wscale optimization
    // that is mathematically identical for plain arrays): with
    // `num_iter = t - average + 1` (= 1 at the first averaged step) and
    // `mu = 1/num_iter`, `avg += (current - avg) * mu`. The accumulator is a
    // PASSIVE observer — it never feeds back into the live `weights`/`intercept`.
    let mut average_coef: Array1<F> = Array1::zeros(n_features);
    let mut average_intercept = F::zero();

    // Build the RNG for shuffling.
    let mut rng = match hyper.random_state {
        Some(seed) => rand::rngs::StdRng::seed_from_u64(seed),
        None => rand::rngs::StdRng::from_os_rng(),
    };

    let mut total_loss = F::zero();

    for _epoch in 0..hyper.max_iter {
        // sklearn shuffles the sample order each epoch only when `shuffle` is
        // set (`_sgd_fast.pyx.tp:579-580`: `if shuffle: dataset.shuffle(seed)`).
        // With `shuffle == false` the indices stay `0..n-1` every epoch, exactly
        // matching sklearn's no-shuffle index order (`:581 for i in range(n)`).
        if hyper.shuffle {
            indices.shuffle(&mut rng);
        }
        let mut epoch_loss = F::zero();

        for &i in &indices {
            t += 1;

            let eta = match hyper.learning_rate {
                LearningRateSchedule::Adaptive => current_eta,
                _ => compute_lr(
                    &hyper.learning_rate,
                    hyper.eta0,
                    hyper.alpha,
                    hyper.power_t,
                    opt_init,
                    t,
                ),
            };

            // Compute prediction: w^T x_i + b.
            let mut y_pred = *intercept;
            let xi = x.row(i);
            for j in 0..n_features {
                y_pred = y_pred + weights[j] * xi[j];
            }

            // Clip the gradient to `[-MAX_DLOSS, MAX_DLOSS]` before forming the
            // update, matching `_sgd_fast.pyx.tp:613-620`.
            let grad = loss_fn
                .gradient(y_binary[i], y_pred)
                .max(-max_dloss)
                .min(max_dloss);
            // Per-sample weight scaling: `update *= class_weight * sample_weight`
            // (`_sgd_fast.pyx.tp:630`). `update = -eta*dloss`, so scaling the
            // update by `w_i` is equivalent to scaling the (clipped) gradient
            // `dloss` by `w_i` BEFORE forming both the weight data term and the
            // (gated) intercept gradient term. This multiplies ONLY the
            // gradient-derived part; the L2 shrink, L1 truncation and the
            // one-class offset below are unaffected. `g` is the scaled gradient.
            let g = grad * sample_w[i];
            // `sumloss` is the SUM (not mean) of per-sample losses over the
            // epoch (`_sgd_fast.pyx.tp:597`), computed from the UNWEIGHTED loss
            // (the weight only multiplies `update`, not `loss.loss(y, p)`).
            epoch_loss = epoch_loss + loss_fn.loss(y_binary[i], y_pred);

            // L2 shrink: scale the whole weight vector by the CLAMPED factor
            // `max(0, 1 - (1-eff)*eta*alpha)` BEFORE the gradient add, mirroring
            // `w.scale(max(0, 1 - (1-l1_ratio)*eta*alpha))`
            // (`_sgd_fast.pyx.tp:632-635`). For pure L2 (`eff=0`) this is
            // `max(0, 1-eta*alpha)`; for L1 (`eff=1`) it is `max(0, 1) = 1`
            // (no L2 shrink); for elasticnet the `(1-eff)` weakens the L2 part.
            let shrink = (F::one() - (F::one() - eff) * eta * hyper.alpha).max(F::zero());
            for j in 0..n_features {
                weights[j] = weights[j] * shrink;
            }
            // Gradient add `w.add(x, update)` with the scaled `update = -eta*g`
            // (`_sgd_fast.pyx.tp:637-638`); `g` is the sample-weighted gradient.
            for j in 0..n_features {
                weights[j] = weights[j] - eta * g * xi[j];
            }
            // The intercept update is gated on `fit_intercept` and is NOT
            // regularized (`intercept_decay=1`, `_sgd_fast.pyx.tp:639-644`:
            // `if fit_intercept == 1: intercept_update = update; if one_class:
            // intercept_update -= 2.*eta*alpha; intercept += intercept_update *
            // intercept_decay`). `update = -eta*g` is the SCALED update (sklearn
            // sets `intercept_update = update` at `:640`, after `update *=
            // class_weight*sample_weight` at `:630`), so the standard path is
            // `intercept -= eta*g`. For the one-class SVM the extra
            // `- 2*eta*alpha` term is added (`:641-642`) and is NOT scaled by the
            // sample weight. When `fit_intercept` is false the intercept is never
            // modified and stays at its init value (`0` clf/reg, `1` one-class).
            if hyper.fit_intercept {
                let two = cst::<F>(2.0);
                let mut intercept_update = -eta * g;
                if hyper.one_class {
                    intercept_update = intercept_update - two * eta * hyper.alpha;
                }
                *intercept = *intercept + intercept_update;
            }

            // L1 cumulative penalty (Tsuruoka truncated gradient), applied AFTER
            // the gradient add only for L1/ElasticNet (`_sgd_fast.pyx.tp:656-658`,
            // `l1penalty` at `:750-778` with `wscale = 1`).
            if apply_l1 {
                u = u + eff * eta * hyper.alpha;
                for j in 0..n_features {
                    let z = weights[j];
                    if weights[j] > F::zero() {
                        weights[j] = (weights[j] - (u + q[j])).max(F::zero());
                    } else if weights[j] < F::zero() {
                        weights[j] = (weights[j] + (u - q[j])).min(F::zero());
                    }
                    q[j] = q[j] + (weights[j] - z);
                }
            }

            // ASGD running-mean accumulation (`_sgd_fast.pyx.tp:646-654`:
            // `if 0 < average <= t: w.add_average(..., t - average + 1);
            // average_intercept += (intercept - average_intercept) /
            // (t - average + 1)`). Performed AFTER the weight/intercept update +
            // L1 truncation, so `weights`/`intercept` hold their final post-step
            // values for this sample. `t` here is the SAME 1-based global step
            // the schedule used above. `num_iter = t - average + 1` is `>= 1`
            // whenever `t >= average`.
            if hyper.average > 0 && t >= hyper.average {
                let num_iter = t - hyper.average + 1;
                let num_iter_f = F::from(num_iter).unwrap_or_else(F::one);
                let mu = F::one() / num_iter_f;
                for j in 0..n_features {
                    average_coef[j] = average_coef[j] + (weights[j] - average_coef[j]) * mu;
                }
                average_intercept = average_intercept + (*intercept - average_intercept) * mu;
            }
        }

        // `epoch_loss` is now the epoch `sumloss` (no mean division).
        total_loss = epoch_loss;

        // Epoch-end stop rule (`_sgd_fast.pyx.tp:678-707`). When early stopping
        // is active, score the CURRENT weights/intercept on the held-out
        // validation set (binary accuracy of the relabeled `{-1,+1}` target,
        // `_stochastic_gradient.py:79`) and run the score-based branch
        // (`best_score` init `-inf`, higher is better, `:678-687`); otherwise
        // the training-loss branch (`sumloss` vs `best_loss`, `:688-695`).
        let should_break = if let (true, Some((x_val, y_val))) = (early_stopping, val_set) {
            let score = binary_accuracy(weights, *intercept, x_val, y_val);
            convergence_tail_score(
                score,
                &mut best_score,
                &mut no_improve_count,
                &mut current_eta,
                tol_active,
                hyper.tol,
                hyper.n_iter_no_change,
                matches!(hyper.learning_rate, LearningRateSchedule::Adaptive),
            )
        } else {
            convergence_tail(
                epoch_loss,
                &mut best_loss,
                &mut no_improve_count,
                &mut current_eta,
                tol_active,
                hyper.tol,
                n_samples,
                hyper.n_iter_no_change,
                matches!(hyper.learning_rate, LearningRateSchedule::Adaptive),
            )
        };
        if should_break {
            break;
        }
    }

    // ASGD finalize: select the averaged weights/intercept when averaging was
    // enabled AND enough steps were taken (`_stochastic_gradient.py:834-836`:
    // `if self.average > 0: if self.average <= self.t_ - 1: coef_ =
    // average_coef`). Here `t` is the returned step counter `= n_iter_ *
    // n_samples` (`initial_t = 0` on the full-fit path), which equals sklearn's
    // `self.t_ - 1` (sklearn inits `self.t_ = 1`, then `self.t_ += n_iter_ *
    // n_samples`). So `average <= self.t_ - 1` maps to `hyper.average <= t`.
    if hyper.average > 0 && hyper.average <= t {
        for j in 0..n_features {
            weights[j] = average_coef[j];
        }
        *intercept = average_intercept;
    }

    (total_loss, t)
}

/// Fitted SGD classifier.
///
/// Holds the learned weight vectors and intercepts. For binary problems
/// there is a single weight vector; for multiclass problems there is one
/// per class (one-vs-all).
///
/// Implements [`Predict`] and [`PartialFit`] to support both inference and
/// online learning.
#[derive(Debug, Clone)]
pub struct FittedSGDClassifier<F> {
    /// Weight matrix: one row per binary sub-problem.
    /// Binary: shape `(1, n_features)`, multiclass: `(n_classes, n_features)`.
    weight_matrix: Vec<Array1<F>>,
    /// Intercept vector, one per sub-problem.
    intercepts: Vec<F>,
    /// Sorted unique class labels.
    classes: Vec<usize>,
    /// Number of features the model was trained on.
    n_features: usize,
    /// The loss function used during training.
    loss: ClassifierLoss,
    /// Hyperparameters for continued training via `partial_fit`.
    hyper: SGDHyper<F>,
    /// Global step counter across all partial_fit calls.
    t: usize,
}

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

    /// Fit the SGD classifier on the given data.
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if `x` and `y` have mismatched
    /// sample counts.
    /// Returns [`FerroError::InsufficientSamples`] if fewer than 2 classes
    /// are present.
    /// Returns [`FerroError::InvalidParameter`] if `eta0` or `alpha` are
    /// not positive.
    fn fit(&self, x: &Array2<F>, y: &Array1<usize>) -> Result<FittedSGDClassifier<F>, FerroError> {
        // Delegate to the sample-weighted path with a uniform `ones(n)` weight
        // vector, so the default `fit` behaviour is byte-identical to the
        // unweighted kernel (`_check_sample_weight` returns `ones` when
        // `sample_weight=None`, `_stochastic_gradient.py:627`).
        let sample_weight = Array1::<F>::from_elem(x.nrows(), F::one());
        self.fit_with_sample_weight(x, y, &sample_weight)
    }
}

impl<F: Float + Send + Sync + ScalarOperand + 'static> SGDClassifier<F> {
    /// Fit the SGD classifier with explicit per-sample weights.
    ///
    /// Mirrors `SGDClassifier.fit(X, y, sample_weight=...)`. The per-sample
    /// weight scales ONLY the gradient-derived part of each update
    /// (`update *= class_weight * sample_weight`, `_sgd_fast.pyx.tp:630`); the
    /// L2 shrink, L1 truncation and one-class offset are unscaled. The
    /// `class_weight` field (via [`compute_class_weight`],
    /// `_stochastic_gradient.py:624`) is combined multiplicatively per sample.
    ///
    /// [`Fit::fit`] delegates here with a uniform `ones(n)` weight vector.
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if `x`/`y`/`sample_weight` have
    /// mismatched sample counts.
    /// Returns [`FerroError::InsufficientSamples`] if fewer than 2 classes
    /// are present.
    /// Returns [`FerroError::InvalidParameter`] if `eta0` or `alpha` are
    /// invalid.
    pub fn fit_with_sample_weight(
        &self,
        x: &Array2<F>,
        y: &Array1<usize>,
        sample_weight: &Array1<F>,
    ) -> Result<FittedSGDClassifier<F>, FerroError> {
        validate_clf_params(
            x,
            y,
            &self.learning_rate,
            self.eta0,
            self.alpha,
            self.l1_ratio,
            self.validation_fraction,
        )?;

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

        // Non-finite input validation (#2263). sklearn `SGDClassifier.fit`
        // -> `self._validate_data(X, y, ...)` (`_stochastic_gradient.py:1476`)
        // keeps the default `force_all_finite=True`, so `check_array` rejects any
        // NaN or +/-inf in X with a `ValueError("Input X contains NaN.")` /
        // `"... contains infinity ..."` BEFORE the SGD kernel. sklearn also
        // validates `sample_weight` via `_check_sample_weight` (default
        // `force_all_finite=True`, `_stochastic_gradient.py:1501`). `y` is
        // `Array1<usize>` (integer class labels), finite by type, so only X +
        // sample_weight need the runtime check. `.iter().any(|v| !v.is_finite())`
        // rejects both NaN and Inf (bounds-safe, no panic, R-CODE-2); the finite
        // path is byte-identical. This is the SGDClassifier fit entry (`Fit::fit`
        // delegates here with unit weights).
        if x.iter().any(|v| !v.is_finite()) {
            return Err(FerroError::InvalidParameter {
                name: "X".into(),
                reason: "Input X contains NaN or infinity.".into(),
            });
        }
        if sample_weight.iter().any(|v| !v.is_finite()) {
            return Err(FerroError::InvalidParameter {
                name: "sample_weight".into(),
                reason: "Input sample_weight contains NaN or infinity.".into(),
            });
        }

        let n_features = x.ncols();
        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: "SGDClassifier requires at least 2 distinct classes".into(),
            });
        }

        let hyper = clf_hyper(self);
        let loss_enum = self.loss;

        // Expanded per-class weights (`_stochastic_gradient.py:624`), aligned to
        // the sorted `classes` (= sklearn `classes_`).
        let expanded = compute_class_weight(&self.class_weight, &classes, &y.to_vec());
        let sw = sample_weight.to_vec();

        let (weight_matrix, intercepts, t) = fit_ova(
            x, y, &classes, n_features, &loss_enum, &hyper, 0, &expanded, &sw,
        )?;

        Ok(FittedSGDClassifier {
            weight_matrix,
            intercepts,
            classes,
            n_features,
            loss: loss_enum,
            hyper,
            t,
        })
    }
}

/// Whether a learning-rate schedule requires `eta0 > 0`.
///
/// Mirrors sklearn `_more_validate_params` (`_stochastic_gradient.py:149-153`):
/// `eta0 > 0` is enforced only for `constant`/`invscaling`/`adaptive`; the
/// `optimal` schedule derives its own initial rate and accepts `eta0 == 0`.
fn schedule_requires_eta0<F: Float>(schedule: &LearningRateSchedule<F>) -> bool {
    matches!(
        schedule,
        LearningRateSchedule::Constant
            | LearningRateSchedule::InvScaling
            | LearningRateSchedule::Adaptive
    )
}

/// Validate classifier input shapes and parameters.
/// Validate `validation_fraction` to the OPEN interval `(0, 1)`.
///
/// Mirrors sklearn `_parameter_constraints["validation_fraction"]`
/// (`_stochastic_gradient.py:525`, `Interval(Real, 0, 1, closed="neither")`):
/// the bounds are both EXCLUSIVE, so `0.0` and `1.0` are invalid. sklearn
/// validates this unconditionally (it is a constructor constraint), independent
/// of `early_stopping`.
fn validate_validation_fraction<F: Float>(validation_fraction: F) -> Result<(), FerroError> {
    if validation_fraction <= F::zero() || validation_fraction >= F::one() {
        return Err(FerroError::InvalidParameter {
            name: "validation_fraction".into(),
            reason: "must be in the open interval (0, 1)".into(),
        });
    }
    Ok(())
}

#[allow(
    clippy::too_many_arguments,
    reason = "threads each validated parameter"
)]
fn validate_clf_params<F: Float>(
    x: &Array2<F>,
    y: &Array1<usize>,
    schedule: &LearningRateSchedule<F>,
    eta0: F,
    alpha: F,
    l1_ratio: F,
    validation_fraction: F,
) -> Result<(), FerroError> {
    let n_samples = x.nrows();
    if n_samples != y.len() {
        return Err(FerroError::ShapeMismatch {
            expected: vec![n_samples],
            actual: vec![y.len()],
            context: "y length must match number of samples in X".into(),
        });
    }
    if n_samples == 0 {
        return Err(FerroError::InsufficientSamples {
            required: 1,
            actual: 0,
            context: "SGDClassifier requires at least one sample".into(),
        });
    }
    if schedule_requires_eta0(schedule) && eta0 <= F::zero() {
        return Err(FerroError::InvalidParameter {
            name: "eta0".into(),
            reason: "must be positive".into(),
        });
    }
    if alpha < F::zero() {
        return Err(FerroError::InvalidParameter {
            name: "alpha".into(),
            reason: "must be non-negative".into(),
        });
    }
    if l1_ratio < F::zero() || l1_ratio > F::one() {
        return Err(FerroError::InvalidParameter {
            name: "l1_ratio".into(),
            reason: "must be in the range [0, 1]".into(),
        });
    }
    validate_validation_fraction(validation_fraction)?;
    Ok(())
}

/// Result type for one-vs-all training: (weight_matrix, intercepts, step_counter).
type OvaResult<F> = (Vec<Array1<F>>, Vec<F>, usize);

/// Number of validation samples for an early-stopping split: `ceil` of
/// `validation_fraction * n` clamped to `[1, n-1]` so that BOTH the train and
/// the validation subset are non-empty.
///
/// sklearn delegates the count to `ShuffleSplit`/`StratifiedShuffleSplit`, which
/// use `ceil(test_size * n)` and raise if either subset is empty
/// (`_stochastic_gradient.py:295-307`). The exact sample SELECTION is numpy-RNG
/// coupled and not cross-impl reproducible (the same barrier as `shuffle`); only
/// the count + non-emptiness are reproduced here.
fn validation_count<F: Float>(validation_fraction: F, n: usize) -> usize {
    let n_f = F::from(n).unwrap_or_else(F::zero);
    let raw = (validation_fraction * n_f).ceil();
    let n_val = raw.to_usize().unwrap_or(1).max(1);
    n_val.min(n.saturating_sub(1))
}

/// Build a seeded, optionally stratified train/validation index partition for
/// early stopping.
///
/// Returns `(train_idx, val_idx)`. The first `n_val` entries of a seeded random
/// permutation form the validation set; for the classifier (`stratify = Some`)
/// the permutation is built per class so the validation set is proportional per
/// class, mirroring sklearn's `StratifiedShuffleSplit`
/// (`_stochastic_gradient.py:280-287`); for the regressor a plain `ShuffleSplit`
/// permutation is used. The RNG is `StdRng::seed_from_u64(random_state ?? 0)`.
///
/// The SELECTION is intentionally NOT identical to sklearn (numpy's
/// Mersenne-Twister permutation differs from `StdRng`); only the deterministic
/// contract — a valid, seeded, stratified-for-classifier, non-empty split — is
/// reproduced. Returns `None` if either subset would be empty (sklearn raises a
/// `ValueError`, `_stochastic_gradient.py:295-307`).
fn make_validation_split<F: Float>(
    n: usize,
    validation_fraction: F,
    random_state: Option<u64>,
    stratify: Option<&[usize]>,
) -> Option<(Vec<usize>, Vec<usize>)> {
    if n < 2 {
        return None;
    }
    let n_val = validation_count(validation_fraction, n);
    if n_val == 0 || n_val >= n {
        return None;
    }
    let mut rng = rand::rngs::StdRng::seed_from_u64(random_state.unwrap_or(0));

    let mut val_mask = vec![false; n];
    match stratify {
        Some(labels) => {
            // Per-class proportional hold-out (`StratifiedShuffleSplit`). For each
            // class, shuffle its member indices and take `round(frac * count)`
            // (at least 1 when the class has >= 2 members) into validation.
            let mut classes: Vec<usize> = labels.to_vec();
            classes.sort_unstable();
            classes.dedup();
            for &c in &classes {
                let mut members: Vec<usize> = (0..n).filter(|&i| labels[i] == c).collect();
                members.shuffle(&mut rng);
                let count = members.len();
                let frac_f = F::from(count).unwrap_or_else(F::zero) * validation_fraction;
                let mut take = frac_f.round().to_usize().unwrap_or(0);
                if take == 0 && count >= 2 {
                    take = 1;
                }
                take = take.min(count.saturating_sub(1)).min(count);
                for &idx in members.iter().take(take) {
                    val_mask[idx] = true;
                }
            }
        }
        None => {
            // Plain shuffle hold-out (`ShuffleSplit`): first `n_val` of a seeded
            // permutation.
            let mut perm: Vec<usize> = (0..n).collect();
            perm.shuffle(&mut rng);
            for &idx in perm.iter().take(n_val) {
                val_mask[idx] = true;
            }
        }
    }

    let val_idx: Vec<usize> = (0..n).filter(|&i| val_mask[i]).collect();
    let train_idx: Vec<usize> = (0..n).filter(|&i| !val_mask[i]).collect();
    if val_idx.is_empty() || train_idx.is_empty() {
        return None;
    }
    Some((train_idx, val_idx))
}

/// Gather the rows of `x` indexed by `idx` into a fresh `Array2`.
fn gather_rows<F: Float>(x: &Array2<F>, idx: &[usize]) -> Array2<F> {
    let n_features = x.ncols();
    let mut out = Array2::<F>::zeros((idx.len(), n_features));
    for (r, &i) in idx.iter().enumerate() {
        let src = x.row(i);
        for j in 0..n_features {
            out[[r, j]] = src[j];
        }
    }
    out
}

/// Gather the entries of `v` indexed by `idx` into a fresh `Array1`.
fn gather<F: Float>(v: &Array1<F>, idx: &[usize]) -> Array1<F> {
    Array1::from_iter(idx.iter().map(|&i| v[i]))
}

/// Train one-vs-all binary classifiers, returning per-class weights, intercepts,
/// and the cumulative step counter.
///
/// `expanded_class_weight[k]` is the weight of `classes[k]` from
/// [`compute_class_weight`] (`_stochastic_gradient.py:624`). `sample_weight[i]`
/// is the user per-sample weight. For each binary subproblem the per-sample
/// weight passed to the kernel is `class_weight_for_sample(i) * sample_weight[i]`
/// where `class_weight_for_sample(i)` is `pos_weight` if sample `i` is the
/// positive class else `neg_weight`, with the sklearn OvA mapping:
/// binary (`_fit_binary`, `:765-766`) `pos = expanded[1]`, `neg = expanded[0]`;
/// multiclass class `k` (`_fit_multiclass`, `:816`) `pos = expanded[k]`,
/// `neg = 1.0`.
#[allow(clippy::too_many_arguments, reason = "threads class + sample weights")]
fn fit_ova<F: Float + Send + Sync + ScalarOperand + 'static>(
    x: &Array2<F>,
    y: &Array1<usize>,
    classes: &[usize],
    n_features: usize,
    loss_enum: &ClassifierLoss,
    hyper: &SGDHyper<F>,
    initial_t: usize,
    expanded_class_weight: &[F],
    sample_weight: &[F],
) -> Result<OvaResult<F>, FerroError> {
    let n_classes = classes.len();
    let mut weight_matrix: Vec<Array1<F>> = Vec::with_capacity(n_classes);
    let mut intercepts: Vec<F> = Vec::with_capacity(n_classes);
    let mut global_t = initial_t;

    // Early-stopping validation split. Computed ONCE over the full multiclass
    // labels (so the hold-out is stratified per class and SHARED by every OvA
    // subproblem, exactly as sklearn precomputes the mask in `_fit_multiclass`,
    // `_stochastic_gradient.py:796`, and reuses it for each binary fit). The
    // split is stratified (`StratifiedShuffleSplit`, `:280-281`). When the split
    // is infeasible (too few samples) it returns `None`, and an empty validation
    // set raises (`:295-307`).
    let split = if hyper.early_stopping {
        match make_validation_split(
            x.nrows(),
            hyper.validation_fraction,
            hyper.random_state,
            Some(&y.to_vec()),
        ) {
            Some(s) => Some(s),
            None => {
                return Err(FerroError::InvalidParameter {
                    name: "validation_fraction".into(),
                    reason: "early_stopping split led to an empty train or validation set; \
                             increase the number of samples or change validation_fraction"
                        .into(),
                });
            }
        }
    } else {
        None
    };

    // Closure: run one OvA binary subproblem (relabel + per-sample weights +
    // optional validation slice + kernel call). `pos`/`neg` are the class-weight
    // mappings for this subproblem.
    let run_subproblem =
        |cls: usize, pos_weight: F, neg_weight: F, w: &mut Array1<F>, b: &mut F, t0: usize| {
            let y_binary: Array1<F> =
                y.mapv(|label| if label == cls { F::one() } else { -F::one() });
            // Per-sample weight = class_weight_for_sample(i) * sample_weight[i]
            // (`_sgd_fast.pyx.tp:599-602,630`). `y_binary[i] > 0` selects pos.
            let sample_w_full: Vec<F> = (0..x.nrows())
                .map(|i| {
                    let cw = if y_binary[i] > F::zero() {
                        pos_weight
                    } else {
                        neg_weight
                    };
                    cw * sample_weight[i]
                })
                .collect();

            if let Some((train_idx, val_idx)) = &split {
                // Train on the train subset only; score on the held-out
                // validation subset (relabeled `{-1,+1}` target, accuracy —
                // `_stochastic_gradient.py:451-454,79`). sklearn does NOT refit
                // on the full data; the train-subset weights are final
                // (`_ValidationScoreCallback` only reads weights, never writes
                // back to the live fit).
                let x_tr = gather_rows(x, train_idx);
                let y_tr = gather(&y_binary, train_idx);
                let sw_tr: Vec<F> = train_idx.iter().map(|&i| sample_w_full[i]).collect();
                let x_val = gather_rows(x, val_idx);
                let y_val = gather(&y_binary, val_idx);
                dispatch_train_binary(
                    &x_tr,
                    &y_tr,
                    w,
                    b,
                    loss_enum,
                    hyper,
                    t0,
                    &sw_tr,
                    Some((&x_val, &y_val)),
                )
            } else {
                dispatch_train_binary(
                    x,
                    &y_binary,
                    w,
                    b,
                    loss_enum,
                    hyper,
                    t0,
                    &sample_w_full,
                    None,
                )
            }
        };

    if n_classes == 2 {
        // Single binary problem: class[0] -> -1, class[1] -> +1.
        // OvA weight mapping (`_fit_binary`, `_stochastic_gradient.py:765-766`):
        // pos_weight = expanded[1], neg_weight = expanded[0].
        let pos_weight = expanded_class_weight[1];
        let neg_weight = expanded_class_weight[0];
        let mut w = Array1::<F>::zeros(n_features);
        let mut b = F::zero();
        let (_, t) = run_subproblem(classes[1], pos_weight, neg_weight, &mut w, &mut b, global_t);
        global_t = t;
        weight_matrix.push(w);
        intercepts.push(b);
    } else {
        // One-vs-all: one binary problem per class. Multiclass mapping
        // (`_fit_multiclass`, `_stochastic_gradient.py:816`): for class k,
        // pos_weight = expanded[k], neg_weight = 1.0.
        for (k, &cls) in classes.iter().enumerate() {
            let pos_weight = expanded_class_weight[k];
            let neg_weight = F::one();
            let mut w = Array1::<F>::zeros(n_features);
            let mut b = F::zero();
            let (_, t) = run_subproblem(cls, pos_weight, neg_weight, &mut w, &mut b, global_t);
            global_t = t;
            weight_matrix.push(w);
            intercepts.push(b);
        }
    }

    Ok((weight_matrix, intercepts, global_t))
}

/// Train one-vs-all using existing weight vectors (for partial_fit).
#[allow(clippy::too_many_arguments)]
fn partial_fit_ova<F: Float + Send + Sync + ScalarOperand + 'static>(
    x: &Array2<F>,
    y: &Array1<usize>,
    classes: &[usize],
    weight_matrix: &mut [Array1<F>],
    intercepts: &mut [F],
    loss_enum: &ClassifierLoss,
    hyper: &SGDHyper<F>,
    initial_t: usize,
) -> usize {
    let n_classes = classes.len();
    let mut global_t = initial_t;
    // `partial_fit` does not (yet) carry `class_weight`/`sample_weight`, so the
    // per-sample weight is uniform `1.0` — the all-ones path is byte-identical to
    // the pre-weighting kernel (`update *= 1*1`, `_sgd_fast.pyx.tp:630`).
    let sample_w: Vec<F> = vec![F::one(); x.nrows()];

    if n_classes == 2 {
        let y_binary: Array1<F> = y.mapv(|label| {
            if label == classes[1] {
                F::one()
            } else {
                -F::one()
            }
        });
        let (_, t) = dispatch_train_binary(
            x,
            &y_binary,
            &mut weight_matrix[0],
            &mut intercepts[0],
            loss_enum,
            hyper,
            global_t,
            &sample_w,
            None,
        );
        global_t = t;
    } else {
        for (idx, &cls) in classes.iter().enumerate() {
            let y_binary: Array1<F> =
                y.mapv(|label| if label == cls { F::one() } else { -F::one() });
            let (_, t) = dispatch_train_binary(
                x,
                &y_binary,
                &mut weight_matrix[idx],
                &mut intercepts[idx],
                loss_enum,
                hyper,
                global_t,
                &sample_w,
                None,
            );
            global_t = t;
        }
    }

    global_t
}

/// Dispatch to the appropriate typed loss training function.
///
/// `sample_w[i] = class_weight_i * sample_weight_i` is the per-sample weight
/// (`_sgd_fast.pyx.tp:599-602,630`), forwarded verbatim to the kernel.
#[allow(
    clippy::too_many_arguments,
    reason = "threads the per-sample weight vector"
)]
fn dispatch_train_binary<F: Float + Send + Sync + ScalarOperand + 'static>(
    x: &Array2<F>,
    y_binary: &Array1<F>,
    w: &mut Array1<F>,
    b: &mut F,
    loss_enum: &ClassifierLoss,
    hyper: &SGDHyper<F>,
    initial_t: usize,
    sample_w: &[F],
    val_set: Option<(&Array2<F>, &Array1<F>)>,
) -> (F, usize) {
    match loss_enum {
        ClassifierLoss::Hinge => train_binary_sgd(
            x, y_binary, w, b, &Hinge, hyper, initial_t, sample_w, val_set,
        ),
        ClassifierLoss::SquaredHinge => train_binary_sgd(
            x,
            y_binary,
            w,
            b,
            &SquaredHinge,
            hyper,
            initial_t,
            sample_w,
            val_set,
        ),
        ClassifierLoss::Perceptron => train_binary_sgd(
            x,
            y_binary,
            w,
            b,
            &Perceptron,
            hyper,
            initial_t,
            sample_w,
            val_set,
        ),
        ClassifierLoss::Log => train_binary_sgd(
            x, y_binary, w, b, &LogLoss, hyper, initial_t, sample_w, val_set,
        ),
        ClassifierLoss::SquaredError => train_binary_sgd(
            x,
            y_binary,
            w,
            b,
            &SquaredError,
            hyper,
            initial_t,
            sample_w,
            val_set,
        ),
        ClassifierLoss::ModifiedHuber => train_binary_sgd(
            x,
            y_binary,
            w,
            b,
            &ModifiedHuber,
            hyper,
            initial_t,
            sample_w,
            val_set,
        ),
    }
}

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

    /// Predict class labels for the given feature matrix.
    ///
    /// For binary classification, uses `sign(w^T x + b)`.
    /// For multiclass, returns the class whose one-vs-all score is highest.
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if the number of features
    /// does not match the fitted model.
    fn predict(&self, x: &Array2<F>) -> Result<Array1<usize>, FerroError> {
        let n_features = x.ncols();
        if n_features != self.n_features {
            return Err(FerroError::ShapeMismatch {
                expected: vec![self.n_features],
                actual: vec![n_features],
                context: "number of features must match fitted model".into(),
            });
        }

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

        if self.classes.len() == 2 {
            // Binary: single weight vector.
            let scores = x.dot(&self.weight_matrix[0]) + self.intercepts[0];
            for i in 0..n_samples {
                predictions[i] = if scores[i] >= F::zero() {
                    self.classes[1]
                } else {
                    self.classes[0]
                };
            }
        } else {
            // Multiclass: one-vs-all, pick highest score.
            for i in 0..n_samples {
                let xi = x.row(i);
                let mut best_class = 0;
                let mut best_score = F::neg_infinity();
                for (c, w) in self.weight_matrix.iter().enumerate() {
                    let score = xi.dot(w) + self.intercepts[c];
                    if score > best_score {
                        best_score = score;
                        best_class = c;
                    }
                }
                predictions[i] = self.classes[best_class];
            }
        }

        Ok(predictions)
    }
}

impl<F: Float + Send + Sync + ScalarOperand + 'static> PartialFit<Array2<F>, Array1<usize>>
    for FittedSGDClassifier<F>
{
    type FitResult = FittedSGDClassifier<F>;
    type Error = FerroError;

    /// Incrementally train the classifier on a new batch of data.
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if `x` and `y` have mismatched
    /// sizes or `x` has the wrong number of features.
    fn partial_fit(
        mut self,
        x: &Array2<F>,
        y: &Array1<usize>,
    ) -> Result<FittedSGDClassifier<F>, FerroError> {
        let n_samples = x.nrows();
        if n_samples != y.len() {
            return Err(FerroError::ShapeMismatch {
                expected: vec![n_samples],
                actual: vec![y.len()],
                context: "y length must match number of samples in X".into(),
            });
        }
        if x.ncols() != self.n_features {
            return Err(FerroError::ShapeMismatch {
                expected: vec![self.n_features],
                actual: vec![x.ncols()],
                context: "number of features must match fitted model".into(),
            });
        }

        // sklearn `SGDClassifier.partial_fit` validates X through
        // `_validate_data(force_all_finite=True)` (`_stochastic_gradient.py:596`),
        // raising `ValueError("Input X contains NaN.")` for any non-finite X
        // BEFORE the kernel. `y` is integer labels (always finite by type).
        // Mirrors the #2263 `fit_with_sample_weight` guard.
        if x.iter().any(|v| !v.is_finite()) {
            return Err(FerroError::InvalidParameter {
                name: "X".into(),
                reason: "Input X contains NaN or infinity.".into(),
            });
        }

        // Use a single-epoch hyper for partial_fit.
        let mut hyper = self.hyper.clone();
        hyper.max_iter = 1;

        let t = partial_fit_ova(
            x,
            y,
            &self.classes,
            &mut self.weight_matrix,
            &mut self.intercepts,
            &self.loss,
            &hyper,
            self.t,
        );
        self.t = t;

        Ok(self)
    }
}

impl<F: Float + Send + Sync + ScalarOperand + 'static> PartialFit<Array2<F>, Array1<usize>>
    for SGDClassifier<F>
{
    type FitResult = FittedSGDClassifier<F>;
    type Error = FerroError;

    /// Initial call to `partial_fit` on an unfitted classifier.
    ///
    /// Equivalent to `fit` but with a single epoch, enabling subsequent
    /// incremental calls.
    ///
    /// # Errors
    ///
    /// Same as [`Fit::fit`].
    fn partial_fit(
        self,
        x: &Array2<F>,
        y: &Array1<usize>,
    ) -> Result<FittedSGDClassifier<F>, FerroError> {
        validate_clf_params(
            x,
            y,
            &self.learning_rate,
            self.eta0,
            self.alpha,
            self.l1_ratio,
            self.validation_fraction,
        )?;

        // sklearn `SGDClassifier.partial_fit` validates X through
        // `_validate_data(force_all_finite=True)` (`_stochastic_gradient.py:596`),
        // raising `ValueError("Input X contains NaN.")` for any non-finite X
        // BEFORE the kernel (`validate_clf_params` checks shape/params, not
        // finiteness). `y` is integer labels (always finite by type).
        if x.iter().any(|v| !v.is_finite()) {
            return Err(FerroError::InvalidParameter {
                name: "X".into(),
                reason: "Input X contains NaN or infinity.".into(),
            });
        }

        let n_features = x.ncols();
        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: "SGDClassifier requires at least 2 distinct classes".into(),
            });
        }

        let mut hyper = clf_hyper(&self);
        hyper.max_iter = 1;
        let loss_enum = self.loss;

        // Initial `partial_fit` does not carry per-sample/class weights here, so
        // the expanded class weights and sample weights are uniform `1.0` — the
        // all-ones path is byte-identical to the pre-weighting kernel.
        let expanded: Vec<F> = vec![F::one(); classes.len()];
        let sw: Vec<F> = vec![F::one(); x.nrows()];

        let (weight_matrix, intercepts, t) = fit_ova(
            x, y, &classes, n_features, &loss_enum, &hyper, 0, &expanded, &sw,
        )?;

        Ok(FittedSGDClassifier {
            weight_matrix,
            intercepts,
            classes,
            n_features,
            loss: loss_enum,
            hyper: clf_hyper(&self),
            t,
        })
    }
}

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

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

// Pipeline integration.
impl<F> PipelineEstimator<F> for SGDClassifier<F>
where
    F: Float + ToPrimitive + FromPrimitive + ScalarOperand + Send + Sync + 'static,
{
    fn fit_pipeline(
        &self,
        x: &Array2<F>,
        y: &Array1<F>,
    ) -> Result<Box<dyn FittedPipelineEstimator<F>>, FerroError> {
        let y_usize: Array1<usize> = y.mapv(|v| v.to_usize().unwrap_or(0));
        let fitted = self.fit(x, &y_usize)?;
        Ok(Box::new(FittedSGDClassifierPipeline(fitted)))
    }
}

/// Wrapper for pipeline integration that converts predictions to float.
struct FittedSGDClassifierPipeline<F>(FittedSGDClassifier<F>)
where
    F: Float + Send + Sync + 'static;

// Safety: inner type fields are Send + Sync.
unsafe impl<F> Send for FittedSGDClassifierPipeline<F> where F: Float + Send + Sync + 'static {}
unsafe impl<F> Sync for FittedSGDClassifierPipeline<F> where F: Float + Send + Sync + 'static {}

impl<F> FittedPipelineEstimator<F> for FittedSGDClassifierPipeline<F>
where
    F: Float + ToPrimitive + FromPrimitive + ScalarOperand + Send + Sync + 'static,
{
    fn predict_pipeline(&self, x: &Array2<F>) -> Result<Array1<F>, FerroError> {
        let preds = self.0.predict(x)?;
        Ok(preds.mapv(|v| F::from_usize(v).unwrap_or_else(F::nan)))
    }
}

// ---------------------------------------------------------------------------
// SGDRegressor
// ---------------------------------------------------------------------------

/// Stochastic Gradient Descent regressor.
///
/// Supports several loss functions for regression, trained using stochastic
/// gradient descent with configurable learning rate schedules.
///
/// # Type Parameters
///
/// - `F`: The floating-point type (`f32` or `f64`).
///
/// # Examples
///
/// ```
/// use ferrolearn_linear::sgd::SGDRegressor;
/// use ferrolearn_core::{Fit, Predict};
/// use ndarray::{array, Array2};
///
/// 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 = SGDRegressor::<f64>::new();
/// let fitted = model.fit(&x, &y).unwrap();
/// let preds = fitted.predict(&x).unwrap();
/// ```
#[derive(Debug, Clone)]
pub struct SGDRegressor<F> {
    /// The loss function to use.
    pub loss: RegressorLoss<F>,
    /// The learning rate schedule.
    pub learning_rate: LearningRateSchedule<F>,
    /// Initial learning rate.
    pub eta0: F,
    /// Regularization strength (`alpha`).
    pub alpha: F,
    /// Regularization penalty (`l2`/`l1`/`elasticnet`). Defaults to `L2`.
    pub penalty: Penalty,
    /// Elastic-net mixing parameter; only used when `penalty == ElasticNet`.
    /// Defaults to `0.15` (sklearn default).
    pub l1_ratio: F,
    /// Maximum number of passes over the training data.
    pub max_iter: usize,
    /// Convergence tolerance.
    pub tol: F,
    /// Optional random seed for sample shuffling.
    pub random_state: Option<u64>,
    /// Power parameter for inverse scaling schedule.
    pub power_t: F,
    /// Whether to shuffle the training data after each epoch. Defaults to
    /// `true` (sklearn `SGDRegressor(shuffle=True)`,
    /// `_stochastic_gradient.py:2038`).
    pub shuffle: bool,
    /// Number of consecutive epochs with no loss improvement (beyond `tol`)
    /// before convergence triggers, or — under the `adaptive` schedule — before
    /// `eta` is divided by 5. Defaults to `5` (sklearn
    /// `_stochastic_gradient.py` `n_iter_no_change=5`, `_sgd_fast.pyx.tp:698`).
    pub n_iter_no_change: usize,
    /// Whether to fit (update) the intercept. Defaults to `true` (sklearn
    /// `SGDRegressor(fit_intercept=True)`, `_stochastic_gradient.py:2031`,
    /// constraint `["boolean"]` at `:86`). When `false` the intercept is never
    /// updated and stays at its init value `0` (`_sgd_fast.pyx.tp:639-644`: the
    /// intercept update is gated on `if fit_intercept == 1`).
    pub fit_intercept: bool,
    /// Averaged-SGD (ASGD) threshold. `0` disables averaging (the default,
    /// matching sklearn `average=False`); `1` averages from the first step
    /// (sklearn `average=True`); `N > 1` begins averaging once the global step
    /// counter `t >= N` (sklearn `average=N`). The averaged weights/intercept
    /// replace the plain ones at fit-end when `average <= self.t_ - 1`
    /// (`_sgd_fast.pyx.tp:646-654`, `_stochastic_gradient.py:834-836`).
    pub average: usize,
    /// Whether to stop training early based on a held-out validation score.
    /// Defaults to `false` (sklearn `SGDRegressor(early_stopping=False)`,
    /// `_stochastic_gradient.py:114`, constraint `["boolean"]` at `:524`). When
    /// `true`, [`validation_fraction`](Self::validation_fraction) of the training
    /// data is held out as a validation set and the epoch-end convergence rule
    /// uses the validation `R^2` instead of the training loss
    /// (`_sgd_fast.pyx.tp:678-687`).
    pub early_stopping: bool,
    /// Fraction of the training data held out as the validation set when
    /// [`early_stopping`](Self::early_stopping) is `true`. Defaults to `0.1`
    /// (sklearn `validation_fraction=0.1`, `_stochastic_gradient.py:115`). Must
    /// lie in the open interval `(0, 1)`
    /// (constraint `Interval(Real, 0, 1, closed="neither")` at `:525`).
    pub validation_fraction: F,
}

impl<F: Float> SGDRegressor<F> {
    /// Create a new `SGDRegressor` with default settings.
    ///
    /// Defaults match scikit-learn's `SGDRegressor.__init__`
    /// (`_stochastic_gradient.py:2042-2068`): `loss = SquaredError`,
    /// `learning_rate = InvScaling`, `eta0 = 0.01`, `alpha = 0.0001`,
    /// `penalty = L2`, `l1_ratio = 0.15`, `max_iter = 1000`, `tol = 1e-3`,
    /// `power_t = 0.25`.
    #[must_use]
    pub fn new() -> Self {
        Self {
            loss: RegressorLoss::SquaredError,
            n_iter_no_change: 5,
            average: 0,
            fit_intercept: true,
            shuffle: true,
            penalty: Penalty::L2,
            l1_ratio: cst(0.15),
            learning_rate: LearningRateSchedule::InvScaling,
            eta0: cst(0.01),
            alpha: cst(0.0001),
            max_iter: 1000,
            tol: cst(1e-3),
            random_state: None,
            power_t: cst(0.25),
            early_stopping: false,
            validation_fraction: cst(0.1),
        }
    }

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

    /// Set the learning rate schedule.
    #[must_use]
    pub fn with_learning_rate(mut self, lr: LearningRateSchedule<F>) -> Self {
        self.learning_rate = lr;
        self
    }

    /// Set the initial learning rate.
    #[must_use]
    pub fn with_eta0(mut self, eta0: F) -> Self {
        self.eta0 = eta0;
        self
    }

    /// Set the regularization strength.
    #[must_use]
    pub fn with_alpha(mut self, alpha: F) -> Self {
        self.alpha = alpha;
        self
    }

    /// Set the regularization penalty (`l2`/`l1`/`elasticnet`).
    #[must_use]
    pub fn with_penalty(mut self, penalty: Penalty) -> Self {
        self.penalty = penalty;
        self
    }

    /// Set the elastic-net mixing parameter (`l1_ratio`, used only when
    /// `penalty == ElasticNet`).
    #[must_use]
    pub fn with_l1_ratio(mut self, l1_ratio: F) -> Self {
        self.l1_ratio = l1_ratio;
        self
    }

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

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

    /// Set the random seed for reproducibility.
    #[must_use]
    pub fn with_random_state(mut self, seed: u64) -> Self {
        self.random_state = Some(seed);
        self
    }

    /// Set the power parameter for inverse scaling.
    #[must_use]
    pub fn with_power_t(mut self, power_t: F) -> Self {
        self.power_t = power_t;
        self
    }

    /// Set whether the training data is shuffled after each epoch.
    ///
    /// Mirrors sklearn's `shuffle` parameter (default `True`,
    /// `_stochastic_gradient.py:2038`). With `false` the samples are visited in
    /// index order `0..n-1` every epoch (`_sgd_fast.pyx.tp:579-581`), making the
    /// fit fully deterministic and cross-impl comparable to sklearn.
    #[must_use]
    pub fn with_shuffle(mut self, shuffle: bool) -> Self {
        self.shuffle = shuffle;
        self
    }

    /// Set the number of consecutive non-improving epochs before convergence
    /// (or, under the `adaptive` schedule, before `eta` is divided by 5).
    ///
    /// Mirrors sklearn's `n_iter_no_change` parameter (default `5`,
    /// `_stochastic_gradient.py`); the epoch-end stop rule at
    /// `_sgd_fast.pyx.tp:698` triggers once `no_improvement_count` reaches it.
    #[must_use]
    pub fn with_n_iter_no_change(mut self, n_iter_no_change: usize) -> Self {
        self.n_iter_no_change = n_iter_no_change;
        self
    }

    /// Set whether the intercept (bias) term is fit.
    ///
    /// Mirrors sklearn's `fit_intercept` parameter (default `True`,
    /// `_stochastic_gradient.py:2031`, constraint `["boolean"]` at `:86`). With
    /// `false` the intercept update is skipped every step
    /// (`_sgd_fast.pyx.tp:639-644`: `if fit_intercept == 1: ... intercept += ...`)
    /// and the intercept stays at its init value `0`.
    #[must_use]
    pub fn with_fit_intercept(mut self, fit_intercept: bool) -> Self {
        self.fit_intercept = fit_intercept;
        self
    }

    /// Set the averaged-SGD (ASGD) threshold.
    ///
    /// Mirrors sklearn's `average` parameter (default `False`,
    /// `_stochastic_gradient.py:2068`). `0` disables averaging (the plain SGD
    /// trajectory, byte-identical to the unaveraged kernel); `with_average(1)`
    /// is sklearn `average=True` (average from the first step); `with_average(N)`
    /// is sklearn `average=N` (begin averaging once the global step counter
    /// `t >= N`). The running mean of the post-update weights/intercept replaces
    /// the plain `coef_`/`intercept_` at fit-end when `average <= self.t_ - 1`
    /// (`_sgd_fast.pyx.tp:646-654`, `_stochastic_gradient.py:834-836`).
    #[must_use]
    pub fn with_average(mut self, average: usize) -> Self {
        self.average = average;
        self
    }

    /// Enable or disable early stopping on a held-out validation score.
    ///
    /// Mirrors sklearn's `early_stopping` parameter (default `False`,
    /// `_stochastic_gradient.py:114`, constraint `["boolean"]` at `:524`). When
    /// enabled, [`with_validation_fraction`](Self::with_validation_fraction) of
    /// the data is held out and the epoch-end convergence is driven by the
    /// validation `R^2` rather than the training loss
    /// (`_sgd_fast.pyx.tp:678-687`).
    #[must_use]
    pub fn with_early_stopping(mut self, early_stopping: bool) -> Self {
        self.early_stopping = early_stopping;
        self
    }

    /// Set the fraction of the training data held out for early-stopping
    /// validation.
    ///
    /// Mirrors sklearn's `validation_fraction` parameter (default `0.1`,
    /// `_stochastic_gradient.py:115`, constraint
    /// `Interval(Real, 0, 1, closed="neither")` at `:525`). Only used when
    /// [`early_stopping`](Self::early_stopping) is `true`; validated to the open
    /// interval `(0, 1)` at fit time.
    #[must_use]
    pub fn with_validation_fraction(mut self, validation_fraction: F) -> Self {
        self.validation_fraction = validation_fraction;
        self
    }
}

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

/// Extract hyperparameter bundle from an `SGDRegressor`.
fn reg_hyper<F: Float>(reg: &SGDRegressor<F>) -> SGDHyper<F> {
    SGDHyper {
        learning_rate: reg.learning_rate,
        eta0: reg.eta0,
        alpha: reg.alpha,
        max_iter: reg.max_iter,
        tol: reg.tol,
        random_state: reg.random_state,
        power_t: reg.power_t,
        penalty: reg.penalty,
        l1_ratio: reg.l1_ratio,
        shuffle: reg.shuffle,
        n_iter_no_change: reg.n_iter_no_change,
        fit_intercept: reg.fit_intercept,
        one_class: false,
        average: reg.average,
        early_stopping: reg.early_stopping,
        validation_fraction: reg.validation_fraction,
    }
}

/// Train a single regressor via SGD, updating `weights` and `intercept`
/// in place. Returns the final loss and step counter.
///
/// `sample_w[i]` is the per-sample weight `sample_weight[i]`, scaling ONLY the
/// gradient-derived part of each update (`update *= class_weight * sample_weight`
/// with `class_weight = 1` for regression, `_sgd_fast.pyx.tp:599-602,630`); the
/// L2 shrink and L1 truncation are unscaled. An all-ones `sample_w` (the default
/// `fit` path) reproduces the byte-identical unweighted behaviour.
#[allow(
    clippy::too_many_arguments,
    reason = "threads the per-sample weight vector"
)]
fn train_regressor_sgd<F, L>(
    x: &Array2<F>,
    y: &Array1<F>,
    weights: &mut Array1<F>,
    intercept: &mut F,
    loss_fn: &L,
    hyper: &SGDHyper<F>,
    initial_t: usize,
    sample_w: &[F],
    val_set: Option<(&Array2<F>, &Array1<F>)>,
) -> (F, usize)
where
    F: Float + ScalarOperand + Send + Sync + 'static,
    L: Loss<F>,
{
    let n_samples = x.nrows();
    let n_features = x.ncols();
    let mut t = initial_t;
    // Epoch-end convergence state, mirroring `_sgd_fast.pyx.tp:525,532-534`:
    // `best_loss = INFINITY`, `no_improvement_count = 0`. `current_eta` carries
    // the adaptive-schedule eta (`eta = eta / 5` decay, `:700`). Under early
    // stopping `best_score = -INFINITY` (higher score is better, `:533`).
    let mut best_loss = F::infinity();
    let mut best_score = F::neg_infinity();
    // Early stopping uses the validation-R^2 branch only when a validation set
    // was split off in the `Fit` path (`_stochastic_gradient.py:79`,
    // `RegressorMixin.score` -> `r2_score`).
    let early_stopping = hyper.early_stopping && val_set.is_some();
    let mut current_eta = hyper.eta0;
    let mut no_improve_count: usize = 0;
    // `tol = None` upstream becomes `-INFINITY`, disabling the stop rule
    // (`tol > -INFINITY` is false, `:690`). ferrolearn encodes that as a
    // finite/-inf `tol`; the criterion is active iff `tol > -inf`.
    let tol_active = hyper.tol > F::neg_infinity();
    let max_dloss = max_dloss::<F>();
    let mut indices: Vec<usize> = (0..n_samples).collect();
    // `optimal_init` (the `optimal` schedule's t0 offset) depends on the loss
    // and alpha, so it is computed once per fit, before the epoch loop
    // (`_sgd_fast.pyx.tp:565-570`).
    let opt_init = optimal_init(loss_fn, hyper.alpha);

    // Effective l1_ratio from the penalty (`_sgd_fast.pyx.tp:558-561`):
    // `L2 -> 0.0`, `L1 -> 1.0`, `ElasticNet -> user l1_ratio`.
    let eff = effective_l1_ratio(hyper.penalty, hyper.l1_ratio);
    let apply_l1 = matches!(hyper.penalty, Penalty::L1 | Penalty::ElasticNet);
    // Tsuruoka cumulative-penalty state (`u` scalar, `q` per-feature), allocated
    // once and persisting for the whole fit (`_sgd_fast.pyx.tp:551-556`).
    let mut u = F::zero();
    let mut q: Array1<F> = Array1::zeros(n_features);

    // Averaged-SGD (ASGD) accumulators (`_sgd_fast.pyx.tp:646-654`). Direct
    // running-mean form of sklearn's lazy `w.add_average` (mathematically
    // identical for plain arrays): once `t >= average`, `avg += (current - avg)
    // / (t - average + 1)`. A passive observer — never fed back into the live
    // `weights`/`intercept` trajectory.
    let mut average_coef: Array1<F> = Array1::zeros(n_features);
    let mut average_intercept = F::zero();

    let mut rng = match hyper.random_state {
        Some(seed) => rand::rngs::StdRng::seed_from_u64(seed),
        None => rand::rngs::StdRng::from_os_rng(),
    };

    let mut total_loss = F::zero();

    for _epoch in 0..hyper.max_iter {
        // sklearn shuffles the sample order each epoch only when `shuffle` is
        // set (`_sgd_fast.pyx.tp:579-580`: `if shuffle: dataset.shuffle(seed)`).
        // With `shuffle == false` the indices stay `0..n-1` every epoch, exactly
        // matching sklearn's no-shuffle index order (`:581 for i in range(n)`).
        if hyper.shuffle {
            indices.shuffle(&mut rng);
        }
        let mut epoch_loss = F::zero();

        for &i in &indices {
            t += 1;

            let eta = match hyper.learning_rate {
                LearningRateSchedule::Adaptive => current_eta,
                _ => compute_lr(
                    &hyper.learning_rate,
                    hyper.eta0,
                    hyper.alpha,
                    hyper.power_t,
                    opt_init,
                    t,
                ),
            };

            let xi = x.row(i);
            let mut y_pred = *intercept;
            for j in 0..n_features {
                y_pred = y_pred + weights[j] * xi[j];
            }

            // Clip the gradient to `[-MAX_DLOSS, MAX_DLOSS]` before forming the
            // update, matching `_sgd_fast.pyx.tp:613-620`.
            let grad = loss_fn
                .gradient(y[i], y_pred)
                .max(-max_dloss)
                .min(max_dloss);
            // Per-sample weight scaling: `update *= class_weight * sample_weight`
            // with `class_weight = 1` for regression (`_sgd_fast.pyx.tp:630`).
            // `g` is the sample-weighted gradient, scaling ONLY the gradient term.
            let g = grad * sample_w[i];
            // `sumloss` is the SUM (not mean) of per-sample losses over the
            // epoch (`_sgd_fast.pyx.tp:597`), computed from the UNWEIGHTED loss.
            epoch_loss = epoch_loss + loss_fn.loss(y[i], y_pred);

            // L2 shrink: clamped multiplicative factor
            // `max(0, 1 - (1-eff)*eta*alpha)` applied to the whole weight vector
            // BEFORE the gradient add (`_sgd_fast.pyx.tp:632-635`); for pure L2
            // (`eff=0`) this is `max(0, 1-eta*alpha)`, for L1 (`eff=1`) it is 1.
            let shrink = (F::one() - (F::one() - eff) * eta * hyper.alpha).max(F::zero());
            for j in 0..n_features {
                weights[j] = weights[j] * shrink;
            }
            // Gradient add `w.add(x, -eta*g)` with the scaled `g`
            // (`_sgd_fast.pyx.tp:637-638`).
            for j in 0..n_features {
                weights[j] = weights[j] - eta * g * xi[j];
            }
            // The intercept update is gated on `fit_intercept` and is NOT
            // regularized (`_sgd_fast.pyx.tp:639-644`: `if fit_intercept == 1:
            // intercept_update = update; ... intercept += intercept_update *
            // intercept_decay`). `update = -eta*g` is the SCALED update. When
            // `fit_intercept` is false the intercept is never modified and stays
            // at its init value `0` (`intercept` enters this fn as `0`).
            if hyper.fit_intercept {
                *intercept = *intercept - eta * g;
            }

            // L1 cumulative penalty (Tsuruoka truncated gradient), applied AFTER
            // the gradient add only for L1/ElasticNet (`_sgd_fast.pyx.tp:656-658`,
            // `l1penalty` at `:750-778` with `wscale = 1`).
            if apply_l1 {
                u = u + eff * eta * hyper.alpha;
                for j in 0..n_features {
                    let z = weights[j];
                    if weights[j] > F::zero() {
                        weights[j] = (weights[j] - (u + q[j])).max(F::zero());
                    } else if weights[j] < F::zero() {
                        weights[j] = (weights[j] + (u - q[j])).min(F::zero());
                    }
                    q[j] = q[j] + (weights[j] - z);
                }
            }

            // ASGD running-mean accumulation (`_sgd_fast.pyx.tp:646-654`),
            // AFTER the weight/intercept update + L1 truncation so
            // `weights`/`intercept` are the final post-step values for this
            // sample. `t` is the SAME 1-based global step the schedule used;
            // `num_iter = t - average + 1` is `>= 1` whenever `t >= average`.
            if hyper.average > 0 && t >= hyper.average {
                let num_iter = t - hyper.average + 1;
                let num_iter_f = F::from(num_iter).unwrap_or_else(F::one);
                let mu = F::one() / num_iter_f;
                for j in 0..n_features {
                    average_coef[j] = average_coef[j] + (weights[j] - average_coef[j]) * mu;
                }
                average_intercept = average_intercept + (*intercept - average_intercept) * mu;
            }
        }

        // `epoch_loss` is now the epoch `sumloss` (no mean division).
        total_loss = epoch_loss;

        // Epoch-end stop rule (`_sgd_fast.pyx.tp:678-707`). When early stopping
        // is active, score the CURRENT weights on the held-out validation set
        // (R^2, `_stochastic_gradient.py:79`) and run the score-based branch
        // (`best_score` init `-inf`, higher is better, `:678-687`); otherwise the
        // training-loss branch (`sumloss` vs `best_loss`, `:688-695`).
        let should_break = if let (true, Some((x_val, y_val))) = (early_stopping, val_set) {
            let score = r2_score(weights, *intercept, x_val, y_val);
            convergence_tail_score(
                score,
                &mut best_score,
                &mut no_improve_count,
                &mut current_eta,
                tol_active,
                hyper.tol,
                hyper.n_iter_no_change,
                matches!(hyper.learning_rate, LearningRateSchedule::Adaptive),
            )
        } else {
            convergence_tail(
                epoch_loss,
                &mut best_loss,
                &mut no_improve_count,
                &mut current_eta,
                tol_active,
                hyper.tol,
                n_samples,
                hyper.n_iter_no_change,
                matches!(hyper.learning_rate, LearningRateSchedule::Adaptive),
            )
        };
        if should_break {
            break;
        }
    }

    // ASGD finalize (`_stochastic_gradient.py:834-836`: averaged coef/intercept
    // chosen when `average <= self.t_ - 1`). `t` here equals sklearn's
    // `self.t_ - 1` on the full-fit path (`initial_t = 0`, sklearn inits
    // `self.t_ = 1` then adds `n_iter_ * n_samples`), so the condition is
    // `hyper.average <= t`.
    if hyper.average > 0 && hyper.average <= t {
        for j in 0..n_features {
            weights[j] = average_coef[j];
        }
        *intercept = average_intercept;
    }

    (total_loss, t)
}

/// Dispatch regressor training to the appropriate typed loss function.
///
/// `sample_w[i] = sample_weight[i]` is forwarded verbatim to the kernel
/// (`_sgd_fast.pyx.tp:630`, `class_weight = 1` for regression).
#[allow(
    clippy::too_many_arguments,
    reason = "threads the per-sample weight vector"
)]
fn dispatch_train_regressor<F: Float + Send + Sync + ScalarOperand + 'static>(
    x: &Array2<F>,
    y: &Array1<F>,
    w: &mut Array1<F>,
    b: &mut F,
    loss_enum: &RegressorLoss<F>,
    hyper: &SGDHyper<F>,
    initial_t: usize,
    sample_w: &[F],
    val_set: Option<(&Array2<F>, &Array1<F>)>,
) -> (F, usize) {
    match loss_enum {
        RegressorLoss::SquaredError => train_regressor_sgd(
            x,
            y,
            w,
            b,
            &SquaredError,
            hyper,
            initial_t,
            sample_w,
            val_set,
        ),
        RegressorLoss::Huber(eps) => train_regressor_sgd(
            x,
            y,
            w,
            b,
            &Huber { epsilon: *eps },
            hyper,
            initial_t,
            sample_w,
            val_set,
        ),
        RegressorLoss::EpsilonInsensitive(eps) => train_regressor_sgd(
            x,
            y,
            w,
            b,
            &EpsilonInsensitive { epsilon: *eps },
            hyper,
            initial_t,
            sample_w,
            val_set,
        ),
        RegressorLoss::SquaredEpsilonInsensitive(eps) => train_regressor_sgd(
            x,
            y,
            w,
            b,
            &SquaredEpsilonInsensitive { epsilon: *eps },
            hyper,
            initial_t,
            sample_w,
            val_set,
        ),
    }
}

/// Fitted SGD regressor.
///
/// Holds the learned weight vector and intercept. Implements [`Predict`]
/// and [`PartialFit`] to support both inference and online learning.
#[derive(Debug, Clone)]
pub struct FittedSGDRegressor<F> {
    /// Learned weight vector (one per feature).
    weights: Array1<F>,
    /// Learned intercept (bias) term.
    intercept: F,
    /// Number of features the model was trained on.
    n_features: usize,
    /// The loss function used during training.
    loss: RegressorLoss<F>,
    /// Hyperparameters for continued training.
    hyper: SGDHyper<F>,
    /// Global step counter.
    t: usize,
}

/// Validate regressor input shapes and parameters.
#[allow(
    clippy::too_many_arguments,
    reason = "threads each validated parameter"
)]
fn validate_reg_params<F: Float>(
    x: &Array2<F>,
    y: &Array1<F>,
    schedule: &LearningRateSchedule<F>,
    eta0: F,
    alpha: F,
    l1_ratio: F,
    loss: &RegressorLoss<F>,
    validation_fraction: F,
) -> Result<(), FerroError> {
    let n_samples = x.nrows();
    if n_samples != y.len() {
        return Err(FerroError::ShapeMismatch {
            expected: vec![n_samples],
            actual: vec![y.len()],
            context: "y length must match number of samples in X".into(),
        });
    }
    if n_samples == 0 {
        return Err(FerroError::InsufficientSamples {
            required: 1,
            actual: 0,
            context: "SGDRegressor requires at least one sample".into(),
        });
    }
    if schedule_requires_eta0(schedule) && eta0 <= F::zero() {
        return Err(FerroError::InvalidParameter {
            name: "eta0".into(),
            reason: "must be positive".into(),
        });
    }
    if alpha < F::zero() {
        return Err(FerroError::InvalidParameter {
            name: "alpha".into(),
            reason: "must be non-negative".into(),
        });
    }
    if l1_ratio < F::zero() || l1_ratio > F::one() {
        return Err(FerroError::InvalidParameter {
            name: "l1_ratio".into(),
            reason: "must be in the range [0, 1]".into(),
        });
    }
    // sklearn `_stochastic_gradient.py:2024`:
    // `"epsilon": [Interval(Real, 0, None, closed="left")]` — epsilon must be
    // `>= 0` (a negative epsilon raises `InvalidParameterError`). ferrolearn
    // carries epsilon inside the loss variant, so the faithful equivalent is to
    // reject a negative epsilon on the variants that have one (SquaredError has
    // none, boundary 0 is valid — closed-left).
    let epsilon = match loss {
        RegressorLoss::Huber(e)
        | RegressorLoss::EpsilonInsensitive(e)
        | RegressorLoss::SquaredEpsilonInsensitive(e) => Some(*e),
        RegressorLoss::SquaredError => None,
    };
    if let Some(e) = epsilon
        && e < F::zero()
    {
        return Err(FerroError::InvalidParameter {
            name: "epsilon".into(),
            reason: "must be in the range [0, inf)".into(),
        });
    }
    validate_validation_fraction(validation_fraction)?;
    Ok(())
}

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

    /// Fit the SGD regressor on the given data.
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if `x` and `y` have mismatched
    /// sample counts.
    /// Returns [`FerroError::InvalidParameter`] if `eta0` or `alpha` are
    /// invalid.
    fn fit(&self, x: &Array2<F>, y: &Array1<F>) -> Result<FittedSGDRegressor<F>, FerroError> {
        // Delegate to the sample-weighted path with a uniform `ones(n)` weight
        // vector — byte-identical to the unweighted kernel
        // (`_check_sample_weight` -> ones, `_stochastic_gradient.py:627`).
        let sample_weight = Array1::<F>::from_elem(x.nrows(), F::one());
        self.fit_with_sample_weight(x, y, &sample_weight)
    }
}

impl<F: Float + Send + Sync + ScalarOperand + 'static> SGDRegressor<F> {
    /// Fit the SGD regressor with explicit per-sample weights.
    ///
    /// Mirrors `SGDRegressor.fit(X, y, sample_weight=...)`. The per-sample weight
    /// scales ONLY the gradient-derived part of each update
    /// (`update *= class_weight * sample_weight` with `class_weight = 1` for
    /// regression, `_sgd_fast.pyx.tp:630`); the L2 shrink and L1 truncation are
    /// unscaled. [`Fit::fit`] delegates here with a uniform `ones(n)` vector.
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if `x`/`y`/`sample_weight` have
    /// mismatched sample counts.
    /// Returns [`FerroError::InvalidParameter`] if `eta0` or `alpha` are
    /// invalid.
    pub fn fit_with_sample_weight(
        &self,
        x: &Array2<F>,
        y: &Array1<F>,
        sample_weight: &Array1<F>,
    ) -> Result<FittedSGDRegressor<F>, FerroError> {
        validate_reg_params(
            x,
            y,
            &self.learning_rate,
            self.eta0,
            self.alpha,
            self.l1_ratio,
            &self.loss,
            self.validation_fraction,
        )?;

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

        // Non-finite input validation (#2263). sklearn `SGDRegressor.fit`
        // -> `self._validate_data(X, y, ...)` (`_stochastic_gradient.py:1476`,
        // the shared `BaseSGD` base path) keeps the default
        // `force_all_finite=True`, so `check_array` rejects any NaN or +/-inf in
        // X OR y with a `ValueError("Input X contains NaN.")` /
        // `"Input y contains NaN."` / `"... contains infinity ..."` BEFORE the
        // SGD kernel. sklearn also validates `sample_weight` via
        // `_check_sample_weight` (default `force_all_finite=True`,
        // `_stochastic_gradient.py:1501`). `y` is `Array1<F>` (float targets), so
        // X, y, AND sample_weight all need the runtime check.
        // `.iter().any(|v| !v.is_finite())` rejects both NaN and Inf (bounds-safe,
        // no panic, R-CODE-2); the finite path is byte-identical. This is the
        // SGDRegressor fit entry (`Fit::fit` delegates here with unit weights).
        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 sample_weight.iter().any(|v| !v.is_finite()) {
            return Err(FerroError::InvalidParameter {
                name: "sample_weight".into(),
                reason: "Input sample_weight contains NaN or infinity.".into(),
            });
        }

        let n_features = x.ncols();
        let hyper = reg_hyper(self);
        let mut w = Array1::<F>::zeros(n_features);
        let mut b = F::zero();
        let sw = sample_weight.to_vec();

        // Early-stopping validation split (`ShuffleSplit`,
        // `_stochastic_gradient.py:282-287`). The hold-out is split off BEFORE
        // training; the kernel trains on the train subset and scores its R^2 on
        // the held-out validation subset each epoch. sklearn does NOT refit on
        // the full data — the train-subset weights are final.
        let t = if hyper.early_stopping {
            let (train_idx, val_idx) = make_validation_split(
                n_samples,
                hyper.validation_fraction,
                hyper.random_state,
                None,
            )
            .ok_or_else(|| FerroError::InvalidParameter {
                name: "validation_fraction".into(),
                reason: "early_stopping split led to an empty train or validation set; \
                         increase the number of samples or change validation_fraction"
                    .into(),
            })?;
            let x_tr = gather_rows(x, &train_idx);
            let y_tr = gather(y, &train_idx);
            let sw_tr: Vec<F> = train_idx.iter().map(|&i| sw[i]).collect();
            let x_val = gather_rows(x, &val_idx);
            let y_val = gather(y, &val_idx);
            let (_, t) = dispatch_train_regressor(
                &x_tr,
                &y_tr,
                &mut w,
                &mut b,
                &self.loss,
                &hyper,
                0,
                &sw_tr,
                Some((&x_val, &y_val)),
            );
            t
        } else {
            let (_, t) =
                dispatch_train_regressor(x, y, &mut w, &mut b, &self.loss, &hyper, 0, &sw, None);
            t
        };

        Ok(FittedSGDRegressor {
            weights: w,
            intercept: b,
            n_features,
            loss: self.loss,
            hyper,
            t,
        })
    }
}

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

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

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

impl<F: Float + Send + Sync + ScalarOperand + 'static> PartialFit<Array2<F>, Array1<F>>
    for FittedSGDRegressor<F>
{
    type FitResult = FittedSGDRegressor<F>;
    type Error = FerroError;

    /// Incrementally train the regressor on a new batch of data.
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if `x` and `y` have mismatched
    /// sizes or `x` has the wrong number of features.
    fn partial_fit(
        mut self,
        x: &Array2<F>,
        y: &Array1<F>,
    ) -> Result<FittedSGDRegressor<F>, FerroError> {
        let n_samples = x.nrows();
        if n_samples != y.len() {
            return Err(FerroError::ShapeMismatch {
                expected: vec![n_samples],
                actual: vec![y.len()],
                context: "y length must match number of samples in X".into(),
            });
        }
        if x.ncols() != self.n_features {
            return Err(FerroError::ShapeMismatch {
                expected: vec![self.n_features],
                actual: vec![x.ncols()],
                context: "number of features must match fitted model".into(),
            });
        }

        // sklearn `SGDRegressor.partial_fit` validates X and the float target y
        // through `_validate_data(force_all_finite=True)`
        // (`_stochastic_gradient.py:1476`), raising `ValueError("Input X contains
        // NaN.")` / `"Input y contains NaN."` BEFORE the kernel. Mirrors the
        // #2263 `fit_with_sample_weight` guard.
        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(),
            });
        }

        let mut hyper = self.hyper.clone();
        hyper.max_iter = 1;
        // `partial_fit` carries no per-sample weight here; uniform `1.0`.
        let sample_w: Vec<F> = vec![F::one(); x.nrows()];

        let (_, t) = dispatch_train_regressor(
            x,
            y,
            &mut self.weights,
            &mut self.intercept,
            &self.loss,
            &hyper,
            self.t,
            &sample_w,
            None,
        );
        self.t = t;

        Ok(self)
    }
}

impl<F: Float + Send + Sync + ScalarOperand + 'static> PartialFit<Array2<F>, Array1<F>>
    for SGDRegressor<F>
{
    type FitResult = FittedSGDRegressor<F>;
    type Error = FerroError;

    /// Initial call to `partial_fit` on an unfitted regressor.
    ///
    /// Equivalent to `fit` but with a single epoch.
    ///
    /// # Errors
    ///
    /// Same as [`Fit::fit`].
    fn partial_fit(
        self,
        x: &Array2<F>,
        y: &Array1<F>,
    ) -> Result<FittedSGDRegressor<F>, FerroError> {
        validate_reg_params(
            x,
            y,
            &self.learning_rate,
            self.eta0,
            self.alpha,
            self.l1_ratio,
            &self.loss,
            self.validation_fraction,
        )?;

        // sklearn `SGDRegressor.partial_fit` validates X and the float target y
        // through `_validate_data(force_all_finite=True)`
        // (`_stochastic_gradient.py:1476`), raising `ValueError("Input X contains
        // NaN.")` / `"Input y contains NaN."` BEFORE the kernel
        // (`validate_reg_params` checks shape/params, not finiteness).
        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(),
            });
        }

        let n_features = x.ncols();
        let mut hyper = reg_hyper(&self);
        hyper.max_iter = 1;
        let mut w = Array1::<F>::zeros(n_features);
        let mut b = F::zero();
        // Initial `partial_fit` carries no per-sample weight here; uniform `1.0`.
        let sample_w: Vec<F> = vec![F::one(); x.nrows()];

        let (_, t) =
            dispatch_train_regressor(x, y, &mut w, &mut b, &self.loss, &hyper, 0, &sample_w, None);

        Ok(FittedSGDRegressor {
            weights: w,
            intercept: b,
            n_features,
            loss: self.loss,
            hyper: reg_hyper(&self),
            t,
        })
    }
}

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

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

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

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

// ---------------------------------------------------------------------------
// SGDOneClassSVM
// ---------------------------------------------------------------------------

/// Linear One-Class SVM trained by Stochastic Gradient Descent.
///
/// Mirrors scikit-learn's
/// [`SGDOneClassSVM`](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.SGDOneClassSVM.html)
/// (`_stochastic_gradient.py:2084-2668`). It solves the linear One-Class SVM
/// primal via the same SGD kernel as [`SGDClassifier`], with the targets fixed
/// to `y = ones(n)`, the Hinge loss (`threshold = 1`), the L2 penalty, and
/// `alpha = nu / 2` (`_stochastic_gradient.py:2479,2588`). The SGD intercept
/// `b` relates to the One-Class offset `rho` by `offset_ = 1 - b`
/// (`_stochastic_gradient.py:2325,2377`), and the per-sample intercept update
/// gains an extra `- 2*eta*alpha` term (`_sgd_fast.pyx.tp:641-642`).
///
/// The decision function is `decision_function(X) = X · coef_ - offset_`
/// (`_stochastic_gradient.py:2622`); `predict` returns `+1` (inlier) where the
/// decision is `>= 0` and `-1` (outlier) otherwise (`:2655-2657`).
///
/// # Type Parameters
///
/// - `F`: The floating-point type (`f32` or `f64`).
///
/// # Examples
///
/// ```
/// use ferrolearn_linear::sgd::SGDOneClassSVM;
/// use ferrolearn_core::{Fit, Predict};
/// use ndarray::{array, Array2};
///
/// let x = Array2::from_shape_vec((4, 2), vec![
///     -1.0, -1.0, -2.0, -1.0, 1.0, 1.0, 2.0, 1.0,
/// ]).unwrap();
///
/// let model = SGDOneClassSVM::<f64>::new()
///     .with_learning_rate(ferrolearn_linear::sgd::LearningRateSchedule::Constant)
///     .with_eta0(0.05)
///     .with_max_iter(10)
///     .with_shuffle(false);
/// let fitted = model.fit(&x, &()).unwrap();
/// let preds = fitted.predict(&x).unwrap();
/// assert_eq!(preds.len(), 4);
/// ```
#[derive(Debug, Clone)]
pub struct SGDOneClassSVM<F> {
    /// The `nu` parameter — an upper bound on the fraction of training errors
    /// and a lower bound on the fraction of support vectors. Must be in
    /// `(0, 1]`. Defaults to `0.5` (`_stochastic_gradient.py:2098-2102,2247`).
    pub nu: F,
    /// Whether to fit (update) the intercept. Defaults to `true`
    /// (`_stochastic_gradient.py:2104-2105,2248`).
    pub fit_intercept: bool,
    /// Maximum number of passes over the training data. Defaults to `1000`
    /// (`_stochastic_gradient.py:2107,2249`).
    pub max_iter: usize,
    /// Convergence tolerance. Defaults to `1e-3`
    /// (`_stochastic_gradient.py:2113,2250`). Set to `F::neg_infinity()` to
    /// disable the early-stop rule (the analog of sklearn's `tol=None`,
    /// `_stochastic_gradient.py:2310`).
    pub tol: F,
    /// Whether to shuffle the training data after each epoch. Defaults to
    /// `true` (`_stochastic_gradient.py:2118,2251`).
    pub shuffle: bool,
    /// The learning rate schedule. Defaults to `Optimal`
    /// (`_stochastic_gradient.py:2132,2254`).
    pub learning_rate: LearningRateSchedule<F>,
    /// Initial learning rate for the `constant`/`invscaling`/`adaptive`
    /// schedules. Defaults to `0.0` (`_stochastic_gradient.py:2145,2255`).
    pub eta0: F,
    /// Power parameter for the inverse-scaling schedule. Defaults to `0.5`
    /// (`_stochastic_gradient.py:2151,2256`).
    pub power_t: F,
    /// Optional random seed for sample shuffling
    /// (`_stochastic_gradient.py:2125`).
    pub random_state: Option<u64>,
    /// Number of consecutive non-improving epochs before convergence (or, under
    /// the `adaptive` schedule, before `eta` is divided by 5). Defaults to `5`
    /// (`_stochastic_gradient.py:2278`).
    pub n_iter_no_change: usize,
}

impl<F: Float> SGDOneClassSVM<F> {
    /// Create a new `SGDOneClassSVM` with default settings.
    ///
    /// Defaults match scikit-learn's `SGDOneClassSVM.__init__`
    /// (`_stochastic_gradient.py:2245-2281`): `nu = 0.5`,
    /// `fit_intercept = true`, `max_iter = 1000`, `tol = 1e-3`,
    /// `shuffle = true`, `learning_rate = Optimal`, `eta0 = 0.0`,
    /// `power_t = 0.5`, `n_iter_no_change = 5`.
    #[must_use]
    pub fn new() -> Self {
        Self {
            nu: cst(0.5),
            fit_intercept: true,
            max_iter: 1000,
            tol: cst(1e-3),
            shuffle: true,
            learning_rate: LearningRateSchedule::Optimal,
            eta0: cst(0.0),
            power_t: cst(0.5),
            random_state: None,
            n_iter_no_change: 5,
        }
    }

    /// Set the `nu` parameter (upper bound on the fraction of training errors).
    #[must_use]
    pub fn with_nu(mut self, nu: F) -> Self {
        self.nu = nu;
        self
    }

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

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

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

    /// Set whether the training data is shuffled after each epoch.
    #[must_use]
    pub fn with_shuffle(mut self, shuffle: bool) -> Self {
        self.shuffle = shuffle;
        self
    }

    /// Set the learning rate schedule.
    #[must_use]
    pub fn with_learning_rate(mut self, lr: LearningRateSchedule<F>) -> Self {
        self.learning_rate = lr;
        self
    }

    /// Set the initial learning rate.
    #[must_use]
    pub fn with_eta0(mut self, eta0: F) -> Self {
        self.eta0 = eta0;
        self
    }

    /// Set the power parameter for inverse scaling.
    #[must_use]
    pub fn with_power_t(mut self, power_t: F) -> Self {
        self.power_t = power_t;
        self
    }

    /// Set the random seed for reproducibility.
    #[must_use]
    pub fn with_random_state(mut self, seed: u64) -> Self {
        self.random_state = Some(seed);
        self
    }

    /// Set the number of consecutive non-improving epochs before convergence.
    #[must_use]
    pub fn with_n_iter_no_change(mut self, n_iter_no_change: usize) -> Self {
        self.n_iter_no_change = n_iter_no_change;
        self
    }

    /// Fit the linear One-Class SVM on `x` (the X-only fit shape).
    ///
    /// This is the inherent entry point mirroring sklearn's `fit(X)`
    /// (`_stochastic_gradient.py:2554-2600`). The [`Fit`] trait impl with a
    /// unit target `()` delegates here.
    ///
    /// # Errors
    ///
    /// - [`FerroError::InvalidParameter`] if `nu` is not in `(0, 1]`
    ///   (`_stochastic_gradient.py:2236`,
    ///   `Interval(Real, 0.0, 1.0, closed="right")`).
    /// - [`FerroError::InvalidParameter`] if `eta0` is not positive for the
    ///   `constant`/`invscaling`/`adaptive` schedules.
    /// - [`FerroError::InsufficientSamples`] if `x` has no rows.
    pub fn fit_one_class(&self, x: &Array2<F>) -> Result<FittedSGDOneClassSVM<F>, FerroError>
    where
        F: Send + Sync + ScalarOperand + 'static,
    {
        // `nu` constraint: `Interval(Real, 0.0, 1.0, closed="right")`, i.e.
        // `0 < nu <= 1` (`_stochastic_gradient.py:2236`).
        if self.nu <= F::zero() || self.nu > F::one() {
            return Err(FerroError::InvalidParameter {
                name: "nu".into(),
                reason: "must be in the range (0, 1]".into(),
            });
        }
        let n_samples = x.nrows();
        if n_samples == 0 {
            return Err(FerroError::InsufficientSamples {
                required: 1,
                actual: 0,
                context: "SGDOneClassSVM requires at least one sample".into(),
            });
        }

        // Non-finite input validation (#2263) — SEPARATE SGD arm. sklearn
        // `SGDOneClassSVM.fit` -> `self._validate_data(X, None, ...)`
        // (`_stochastic_gradient.py:2392`) keeps the default
        // `force_all_finite=True`, so `check_array` rejects any NaN or +/-inf in
        // X with a `ValueError("Input X contains NaN.")` / `"... contains
        // infinity ..."` BEFORE the SGD kernel. This is an X-only fit (no `y`, no
        // `sample_weight` argument), so X is the only runtime check.
        // `.iter().any(|v| !v.is_finite())` rejects both NaN and Inf (bounds-safe,
        // no panic, R-CODE-2); the finite path is byte-identical.
        if x.iter().any(|v| !v.is_finite()) {
            return Err(FerroError::InvalidParameter {
                name: "X".into(),
                reason: "Input X contains NaN or infinity.".into(),
            });
        }

        // `eta0 > 0` is required for the constant/invscaling/adaptive schedules
        // (mirrors `_more_validate_params`, `_stochastic_gradient.py:149-153`);
        // the `optimal` schedule accepts `eta0 == 0`.
        if schedule_requires_eta0(&self.learning_rate) && self.eta0 <= F::zero() {
            return Err(FerroError::InvalidParameter {
                name: "eta0".into(),
                reason: "must be positive".into(),
            });
        }

        let n_features = x.ncols();
        // sklearn: `alpha = self.nu / 2` (`_stochastic_gradient.py:2588`),
        // `penalty="l2"`, `l1_ratio=0`, `loss="hinge"` (`:2262-2265`).
        let two = cst::<F>(2.0);
        let alpha = self.nu / two;
        let hyper = SGDHyper {
            learning_rate: self.learning_rate,
            eta0: self.eta0,
            alpha,
            max_iter: self.max_iter,
            tol: self.tol,
            random_state: self.random_state,
            power_t: self.power_t,
            penalty: Penalty::L2,
            l1_ratio: F::zero(),
            shuffle: self.shuffle,
            n_iter_no_change: self.n_iter_no_change,
            fit_intercept: self.fit_intercept,
            one_class: true,
            // sklearn's `SGDOneClassSVM` has no `average` parameter — averaging is
            // always off on the one-class path (`_stochastic_gradient.py:2245-2281`).
            average: 0,
            // The one-class SVM exposes no `early_stopping`/`validation_fraction`
            // (`_stochastic_gradient.py:2245-2281`); the early-stop score branch
            // is always off, leaving the one-class trajectory byte-identical.
            early_stopping: false,
            validation_fraction: cst(0.1),
        };

        // `y = np.ones(n_samples)` (`_stochastic_gradient.py:2289`).
        let y_ones: Array1<F> = Array1::from_elem(n_samples, F::one());
        let mut w = Array1::<F>::zeros(n_features);
        // The One-Class offset is initialized to 0, so the SGD intercept starts
        // at `b = 1 - offset_ = 1` (`_stochastic_gradient.py:2238,2325`). This
        // differs from the classifier/regressor paths, which start at `b = 0`.
        let mut b = F::one();
        // One-Class SVM fit has no per-sample weighting here; the per-sample
        // weight is uniform `1.0` (byte-identical to the pre-weighting kernel).
        let sample_w: Vec<F> = vec![F::one(); n_samples];

        let (_, _t) = train_binary_sgd(
            x, &y_ones, &mut w, &mut b, &Hinge, &hyper, 0, &sample_w, None,
        );

        // `offset_ = 1 - intercept` (`_stochastic_gradient.py:2377`).
        let offset = F::one() - b;

        Ok(FittedSGDOneClassSVM {
            coef: w,
            offset,
            n_features,
        })
    }
}

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

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

    /// Fit the linear One-Class SVM. The target `y` is ignored (present for API
    /// consistency, mirroring sklearn's `fit(X, y=None)`,
    /// `_stochastic_gradient.py:2554`); the fit uses `y = ones(n)` internally.
    ///
    /// # Errors
    ///
    /// See [`SGDOneClassSVM::fit_one_class`].
    fn fit(&self, x: &Array2<F>, _y: &()) -> Result<FittedSGDOneClassSVM<F>, FerroError> {
        self.fit_one_class(x)
    }
}

/// Fitted linear One-Class SVM.
///
/// Holds the learned weight vector `coef_` and the One-Class offset `offset_`
/// (`_stochastic_gradient.py:2177-2182`). Implements [`Predict`] (returning
/// `+1`/`-1` inlier/outlier labels) and exposes [`decision_function`],
/// [`score_samples`], [`coef`], and [`offset`].
///
/// [`decision_function`]: FittedSGDOneClassSVM::decision_function
/// [`score_samples`]: FittedSGDOneClassSVM::score_samples
/// [`coef`]: FittedSGDOneClassSVM::coef
/// [`offset`]: FittedSGDOneClassSVM::offset
#[derive(Debug, Clone)]
pub struct FittedSGDOneClassSVM<F> {
    /// Weight vector (`coef_`, shape `(n_features,)`).
    coef: Array1<F>,
    /// The One-Class offset (`offset_`), a scalar: `offset_ = 1 - intercept`.
    offset: F,
    /// Number of features the model was trained on.
    n_features: usize,
}

impl<F: Float + Send + Sync + ScalarOperand + 'static> FittedSGDOneClassSVM<F> {
    /// The learned weight vector (`coef_`).
    #[must_use]
    pub fn coef(&self) -> &Array1<F> {
        &self.coef
    }

    /// The One-Class offset (`offset_`).
    ///
    /// Satisfies `decision_function = score_samples - offset_`
    /// (`_stochastic_gradient.py:2182`).
    #[must_use]
    pub fn offset(&self) -> F {
        self.offset
    }

    /// Signed distance to the separating hyperplane:
    /// `decision_function(X) = X · coef_ - offset_`
    /// (`_stochastic_gradient.py:2622`).
    ///
    /// Positive for an inlier, negative for an outlier.
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if the number of features does
    /// not match the fitted model.
    pub fn decision_function(&self, x: &Array2<F>) -> Result<Array1<F>, FerroError> {
        let n_features = x.ncols();
        if n_features != self.n_features {
            return Err(FerroError::ShapeMismatch {
                expected: vec![self.n_features],
                actual: vec![n_features],
                context: "number of features must match fitted model".into(),
            });
        }
        Ok(x.dot(&self.coef) - self.offset)
    }

    /// Raw scoring function of the samples:
    /// `score_samples(X) = decision_function(X) + offset_ = X · coef_`
    /// (`_stochastic_gradient.py:2639`).
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if the number of features does
    /// not match the fitted model.
    pub fn score_samples(&self, x: &Array2<F>) -> Result<Array1<F>, FerroError> {
        Ok(self.decision_function(x)? + self.offset)
    }
}

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

    /// Return labels (`+1` inlier, `-1` outlier) for the given feature matrix.
    ///
    /// Mirrors `_stochastic_gradient.py:2655-2657`:
    /// `y = (decision_function(X) >= 0); y[y == 0] = -1`, i.e. `+1` where the
    /// decision is `>= 0` and `-1` otherwise.
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if the number of features does
    /// not match the fitted model.
    fn predict(&self, x: &Array2<F>) -> Result<Array1<isize>, FerroError> {
        let decisions = self.decision_function(x)?;
        Ok(decisions.mapv(|d| if d >= F::zero() { 1 } else { -1 }))
    }
}

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

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

    // -----------------------------------------------------------------------
    // Early-stopping validation-score helpers (REQ-13)
    //
    // These pin the DETERMINISTIC per-epoch validation score against the live
    // sklearn 1.5.2 oracle (`sklearn.metrics.r2_score` / `accuracy_score`, the
    // regressor/classifier `_ValidationScoreCallback`,
    // `_stochastic_gradient.py:79`). The validation-set SELECTION is numpy-RNG
    // coupled and NOT verified here (same barrier as `shuffle`); the SCORE math
    // on a GIVEN (weights, intercept, val set) is fully deterministic.
    // -----------------------------------------------------------------------

    #[test]
    fn test_validation_r2_matches_sklearn() {
        // Oracle:
        //   python3 -c "import numpy as np; from sklearn.metrics import r2_score; \
        //     X=np.array([[1.,2.],[3.,1.],[0.,4.]]); w=np.array([1.,-2.]); b=0.5; \
        //     y=np.array([-2.,4.,-7.]); print(r2_score(y, X@w+b))"
        //   -> 0.8887362637362637
        let weights = array![1.0_f64, -2.0];
        let x_val = array![[1.0_f64, 2.0], [3.0, 1.0], [0.0, 4.0]];
        let y_val = array![-2.0_f64, 4.0, -7.0];
        let got = r2_score(&weights, 0.5, &x_val, &y_val);
        assert!(
            (got - 0.8887362637362637).abs() < 1e-12,
            "r2 {got} != sklearn 0.8887362637362637"
        );
    }

    #[test]
    fn test_validation_r2_constant_y_edge_cases() {
        // sklearn `r2_score` with SS_tot == 0 (constant y_val): perfect const
        // prediction -> 1.0, imperfect -> 0.0.
        //   python3 -c "from sklearn.metrics import r2_score; \
        //     print(r2_score([5,5,5],[5,5,5]), r2_score([5,5,5],[4,4,4]))"
        //   -> 1.0 0.0
        let weights = array![0.0_f64, 0.0];
        let x = array![[1.0_f64, 2.0], [3.0, 1.0], [0.0, 4.0]];
        let y_const = array![5.0_f64, 5.0, 5.0];
        // weights=0, intercept=5 -> all predictions 5 -> perfect -> 1.0.
        let got_perfect = r2_score(&weights, 5.0, &x, &y_const);
        assert!(
            (got_perfect - 1.0).abs() < 1e-12,
            "perfect-const r2 {got_perfect} != 1.0"
        );
        // weights=0, intercept=4 -> all predictions 4, y all 5 -> imperfect -> 0.0.
        let got_imperfect = r2_score(&weights, 4.0, &x, &y_const);
        assert!(
            got_imperfect.abs() < 1e-12,
            "imperfect-const r2 {got_imperfect} != 0.0"
        );
    }

    #[test]
    fn test_validation_binary_accuracy_matches_sklearn() {
        // Oracle:
        //   python3 -c "import numpy as np; from sklearn.metrics import accuracy_score; \
        //     X=np.array([[1.,1.],[2.,2.],[0.,0.],[3.,1.]]); w=np.array([1.,1.]); b=-3.; \
        //     dec=X@w+b; pred=np.where(dec>=0,1.,-1.); y=np.array([1.,1.,-1.,1.]); \
        //     print(accuracy_score(y,pred))"
        //   -> 0.75  (dec=[-1,1,-3,1] -> pred=[-1,1,-1,1], y=[1,1,-1,1] -> 3/4)
        let weights = array![1.0_f64, 1.0];
        let x_val = array![[1.0_f64, 1.0], [2.0, 2.0], [0.0, 0.0], [3.0, 1.0]];
        let y_val = array![1.0_f64, 1.0, -1.0, 1.0];
        let got = binary_accuracy(&weights, -3.0, &x_val, &y_val);
        assert!((got - 0.75).abs() < 1e-12, "accuracy {got} != sklearn 0.75");
    }

    // -----------------------------------------------------------------------
    // Loss function tests
    // -----------------------------------------------------------------------

    #[test]
    fn test_hinge_loss_correct_side() {
        let h = Hinge;
        // y=1, pred=2 => margin=2 >= 1 => loss=0
        assert!((Loss::<f64>::loss(&h, 1.0, 2.0) - 0.0).abs() < 1e-10);
        assert!((Loss::<f64>::gradient(&h, 1.0, 2.0) - 0.0).abs() < 1e-10);
    }

    #[test]
    fn test_hinge_loss_wrong_side() {
        let h = Hinge;
        // y=1, pred=-0.5 => margin=-0.5 < 1 => loss=1.5
        assert!((Loss::<f64>::loss(&h, 1.0, -0.5) - 1.5).abs() < 1e-10);
        assert!((Loss::<f64>::gradient(&h, 1.0, -0.5) - (-1.0)).abs() < 1e-10);
    }

    #[test]
    fn test_log_loss_zero_pred() {
        let l = LogLoss;
        // y=1, pred=0 => loss=log(1+exp(0))=log(2)
        let loss = Loss::<f64>::loss(&l, 1.0, 0.0);
        assert!((loss - 2.0_f64.ln()).abs() < 1e-10);
    }

    #[test]
    fn test_log_loss_large_correct() {
        let l = LogLoss;
        // y=1, pred=20 => very small loss
        let loss = Loss::<f64>::loss(&l, 1.0, 20.0);
        assert!(loss < 1e-5);
    }

    #[test]
    fn test_squared_error_loss() {
        let s = SquaredError;
        assert!((Loss::<f64>::loss(&s, 3.0, 1.0) - 2.0).abs() < 1e-10);
        assert!((Loss::<f64>::gradient(&s, 3.0, 1.0) - (-2.0)).abs() < 1e-10);
    }

    #[test]
    fn test_modified_huber_loss() {
        let mh = ModifiedHuber;
        // y=1, pred=2 => z=2 >= 1 => loss=0
        assert!((Loss::<f64>::loss(&mh, 1.0, 2.0)).abs() < 1e-10);
        // y=1, pred=0.5 => z=0.5 => loss=(1-0.5)^2=0.25
        assert!((Loss::<f64>::loss(&mh, 1.0, 0.5) - 0.25).abs() < 1e-10);
        // y=1, pred=-2 => z=-2 < -1 => loss=-4*(-2)=8
        assert!((Loss::<f64>::loss(&mh, 1.0, -2.0) - 8.0).abs() < 1e-10);
    }

    #[test]
    fn test_huber_loss_quadratic_region() {
        let h = Huber { epsilon: 1.0_f64 };
        // |y - p| = 0.5 <= 1.0 => quadratic
        assert!((Loss::<f64>::loss(&h, 1.0, 0.5) - 0.125).abs() < 1e-10);
    }

    #[test]
    fn test_huber_loss_linear_region() {
        let h = Huber { epsilon: 1.0_f64 };
        // |y - p| = 3 > 1 => linear: 1*(3 - 0.5) = 2.5
        assert!((Loss::<f64>::loss(&h, 3.0, 0.0) - 2.5).abs() < 1e-10);
    }

    #[test]
    fn test_epsilon_insensitive_inside() {
        let ei = EpsilonInsensitive { epsilon: 0.1_f64 };
        // |y - p| = 0.05 <= 0.1 => loss=0
        assert!((Loss::<f64>::loss(&ei, 1.0, 0.95)).abs() < 1e-10);
    }

    #[test]
    fn test_epsilon_insensitive_outside() {
        let ei = EpsilonInsensitive { epsilon: 0.1_f64 };
        // |y - p| = 0.5 > 0.1 => loss=0.4
        assert!((Loss::<f64>::loss(&ei, 1.0, 0.5) - 0.4).abs() < 1e-10);
    }

    // -----------------------------------------------------------------------
    // Learning rate tests
    // -----------------------------------------------------------------------

    #[test]
    fn test_constant_lr() {
        let lr: LearningRateSchedule<f64> = LearningRateSchedule::Constant;
        // optimal_init is ignored by the Constant schedule.
        assert!((compute_lr(&lr, 0.1, 0.01, 0.25, 1.0, 100) - 0.1).abs() < 1e-10);
    }

    #[test]
    fn test_optimal_lr() {
        // sklearn `_sgd_fast.pyx.tp:592`: `eta = 1/(alpha*(optimal_init+t-1))`.
        // At t=1 this equals `initial_eta0 = 1/(alpha*optimal_init)`. With
        // alpha=0.01 and Hinge `dloss(1,-typw)` (= -1, |·|=1), the live oracle
        // gives `typw = sqrt(1/sqrt(0.01)) = 3.1622776601683795` and
        // `optimal_init = 1/(typw*0.01) = 31.62277660168379` (computed by
        // `_sgd_fast`'s init block), so eta@t=1 == typw.
        let lr: LearningRateSchedule<f64> = LearningRateSchedule::Optimal;
        const OPTIMAL_INIT: f64 = 31.62277660168379;
        const ETA_T1: f64 = 3.1622776601683795; // = typw = initial_eta0
        assert!((compute_lr(&lr, 0.0, 0.01, 0.5, OPTIMAL_INIT, 1) - ETA_T1).abs() < 1e-9);
        // At t=10: eta = 1/(0.01*(31.62277660168379 + 9)).
        let expected_t10 = 1.0 / (0.01 * (OPTIMAL_INIT + 9.0));
        assert!((compute_lr(&lr, 0.0, 0.01, 0.5, OPTIMAL_INIT, 10) - expected_t10).abs() < 1e-10);
    }

    #[test]
    fn test_optimal_init_matches_oracle() {
        // optimal_init derivation matches the live sklearn oracle (Hinge):
        // python3 -c "import numpy as np; from sklearn.linear_model._sgd_fast \
        //   import Hinge; a=0.01; typw=np.sqrt(1/np.sqrt(a)); \
        //   e0=typw/max(1.0,abs(Hinge(1.0).py_dloss(1.0,-typw))); print(1/(e0*a))"
        // -> 31.62277660168379
        const SK_OPTIMAL_INIT: f64 = 31.62277660168379;
        let got = optimal_init(&Hinge, 0.01_f64);
        assert!((got - SK_OPTIMAL_INIT).abs() < 1e-9, "got {got}");
        // alpha == 0 returns 1.0 (schedule guarded / unused upstream).
        assert!((optimal_init(&Hinge, 0.0_f64) - 1.0).abs() < 1e-12);
    }

    #[test]
    fn test_invscaling_lr() {
        let lr: LearningRateSchedule<f64> = LearningRateSchedule::InvScaling;
        // eta = 0.1 / 10^0.5 = 0.1 / 3.162... ~= 0.0316...
        let result = compute_lr(&lr, 0.1, 0.01, 0.5, 1.0, 10);
        let expected = 0.1 / 10.0_f64.sqrt();
        assert!((result - expected).abs() < 1e-10);
    }

    #[test]
    fn test_adaptive_lr_returns_eta0() {
        let lr: LearningRateSchedule<f64> = LearningRateSchedule::Adaptive;
        assert!((compute_lr(&lr, 0.05, 0.01, 0.25, 1.0, 100) - 0.05).abs() < 1e-10);
    }

    // -----------------------------------------------------------------------
    // SGDClassifier tests
    // -----------------------------------------------------------------------

    #[test]
    fn test_sgd_classifier_binary() {
        // Well-separated clusters centered near origin for SGD stability.
        let x = Array2::from_shape_vec(
            (8, 2),
            vec![
                -2.0, -2.0, -1.5, -2.0, -2.0, -1.5, -1.5, -1.5, 2.0, 2.0, 1.5, 2.0, 2.0, 1.5, 1.5,
                1.5,
            ],
        )
        .unwrap();
        let y = array![0, 0, 0, 0, 1, 1, 1, 1];

        let clf = SGDClassifier::<f64>::new()
            .with_loss(ClassifierLoss::Log)
            .with_random_state(42)
            .with_max_iter(1000)
            .with_eta0(0.01);
        let fitted = clf.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 >= 6 correct, got {correct}");
    }

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

        let clf = SGDClassifier::<f64>::new()
            .with_loss(ClassifierLoss::Log)
            .with_random_state(42)
            .with_max_iter(500);
        let fitted = clf.fit(&x, &y).unwrap();
        let preds = fitted.predict(&x).unwrap();

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

    #[test]
    fn test_sgd_classifier_multiclass() {
        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,
            ],
        )
        .unwrap();
        let y = array![0, 0, 0, 1, 1, 1, 2, 2, 2];

        let clf = SGDClassifier::<f64>::new()
            .with_random_state(42)
            .with_max_iter(1000)
            .with_eta0(0.01);
        let fitted = clf.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 >= 6 correct for multiclass, got {correct}"
        );
    }

    #[test]
    fn test_sgd_classifier_shape_mismatch_fit() {
        let x = Array2::from_shape_vec((3, 1), vec![1.0, 2.0, 3.0]).unwrap();
        let y = array![0, 1]; // Wrong length
        let clf = SGDClassifier::<f64>::new();
        assert!(clf.fit(&x, &y).is_err());
    }

    #[test]
    fn test_sgd_classifier_shape_mismatch_predict() {
        let x =
            Array2::from_shape_vec((4, 2), vec![1.0, 1.0, 2.0, 2.0, 8.0, 8.0, 9.0, 9.0]).unwrap();
        let y = array![0, 0, 1, 1];
        let clf = SGDClassifier::<f64>::new().with_random_state(42);
        let fitted = clf.fit(&x, &y).unwrap();

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

    #[test]
    fn test_sgd_classifier_single_class_error() {
        let x = Array2::from_shape_vec((3, 1), vec![1.0, 2.0, 3.0]).unwrap();
        let y = array![0, 0, 0];
        let clf = SGDClassifier::<f64>::new();
        assert!(clf.fit(&x, &y).is_err());
    }

    #[test]
    fn test_sgd_classifier_invalid_eta0() {
        let x = Array2::from_shape_vec((4, 1), vec![1.0, 2.0, 3.0, 4.0]).unwrap();
        let y = array![0, 0, 1, 1];
        // sklearn enforces `eta0 > 0` only for constant/invscaling/adaptive
        // (`_stochastic_gradient.py:149-153`); under the default `optimal`
        // schedule `eta0 = 0.0` is valid, so use `constant` to hit the reject.
        let clf = SGDClassifier::<f64>::new()
            .with_learning_rate(LearningRateSchedule::Constant)
            .with_eta0(0.0);
        assert!(clf.fit(&x, &y).is_err());
    }

    #[test]
    fn test_sgd_classifier_optimal_eta0_zero_ok() {
        // Default `optimal` schedule with `eta0 = 0.0` (sklearn default) must
        // NOT be rejected by validation (`_stochastic_gradient.py:149-153`).
        let x = Array2::from_shape_vec((4, 1), vec![-2.0, -1.0, 1.0, 2.0]).unwrap_or_default();
        let y = array![0, 0, 1, 1];
        let clf = SGDClassifier::<f64>::new().with_random_state(42);
        assert!((clf.eta0 - 0.0).abs() < 1e-12);
        assert!(clf.fit(&x, &y).is_ok());
    }

    #[test]
    fn test_sgd_classifier_invalid_alpha() {
        let x = Array2::from_shape_vec((4, 1), vec![1.0, 2.0, 3.0, 4.0]).unwrap();
        let y = array![0, 0, 1, 1];
        let clf = SGDClassifier::<f64>::new().with_alpha(-1.0);
        assert!(clf.fit(&x, &y).is_err());
    }

    #[test]
    fn test_sgd_classifier_has_coefficients() {
        let x =
            Array2::from_shape_vec((4, 2), vec![1.0, 1.0, 2.0, 2.0, 8.0, 8.0, 9.0, 9.0]).unwrap();
        let y = array![0, 0, 1, 1];
        let clf = SGDClassifier::<f64>::new().with_random_state(42);
        let fitted = clf.fit(&x, &y).unwrap();
        assert_eq!(fitted.coefficients().len(), 2);
    }

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

        let clf = SGDClassifier::<f64>::new().with_random_state(42);
        let fitted = clf.partial_fit(&x, &y).unwrap();
        let fitted = fitted.partial_fit(&x, &y).unwrap();
        let preds = fitted.predict(&x).unwrap();
        assert_eq!(preds.len(), 4);
    }

    #[test]
    fn test_sgd_classifier_partial_fit_chain() {
        // Test the chaining pattern:
        // model.partial_fit(&b1, &y1)?.partial_fit(&b2, &y2)?.predict(&x)?
        let x1 =
            Array2::from_shape_vec((4, 2), vec![1.0, 1.0, 2.0, 2.0, 8.0, 8.0, 9.0, 9.0]).unwrap();
        let y1 = array![0, 0, 1, 1];
        let x2 =
            Array2::from_shape_vec((4, 2), vec![0.5, 0.5, 1.5, 1.5, 7.5, 7.5, 8.5, 8.5]).unwrap();
        let y2 = array![0, 0, 1, 1];

        let clf = SGDClassifier::<f64>::new().with_random_state(42);
        let preds = clf
            .partial_fit(&x1, &y1)
            .unwrap()
            .partial_fit(&x2, &y2)
            .unwrap()
            .predict(&x1)
            .unwrap();
        assert_eq!(preds.len(), 4);
    }

    #[test]
    fn test_sgd_classifier_partial_fit_shape_mismatch() {
        let x =
            Array2::from_shape_vec((4, 2), vec![1.0, 1.0, 2.0, 2.0, 8.0, 8.0, 9.0, 9.0]).unwrap();
        let y = array![0, 0, 1, 1];
        let clf = SGDClassifier::<f64>::new().with_random_state(42);
        let fitted = clf.partial_fit(&x, &y).unwrap();

        let x_bad = Array2::from_shape_vec((2, 3), vec![1.0; 6]).unwrap();
        let y_bad = array![0, 1];
        assert!(fitted.partial_fit(&x_bad, &y_bad).is_err());
    }

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

        let clf = SGDClassifier::<f64>::new()
            .with_loss(ClassifierLoss::ModifiedHuber)
            .with_random_state(42)
            .with_max_iter(500);
        let fitted = clf.fit(&x, &y).unwrap();
        let preds = fitted.predict(&x).unwrap();
        assert_eq!(preds.len(), 6);
    }

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

        let clf = SGDClassifier::<f64>::new()
            .with_loss(ClassifierLoss::SquaredError)
            .with_random_state(42)
            .with_max_iter(500);
        let fitted = clf.fit(&x, &y).unwrap();
        let preds = fitted.predict(&x).unwrap();
        assert_eq!(preds.len(), 6);
    }

    #[test]
    fn test_sgd_classifier_pipeline() {
        let x = Array2::from_shape_vec((6, 1), vec![-3.0, -2.0, -1.0, 1.0, 2.0, 3.0]).unwrap();
        let y = Array1::from_vec(vec![0.0, 0.0, 0.0, 1.0, 1.0, 1.0]);

        let clf = SGDClassifier::<f64>::new().with_random_state(42);
        let fitted = clf.fit_pipeline(&x, &y).unwrap();
        let preds = fitted.predict_pipeline(&x).unwrap();
        assert_eq!(preds.len(), 6);
    }

    #[test]
    fn test_sgd_classifier_constant_lr() -> Result<(), FerroError> {
        let x = Array2::from_shape_vec((4, 1), vec![-2.0, -1.0, 1.0, 2.0]).unwrap_or_default();
        let y = array![0, 0, 1, 1];

        // The `constant` schedule requires `eta0 > 0`; the default is now 0.0
        // (sklearn `optimal` default), so set it explicitly.
        let clf = SGDClassifier::<f64>::new()
            .with_learning_rate(LearningRateSchedule::Constant)
            .with_eta0(0.01)
            .with_random_state(42)
            .with_max_iter(200);
        let fitted = clf.fit(&x, &y)?;
        assert_eq!(fitted.predict(&x)?.len(), 4);
        Ok(())
    }

    #[test]
    fn test_sgd_classifier_f32() {
        let x = Array2::from_shape_vec((4, 1), vec![-2.0f32, -1.0, 1.0, 2.0]).unwrap();
        let y = array![0_usize, 0, 1, 1];

        let clf = SGDClassifier::<f32>::new()
            .with_random_state(42)
            .with_max_iter(200);
        let fitted = clf.fit(&x, &y).unwrap();
        assert_eq!(fitted.predict(&x).unwrap().len(), 4);
    }

    // -----------------------------------------------------------------------
    // SGDRegressor tests
    // -----------------------------------------------------------------------

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

        let model = SGDRegressor::<f64>::new()
            .with_random_state(42)
            .with_max_iter(2000)
            .with_eta0(0.01)
            .with_alpha(0.0);
        let fitted = model.fit(&x, &y).unwrap();
        let preds = fitted.predict(&x).unwrap();

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

    #[test]
    fn test_sgd_regressor_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]; // Wrong length
        let model = SGDRegressor::<f64>::new();
        assert!(model.fit(&x, &y).is_err());
    }

    #[test]
    fn test_sgd_regressor_predict_shape_mismatch() {
        let x =
            Array2::from_shape_vec((4, 2), vec![1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 4.0, 4.0]).unwrap();
        let y = array![1.0, 2.0, 3.0, 4.0];
        let model = SGDRegressor::<f64>::new().with_random_state(42);
        let fitted = model.fit(&x, &y).unwrap();

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

    #[test]
    fn test_sgd_regressor_invalid_eta0() {
        let x = Array2::from_shape_vec((3, 1), vec![1.0, 2.0, 3.0]).unwrap();
        let y = array![1.0, 2.0, 3.0];
        let model = SGDRegressor::<f64>::new().with_eta0(-0.1);
        assert!(model.fit(&x, &y).is_err());
    }

    #[test]
    fn test_sgd_regressor_has_coefficients() {
        let x =
            Array2::from_shape_vec((4, 2), vec![1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 4.0, 4.0]).unwrap();
        let y = array![1.0, 2.0, 3.0, 4.0];
        let model = SGDRegressor::<f64>::new().with_random_state(42);
        let fitted = model.fit(&x, &y).unwrap();
        assert_eq!(fitted.coefficients().len(), 2);
    }

    #[test]
    fn test_sgd_regressor_partial_fit() {
        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 = SGDRegressor::<f64>::new().with_random_state(42);
        let fitted = model.partial_fit(&x, &y).unwrap();
        let fitted = fitted.partial_fit(&x, &y).unwrap();
        let preds = fitted.predict(&x).unwrap();
        assert_eq!(preds.len(), 4);
    }

    #[test]
    fn test_sgd_regressor_partial_fit_chain() {
        let x1 = Array2::from_shape_vec((3, 1), vec![1.0, 2.0, 3.0]).unwrap();
        let y1 = array![2.0, 4.0, 6.0];
        let x2 = Array2::from_shape_vec((3, 1), vec![4.0, 5.0, 6.0]).unwrap();
        let y2 = array![8.0, 10.0, 12.0];

        let model = SGDRegressor::<f64>::new().with_random_state(42);
        let preds = model
            .partial_fit(&x1, &y1)
            .unwrap()
            .partial_fit(&x2, &y2)
            .unwrap()
            .predict(&x1)
            .unwrap();
        assert_eq!(preds.len(), 3);
    }

    #[test]
    fn test_sgd_regressor_partial_fit_shape_mismatch() {
        let x = Array2::from_shape_vec((3, 2), vec![1.0, 1.0, 2.0, 2.0, 3.0, 3.0]).unwrap();
        let y = array![1.0, 2.0, 3.0];
        let model = SGDRegressor::<f64>::new().with_random_state(42);
        let fitted = model.partial_fit(&x, &y).unwrap();

        let x_bad = Array2::from_shape_vec((2, 3), vec![1.0; 6]).unwrap();
        let y_bad = array![1.0, 2.0];
        assert!(fitted.partial_fit(&x_bad, &y_bad).is_err());
    }

    #[test]
    fn test_sgd_regressor_huber_loss() {
        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 = SGDRegressor::<f64>::new()
            .with_loss(RegressorLoss::Huber(1.35))
            .with_random_state(42)
            .with_max_iter(500);
        let fitted = model.fit(&x, &y).unwrap();
        assert_eq!(fitted.predict(&x).unwrap().len(), 4);
    }

    #[test]
    fn test_sgd_regressor_epsilon_insensitive() {
        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 = SGDRegressor::<f64>::new()
            .with_loss(RegressorLoss::EpsilonInsensitive(0.1))
            .with_random_state(42)
            .with_max_iter(500);
        let fitted = model.fit(&x, &y).unwrap();
        assert_eq!(fitted.predict(&x).unwrap().len(), 4);
    }

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

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

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

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

    #[test]
    fn test_sgd_regressor_empty_data() {
        let x = Array2::<f64>::zeros((0, 2));
        let y = Array1::<f64>::zeros(0);
        let model = SGDRegressor::<f64>::new();
        assert!(model.fit(&x, &y).is_err());
    }

    #[test]
    fn test_sgd_classifier_empty_data() {
        let x = Array2::<f64>::zeros((0, 2));
        let y = Array1::<usize>::zeros(0);
        let clf = SGDClassifier::<f64>::new();
        assert!(clf.fit(&x, &y).is_err());
    }

    #[test]
    fn test_sgd_classifier_default() {
        // sklearn `SGDClassifier()` defaults (live oracle / `:1242-1244`):
        // learning_rate='optimal', eta0=0.0, alpha=0.0001, power_t=0.5.
        let clf = SGDClassifier::<f64>::default();
        assert!(matches!(clf.learning_rate, LearningRateSchedule::Optimal));
        assert!((clf.eta0 - 0.0).abs() < 1e-12);
        assert!((clf.alpha - 0.0001).abs() < 1e-12);
        assert!((clf.power_t - 0.5).abs() < 1e-12);
    }

    #[test]
    fn test_sgd_regressor_default() {
        let model = SGDRegressor::<f64>::default();
        assert!(model.eta0 > 0.0);
        assert!(model.alpha >= 0.0);
    }
}