//! ElasticNet regression (combined L1 and L2 regularization).
//!
//! This module provides [`ElasticNet`], which fits a linear model with a
//! blended L1/L2 regularization penalty using coordinate descent with
//! soft-thresholding:
//!
//! ```text
//! minimize (1/(2n)) * ||X @ w - y||^2
//! + alpha * l1_ratio * ||w||_1
//! + (alpha/2) * (1 - l1_ratio) * ||w||_2^2
//! ```
//!
//! When `l1_ratio = 1`, ElasticNet is equivalent to Lasso. When
//! `l1_ratio = 0`, it is equivalent to Ridge. Intermediate values produce
//! solutions that are both sparse (L1) and small in magnitude (L2).
//!
//! ## REQ status (per `.design/linear/elastic_net.md`, mirrors `sklearn/linear_model/_coordinate_descent.py` @ 1.5.2)
//!
//! Mirrors `sklearn.linear_model.ElasticNet` (`_coordinate_descent.py:729`). CD with the L1/L2
//! split `soft_threshold(Xⱼᵀr/n, α·l1_ratio)/(XⱼᵀXⱼ/n + α·(1−l1_ratio))` ≡ sklearn's
//! `l1_reg=α·l1_ratio·n` / `l2_reg=α·(1−l1_ratio)·n`, stopped on sklearn's relative-change +
//! dual-gap criterion (REQ-13). coef_/intercept_ AND `n_iter_` match the live oracle exactly
//! (coef_ ≤1e-7); default `l1_ratio=0.5` matches sklearn.
//!
//! | REQ | Status | Evidence |
//! |---|---|---|
//! | REQ-1 (CD ElasticNet fit, L1/L2 split) | SHIPPED | `Fit for ElasticNet`; converged coef/intercept match oracle <1e-5 over alpha∈{0.01,0.1,1}×l1_ratio∈{0.1,0.5,0.9}. Consumers: `RsElasticNet` in `ferrolearn-python`, `ElasticNetCV`. |
//! | REQ-2 (predict) | SHIPPED | `Predict for FittedElasticNet`. |
//! | REQ-3 (fit_intercept incl. false) | SHIPPED | centering. |
//! | REQ-4 (l1_ratio mixing; =1→Lasso, =0→L2) | SHIPPED | l1_ratio=1 ≡ Lasso; l1_ratio=0 ≡ sklearn ElasticNet L2; both match oracle. |
//! | REQ-5 (L1 sparsity) | SHIPPED | exact-zero support set bit-identical to sklearn. |
//! | REQ-6 (HasCoefficients) | SHIPPED | `HasCoefficients for FittedElasticNet`. |
//! | REQ-7 (alpha/l1_ratio validation; l1_ratio∈[0,1]) | SHIPPED | matches sklearn's `_parameter_constraints` (l1_ratio=0 accepted by the class; the auto-grid error is owned by elastic_net_cv). |
//! | REQ-8 (positive=True) | SHIPPED | `positive` field + `with_positive` builder; CD loop branches on `self.positive` to `soft_threshold_positive(rho_j, alpha_l1) / denominators[j]` (non-negative soft-threshold, L2 in the denominator unchanged), mirroring sklearn's `positive` param (`_coordinate_descent.py:800`) clip `if positive and tmp < 0: w[ii] = 0.0` (`_cd_fast.pyx:191-195`). Oracle test `elasticnet_positive_matches_sklearn` → coef `[1.13685345, 0.0]`, intercept `-5.96023707` (live sklearn 1.5.2, differs from unconstrained `[0.9081389, -1.7687475]`); `elasticnet_positive_false_unchanged` regression guard. |
//! | REQ-12 (n_iter_ / dual_gap_ attrs) | SHIPPED | `FittedElasticNet<F>` carries `n_iter`/`dual_gap` fields + `n_iter()`/`dual_gap()` getters, mirroring sklearn `ElasticNet.n_iter_` (`_coordinate_descent.py:827`) and `dual_gap_` (`:831`). `fn enet_dual_gap` computes the duality gap on the CD design (centered/raw) using sklearn's `_cd_fast.pyx:216-247` formula (`l1_reg = α·l1_ratio·n`, `beta = α·(1−l1_ratio)·n`, the `XtA = XᵀR − beta·w` term + `0.5·beta·(1+const²)·‖w‖²`) with a final `/n` mapping to the `(1/2n)` objective; reduces to `lasso_dual_gap` when `l1_ratio = 1` (`beta = 0`). With REQ-13's dual-gap stopping criterion now landed, `n_iter_`'s VALUE matches sklearn exactly (`n_iter_ == 16` at alpha=0.3, `== 19` at alpha=0.1 on the fixture); `dual_gap_` matches sklearn's formula/value (`0.00010575563` at `alpha=0.3, l1_ratio=0.5`). Verification: `cargo test -p ferrolearn-linear --lib elastic_net` (`enet_dual_gap_formula_matches_numpy`, `enet_fitted_dual_gap_and_n_iter`, `enet_fields_dont_change_coef`, `enet_dual_gap_stopping_matches_sklearn_coef_and_niter`). |
//! | REQ-13 (dual-gap stopping criterion) | SHIPPED | `Fit::fit for ElasticNet` now uses sklearn's two-level criterion (`_cd_fast.pyx:167-249`): `tol_scaled = tol·(target·target)` (`:167-168`), per sweep track `w_max`/`d_w_max`, gate on `w_max==0 || d_w_max/w_max < tol || last_iter` (`:207-211`), and inside the gate break only when the UN-normalized gap `enet_dual_gap(...)·n < tol_scaled` (`:249`) — `enet_dual_gap` already carries the L2/beta term. Matches sklearn's `coef_` to ≤1e-7 and `n_iter_` exactly (16 at alpha=0.3, 19 at alpha=0.1). Verification: `cargo test -p ferrolearn-linear --lib elastic_net` (`enet_dual_gap_stopping_matches_sklearn_coef_and_niter`, `enet_dual_gap_stopping_second_alpha`). |
//! | REQ-10 (selection='random' + random_state) | SHIPPED | Reuses `pub enum CoordSelection { Cyclic, Random }` from `lasso.rs` + `pub selection`/`pub random_state` fields on `ElasticNet` with `with_selection`/`with_random_state` builders, mirroring sklearn `ElasticNet(selection=..., random_state=...)` (`_coordinate_descent.py` `__init__`). `Fit::fit`'s CD loop visits `0..n_features` in order for `Cyclic` (BYTE-IDENTICAL to the prior cyclic path, so coef_/`n_iter_`/dual-gap stay unchanged) and shuffles a reused index `Vec` each sweep for `Random` via `StdRng::seed_from_u64(random_state.unwrap_or(0))` (sklearn `_cd_fast.pyx` `enet_coordinate_descent` `random` branch picks `ii` instead of `f_iter`); per-coordinate update math + dual-gap stopping (REQ-13) are unchanged. The ElasticNet optimum is unique, so `Random` converges to the same optimum (≈3e-4 from cyclic due to stopping-within-tol). Exact bit-match to sklearn's `selection='random'` is numpy-MT19937-RNG-blocked (Rust `StdRng` ≠ numpy MT), so the random path verifies convergence-to-the-unique-optimum, not bitwise sklearn parity; the cyclic default IS bit-exact. Verification: `cargo test -p ferrolearn-linear --lib elastic_net` (`enet_selection_cyclic_default_unchanged`, `enet_selection_random_converges_to_optimum`). |
//! | REQ-11 (precompute/Gram) | SHIPPED | `pub precompute: bool` field (default `false`) on `ElasticNet` + `with_precompute` builder, mirroring sklearn `ElasticNet(precompute=False)` (`_coordinate_descent.py:774`). When `true`, `Fit::fit` runs CD on the precomputed `Q = Xcᵀ Xc` / `q = Xcᵀ yc` with an incrementally-maintained `H = Q·w` (sklearn `_cd_fast.pyx enet_coordinate_descent_gram`); `tmp = (q[j]−H[j])/n + col_norms[j]·w[j] ≡` the direct path's `rho_j + (XⱼᵀXⱼ/n)·w_old` since `Xⱼᵀr = q[j]−(Q·w)[j]`, then `soft_threshold(tmp, α·l1_ratio)/(col_norm + α·(1−l1_ratio))` keeps the L2 term in the denominator, so it reaches the SAME unique optimum (to ~1e-10 fp reassociation) with the SAME coordinate order + dual-gap stopping. `precompute=false` (default) is the byte-identical direct path. Verification: `cargo test -p ferrolearn-linear --lib elastic_net` (`enet_precompute_matches_sklearn`, `enet_precompute_default_false_unchanged`, `enet_precompute_equals_direct`). |
//! | REQ-9 (warm_start) | SHIPPED | `ElasticNet<F>` carries `pub warm_start: bool` (default `false`) + `pub coef_init: Option<Array1<F>>` (default `None`) with `with_warm_start`/`with_coef_init` builders, mirroring sklearn `ElasticNet(warm_start=False)` (`_coordinate_descent.py:795`). R-DEV-4 adaptation: ferrolearn estimators are immutable value types — there is no mutable `self.coef_` carried across repeated `.fit()` calls like sklearn's mutable estimator (`_coordinate_descent.py:1062-1063` reuses `self.coef_` when `warm_start`), so the prior coefficient vector is supplied EXPLICITLY via `coef_init` (sklearn's path solver seeds the same way: `_coordinate_descent.py:648-651`, `coef_ = np.zeros(...)` when `coef_init is None` else `np.asfortranarray(coef_init, ...)`). In `Fit::fit`, when `warm_start && coef_init.is_some()` the init vector is length-validated (`ShapeMismatch` on mismatch) and `w` is cloned from it (the direct path also seeds `residual = y_work − X_work·w`; the Gram path's `H = Q·w` already derives from the actual `w`); otherwise `w = zeros` — BYTE-IDENTICAL to the cold path. The numerics are identical, only the CD start point changes, so warm-from-converged reaches the same unique optimum in fewer sweeps. Verification (live sklearn 1.5.2 oracle, R-CHAR-3): cold `ElasticNet(alpha=0.5, l1_ratio=0.5)` → coef `[0.7643620892, 1.2564536255]`, `n_iter_=14`; warm (refit from converged coef) → coef `[0.7642996441, 1.2564980309]`, `n_iter_=1`. Tests `enet_warm_start_from_converged_matches_sklearn`, `enet_warm_start_default_unchanged`, `enet_warm_start_none_coef_init_equals_cold`, `enet_warm_start_coef_init_wrong_len_errors`. |
//! | REQ-14..15 NOT-STARTED | MultiTaskElasticNet (#418), ferray substrate (#419). |
//! | REQ-16 (non-finite input rejected) | SHIPPED | `Fit::fit for ElasticNet` rejects any NaN/+/-inf in X or y BEFORE coordinate descent with `FerroError::InvalidParameter`, mirroring sklearn's `_validate_data(force_all_finite=True)` (`_coordinate_descent.py:980`, default `force_all_finite=True` → `check_array` raises `ValueError("Input X contains NaN.")` / `"... contains infinity ..."`). `.iter().any(|v| !v.is_finite())` catches both NaN and Inf; the finite path is byte-identical (the guard never fires on finite input). Verified vs the live sklearn 1.5.2 oracle (R-CHAR-3): `ElasticNet().fit` raises `ValueError` for NaN/+inf/-inf in X and NaN/inf in y (`tests/divergence_linear_nonfinite.rs::enet_*`). Non-test consumer: the existing `Fit::fit` / `RsElasticNet` / `ElasticNetCV` consumers. (#2256) |
//!
//! acto-critic: NO DIVERGENCE FOUND — coef/intercept grid parity, l1_ratio=1↔Lasso, l1_ratio=0↔L2,
//! sparsity support, default l1_ratio, and a badly-scaled-feature stress all match the live oracle.
//! Two states only per goal.md R-DEFER-2.
//!
//! # Examples
//!
//! ```
//! use ferrolearn_linear::ElasticNet;
//! use ferrolearn_core::{Fit, Predict};
//! use ndarray::{array, Array1, Array2};
//!
//! let model = ElasticNet::<f64>::new()
//! .with_alpha(0.1)
//! .with_l1_ratio(0.5);
//! let x = Array2::from_shape_vec((4, 1), vec![1.0, 2.0, 3.0, 4.0]).unwrap();
//! let y = array![2.0, 4.0, 6.0, 8.0];
//!
//! let fitted = model.fit(&x, &y).unwrap();
//! let preds = fitted.predict(&x).unwrap();
//! ```
use crate::lasso::CoordSelection;
use ferrolearn_core::error::FerroError;
use ferrolearn_core::introspection::HasCoefficients;
use ferrolearn_core::pipeline::{FittedPipelineEstimator, PipelineEstimator};
use ferrolearn_core::traits::{Fit, Predict};
use ndarray::{Array1, Array2, Axis, ScalarOperand};
use num_traits::{Float, FromPrimitive};
use rand::SeedableRng;
use rand::seq::SliceRandom;
/// ElasticNet regression (L1 + L2 regularized least squares).
///
/// Minimizes a combination of L1 and L2 penalties controlled by
/// `alpha` and `l1_ratio`. Uses coordinate descent with soft-thresholding
/// to handle the non-smooth L1 component.
///
/// # Type Parameters
///
/// - `F`: The floating-point type (`f32` or `f64`).
#[derive(Debug, Clone)]
pub struct ElasticNet<F> {
/// Overall regularization strength. Larger values enforce stronger
/// regularization.
pub alpha: F,
/// Mix between L1 and L2 regularization.
/// - `l1_ratio = 1.0` → pure Lasso (L1 only)
/// - `l1_ratio = 0.0` → pure Ridge (L2 only)
/// - `0.0 < l1_ratio < 1.0` → ElasticNet blend
pub l1_ratio: F,
/// Maximum number of coordinate descent iterations.
pub max_iter: usize,
/// Convergence tolerance on the maximum coefficient change per pass.
pub tol: F,
/// Whether to fit an intercept (bias) term.
pub fit_intercept: bool,
/// When `true`, constrain coefficients to be non-negative.
pub positive: bool,
/// When `true`, run coordinate descent on the precomputed Gram matrix
/// `Q = Xcᵀ Xc` and `q = Xcᵀ yc` instead of recomputing residuals each
/// pass.
///
/// Mirrors sklearn `ElasticNet(precompute=False)` (`_coordinate_descent.py:774`);
/// the Gram path runs sklearn's `enet_coordinate_descent_gram`
/// (`_cd_fast.pyx`) with the ElasticNet L2 term `α·(1−l1_ratio)` in the
/// denominator. Reaches the same unique optimum (differing only at
/// floating-point reassociation level, ~1e-10).
pub precompute: bool,
/// Order in which coordinates are visited each coordinate-descent sweep.
///
/// Mirrors sklearn `ElasticNet(selection=...)` (default `Cyclic`).
pub selection: CoordSelection,
/// Seed for the RNG used when `selection == CoordSelection::Random`.
///
/// Mirrors sklearn `ElasticNet(random_state=...)` (default `None`). `None`
/// falls back to seed `0`.
pub random_state: Option<u64>,
/// When `true`, initialize coordinate descent from [`ElasticNet::coef_init`]
/// (the prior solution) instead of zeros.
///
/// Mirrors sklearn `ElasticNet(warm_start=False)`
/// (`_coordinate_descent.py:795`), which "reuse[s] the solution of the
/// previous call to fit as initialization" (`:796`). In sklearn the prior
/// solution is the mutable estimator's own `self.coef_`, reused when
/// `warm_start` is set (`_coordinate_descent.py:1062-1063`: `if not
/// self.warm_start or not hasattr(self, "coef_"): coef_ = np.zeros(...)`).
///
/// R-DEV-4 adaptation: ferrolearn estimators are immutable value types —
/// there is no mutable `self.coef_` carried across repeated `.fit()` calls.
/// So the prior coefficient vector is supplied EXPLICITLY through
/// [`ElasticNet::coef_init`] rather than read off the estimator. The numerics
/// are identical: CD starts from `coef_init` instead of zeros.
pub warm_start: bool,
/// Explicit coordinate-descent initialization vector used when
/// [`ElasticNet::warm_start`] is `true` (the R-DEV-4 stand-in for sklearn's
/// reused `self.coef_`).
///
/// Mirrors the `coef_init` seed fed to the path solver
/// (`_coordinate_descent.py:648-651`: `coef_ = np.zeros(...)` when
/// `coef_init is None`, else `coef_ = np.asfortranarray(coef_init, ...)`).
/// `None` (the default) — or `warm_start == false` — initializes `w` to
/// zeros, the byte-identical cold-start path. When `Some`, its length must
/// equal `n_features`.
pub coef_init: Option<Array1<F>>,
}
impl<F: Float + FromPrimitive> ElasticNet<F> {
/// Create a new `ElasticNet` with default settings.
///
/// Defaults: `alpha = 1.0`, `l1_ratio = 0.5`, `max_iter = 1000`,
/// `tol = 1e-4`, `fit_intercept = true`.
#[must_use]
pub fn new() -> Self {
Self {
alpha: F::one(),
l1_ratio: F::from(0.5).unwrap(),
max_iter: 1000,
tol: F::from(1e-4).unwrap(),
fit_intercept: true,
positive: false,
precompute: false,
selection: CoordSelection::Cyclic,
random_state: None,
warm_start: false,
coef_init: None,
}
}
/// Set the overall regularization strength.
#[must_use]
pub fn with_alpha(mut self, alpha: F) -> Self {
self.alpha = alpha;
self
}
/// Set the L1/L2 mixing ratio.
///
/// Must be in `[0.0, 1.0]`. Values outside this range will be rejected
/// at fit time.
#[must_use]
pub fn with_l1_ratio(mut self, l1_ratio: F) -> Self {
self.l1_ratio = l1_ratio;
self
}
/// Set the maximum number of coordinate descent iterations.
#[must_use]
pub fn with_max_iter(mut self, max_iter: usize) -> Self {
self.max_iter = max_iter;
self
}
/// Set the convergence tolerance on maximum coefficient change.
#[must_use]
pub fn with_tol(mut self, tol: F) -> Self {
self.tol = tol;
self
}
/// Set whether to fit an intercept term.
#[must_use]
pub fn with_fit_intercept(mut self, fit_intercept: bool) -> Self {
self.fit_intercept = fit_intercept;
self
}
/// Set whether to constrain coefficients to be non-negative.
///
/// Mirrors `sklearn.linear_model.ElasticNet(positive=...)`.
#[must_use]
pub fn with_positive(mut self, positive: bool) -> Self {
self.positive = positive;
self
}
/// Set whether to run coordinate descent on the precomputed Gram matrix.
///
/// Mirrors `sklearn.linear_model.ElasticNet(precompute=...)`
/// (`_coordinate_descent.py:774`); `true` selects sklearn's
/// `enet_coordinate_descent_gram` (`_cd_fast.pyx`), with the ElasticNet
/// L2 term `α·(1−l1_ratio)` in the per-coordinate denominator.
#[must_use]
pub fn with_precompute(mut self, precompute: bool) -> Self {
self.precompute = precompute;
self
}
/// Set the coordinate-selection order for coordinate descent.
///
/// Mirrors `sklearn.linear_model.ElasticNet(selection=...)`.
#[must_use]
pub fn with_selection(mut self, selection: CoordSelection) -> Self {
self.selection = selection;
self
}
/// Set the RNG seed used when `selection == CoordSelection::Random`.
///
/// Mirrors `sklearn.linear_model.ElasticNet(random_state=...)`.
#[must_use]
pub fn with_random_state(mut self, seed: u64) -> Self {
self.random_state = Some(seed);
self
}
/// Enable/disable warm-start coordinate-descent initialization.
///
/// Mirrors `sklearn.linear_model.ElasticNet(warm_start=...)`
/// (`_coordinate_descent.py:795`): when `true`, "reuse the solution of the
/// previous call to fit as initialization". R-DEV-4: ferrolearn estimators
/// are immutable value types with no mutable `self.coef_` carried across
/// `.fit()` calls, so the prior solution is supplied explicitly via
/// [`ElasticNet::with_coef_init`]; `warm_start` only gates whether that
/// vector (when present) is used instead of zeros.
#[must_use]
pub fn with_warm_start(mut self, warm_start: bool) -> Self {
self.warm_start = warm_start;
self
}
/// Provide the explicit coordinate-descent initialization vector used when
/// [`ElasticNet::warm_start`] is `true`.
///
/// R-DEV-4 adaptation of sklearn's reused `self.coef_`
/// (`_coordinate_descent.py:1062-1063`, seeded into the path solver's
/// `coef_init` at `:648-651`): because ferrolearn estimators are immutable
/// value types, the prior coefficient vector is passed in explicitly rather
/// than read off a mutated estimator. Its length must equal `n_features` at
/// fit time, else [`Fit::fit`] returns [`FerroError::ShapeMismatch`].
#[must_use]
pub fn with_coef_init(mut self, coef: Array1<F>) -> Self {
self.coef_init = Some(coef);
self
}
}
impl<F: Float + FromPrimitive> Default for ElasticNet<F> {
fn default() -> Self {
Self::new()
}
}
/// Fitted ElasticNet regression model.
///
/// Stores the learned (potentially sparse) coefficients and intercept.
/// Implements [`Predict`] and [`HasCoefficients`].
#[derive(Debug, Clone)]
pub struct FittedElasticNet<F> {
/// Learned coefficient vector (some may be exactly zero when L1 > 0).
coefficients: Array1<F>,
/// Learned intercept (bias) term.
intercept: F,
/// Number of full coordinate-descent sweeps performed before
/// convergence/break (mirrors sklearn `ElasticNet.n_iter_`).
n_iter: usize,
/// Duality gap at the returned solution (mirrors sklearn `ElasticNet.dual_gap_`).
dual_gap: F,
}
impl<F: Float> FittedElasticNet<F> {
/// Returns the intercept (bias) term learned during fitting.
pub fn intercept(&self) -> F {
self.intercept
}
/// Number of coordinate-descent sweeps run by the solver.
///
/// Mirrors sklearn's `ElasticNet.n_iter_` attribute
/// (`_coordinate_descent.py:827`). ferrolearn uses sklearn's relative-change
/// and dual-gap stopping criterion (REQ-13, `_cd_fast.pyx:167-249`), so this
/// 1-based count matches sklearn's `n_iter_` value exactly at the same
/// optimum.
#[must_use]
pub fn n_iter(&self) -> usize {
self.n_iter
}
/// Duality gap at the returned solution, on the `(1/2n)`-scaled objective.
///
/// Mirrors sklearn's `ElasticNet.dual_gap_` attribute
/// (`_coordinate_descent.py:831`); computed by [`enet_dual_gap`] on the same
/// (centered/raw) design the coordinate descent solved.
#[must_use]
pub fn dual_gap(&self) -> F {
self.dual_gap
}
}
/// ElasticNet duality gap on the `(1/2n)`-scaled objective, mirroring sklearn's
/// `enet_coordinate_descent` gap (`_cd_fast.pyx:216-247`) with the final `/n`
/// mapping sklearn's un-normalized `(1/2)||y−Xw||² + l1_reg·||w||₁ +
/// (1/2)·l2_reg·||w||²` (`l1_reg = alpha·l1_ratio·n`, `l2_reg =
/// alpha·(1−l1_ratio)·n`, `_coordinate_descent.py:655-656`) back to ferrolearn's
/// `(1/2n)` scaling. Reduces to the Lasso gap when `l1_ratio = 1` (`beta = 0`).
///
/// `xc`/`yc` are the design the coordinate descent actually solved on
/// (centered when `fit_intercept`, raw otherwise); `w` is the fitted coef.
pub(crate) fn enet_dual_gap<F>(
xc: &Array2<F>,
yc: &Array1<F>,
w: &Array1<F>,
alpha: F,
l1_ratio: F,
) -> F
where
F: Float + ScalarOperand + 'static,
{
let n = xc.nrows();
let n_f = F::from(n).unwrap_or_else(F::one);
// R = yc − Xc·w
let residual = yc - &xc.dot(w);
// l1_reg = alpha · l1_ratio · n ; beta = alpha · (1 − l1_ratio) · n.
let l1_reg = alpha * l1_ratio * n_f;
let beta = alpha * (F::one() - l1_ratio) * n_f;
// XtA = Xcᵀ·R − beta·w ; dual_norm_XtA = max(|XtA[j]|).
let xt_a = xc.t().dot(&residual) - &(w * beta);
let dual_norm_xt_a = xt_a.iter().fold(F::zero(), |acc, &v| acc.max(v.abs()));
let r_norm2 = residual.dot(&residual);
let w_norm2 = w.dot(w);
let (const_factor, mut gap) = if dual_norm_xt_a > l1_reg {
let c = l1_reg / dual_norm_xt_a;
let half = F::from(0.5).unwrap_or_else(F::one);
(c, half * (r_norm2 + r_norm2 * c * c))
} else {
(F::one(), r_norm2)
};
// l1_norm = ‖w‖₁
let l1_norm = w.iter().fold(F::zero(), |acc, &wj| acc + wj.abs());
// R · yc
let r_dot_y = residual.dot(yc);
let half = F::from(0.5).unwrap_or_else(F::one);
gap = gap + l1_reg * l1_norm - const_factor * r_dot_y
+ half * beta * (F::one() + const_factor * const_factor) * w_norm2;
gap / n_f
}
/// Soft-thresholding operator used in coordinate descent for L1 penalty.
///
/// Returns `sign(x) * max(|x| - threshold, 0)`.
#[inline]
fn soft_threshold<F: Float>(x: F, threshold: F) -> F {
if x > threshold {
x - threshold
} else if x < -threshold {
x + threshold
} else {
F::zero()
}
}
/// Non-negative soft-thresholding operator for `positive=True` ElasticNet.
///
/// Returns `max(x - threshold, 0)`, dropping the negative branch so the
/// coordinate is never negative. Mirrors sklearn `_cd_fast.pyx:191-195`
/// (`if positive and tmp < 0: w[ii] = 0.0`); the L2 term lives in the
/// denominator and is unaffected.
#[inline]
fn soft_threshold_positive<F: Float>(x: F, threshold: F) -> F {
let z = x - threshold;
if z > F::zero() { z } else { F::zero() }
}
impl<F: Float + Send + Sync + ScalarOperand + FromPrimitive + 'static> Fit<Array2<F>, Array1<F>>
for ElasticNet<F>
{
type Fitted = FittedElasticNet<F>;
type Error = FerroError;
/// Fit the ElasticNet model using coordinate descent.
///
/// Centers the data if `fit_intercept` is `true`, then alternates
/// coordinate updates using the soft-threshold rule with L2 scaling.
///
/// # Errors
///
/// - [`FerroError::ShapeMismatch`] if `x` and `y` have different numbers
/// of samples.
/// - [`FerroError::InvalidParameter`] if `alpha` is negative, `l1_ratio`
/// is outside `[0, 1]`, or `tol` is non-positive.
/// - [`FerroError::InsufficientSamples`] if `n_samples == 0`.
fn fit(&self, x: &Array2<F>, y: &Array1<F>) -> Result<FittedElasticNet<F>, FerroError> {
let (n_samples, n_features) = x.dim();
if n_samples != y.len() {
return Err(FerroError::ShapeMismatch {
expected: vec![n_samples],
actual: vec![y.len()],
context: "y length must match number of samples in X".into(),
});
}
if self.alpha < F::zero() {
return Err(FerroError::InvalidParameter {
name: "alpha".into(),
reason: "must be non-negative".into(),
});
}
if self.l1_ratio < F::zero() || self.l1_ratio > F::one() {
return Err(FerroError::InvalidParameter {
name: "l1_ratio".into(),
reason: "must be in [0, 1]".into(),
});
}
if n_samples == 0 {
return Err(FerroError::InsufficientSamples {
required: 1,
actual: 0,
context: "ElasticNet requires at least one sample".into(),
});
}
// sklearn `ElasticNet.fit` -> `self._validate_data(X, y, ...)`
// (`_coordinate_descent.py:980`); the call keeps the default
// `force_all_finite=True`, so `check_array` rejects any NaN or +/-inf in
// X OR y with a `ValueError` BEFORE coordinate descent runs.
// `.iter().any(|v| !v.is_finite())` rejects both NaN and Inf (bounds-safe,
// no panic, R-CODE-2), matching the crate idiom (`multi_task_lasso.rs`).
// (#2256)
if x.iter().any(|v| !v.is_finite()) {
return Err(FerroError::InvalidParameter {
name: "X".into(),
reason: "Input X contains NaN or infinity.".into(),
});
}
if y.iter().any(|v| !v.is_finite()) {
return Err(FerroError::InvalidParameter {
name: "y".into(),
reason: "Input y contains NaN or infinity.".into(),
});
}
let n_f = F::from(n_samples).ok_or_else(|| FerroError::NumericalInstability {
message: "failed to convert n_samples to float".into(),
})?;
// Center data when fitting intercept.
let (x_work, y_work, x_mean, y_mean) = if self.fit_intercept {
let x_mean = x
.mean_axis(Axis(0))
.ok_or_else(|| FerroError::NumericalInstability {
message: "failed to compute column means".into(),
})?;
let y_mean = y.mean().ok_or_else(|| FerroError::NumericalInstability {
message: "failed to compute target mean".into(),
})?;
let x_c = x - &x_mean;
let y_c = y - y_mean;
(x_c, y_c, Some(x_mean), Some(y_mean))
} else {
(x.clone(), y.clone(), None, None)
};
// Precompute per-column X_j^T X_j / n (used as denominator).
let col_norms: Vec<F> = (0..n_features)
.map(|j| {
let col = x_work.column(j);
col.dot(&col) / n_f
})
.collect();
// L1 and L2 penalty strengths split from alpha/l1_ratio.
let alpha_l1 = self.alpha * self.l1_ratio;
let alpha_l2 = self.alpha * (F::one() - self.l1_ratio);
// Effective denominator per column: (X_j^T X_j / n) + alpha_l2.
let denominators: Vec<F> = col_norms.iter().map(|&cn| cn + alpha_l2).collect();
// Initialize coefficients. Cold start (default) is zeros; warm start
// reuses the explicit `coef_init` (the R-DEV-4 stand-in for sklearn's
// reused mutable `self.coef_`, `_coordinate_descent.py:1062-1063`/
// `:648-651`). `warm_start == false` or `coef_init == None` is the
// byte-identical zeros path.
let mut w = if self.warm_start
&& let Some(coef) = &self.coef_init
{
if coef.len() != n_features {
return Err(FerroError::ShapeMismatch {
expected: vec![n_features],
actual: vec![coef.len()],
context: "coef_init length must equal number of features".into(),
});
}
coef.clone()
} else {
Array1::<F>::zeros(n_features)
};
// Keep the (centered/raw) target for the final dual-gap computation;
// the CD loop consumes a working copy into `residual`.
let target = y_work.clone();
let mut residual = y_work;
// sklearn's stopping criterion (`_cd_fast.pyx:144-249`):
// - `d_w_tol = tol` is the UN-scaled relative-change gate (`:144`);
// - `tol_scaled = tol · (target·target)` is the gap threshold,
// `tol *= dot(y, y)` at `:167-168` (`target` is the centered/raw
// target the CD actually solves on).
let d_w_tol = self.tol;
let tol_scaled = self.tol * target.dot(&target);
// For `selection == Random`, build the RNG ONCE before the sweep loop
// and reuse a reusable index buffer; each sweep shuffles the visiting
// order (sklearn `_cd_fast.pyx` `enet_coordinate_descent` `random`
// branch picks `ii` via `rand_int` instead of the cyclic `f_iter`).
// `Cyclic` keeps the byte-identical `0..n_features` order, so the
// per-coordinate update math AND the dual-gap stopping criterion
// (REQ-13) stay unchanged.
let mut rng = rand::rngs::StdRng::seed_from_u64(self.random_state.unwrap_or(0));
let mut order: Vec<usize> = (0..n_features).collect();
let mut n_iter = 0_usize;
// REQ-11: Gram (precompute) coordinate-descent path. Mirrors sklearn's
// `enet_coordinate_descent_gram` (`_cd_fast.pyx`): run CD on the
// precomputed `Q = Xcᵀ Xc` and `q = Xcᵀ yc`, maintaining `H = Q·w`
// incrementally instead of recomputing residuals each sweep.
// Algebraically identical to the direct path (`Xⱼᵀr = q[j] − (Q·w)[j]`),
// so it converges to the same unique optimum (to fp reassociation,
// ~1e-10). Keeps the SAME `(1/n)` normalization, the same L1 threshold
// `alpha_l1` and L2-in-the-denominator `alpha_l2`, the same coordinate
// visiting order, and the same dual-gap stopping criterion as the direct
// path so `n_iter_` matches.
if self.precompute {
// Q = Xcᵀ Xc (n_features × n_features); q = Xcᵀ yc (here `residual`
// still equals the centered/raw target — it is not yet adjusted for
// a warm-start `w` since the Gram path tracks `H = Q·w` instead).
let gram = x_work.t().dot(&x_work);
let q = x_work.t().dot(&residual);
// H = Q·w (zeros for a cold start where `w == 0`; the actual `Q·w`
// for a warm start, so `tmp = (q[j] − H[j])/n + col_norms[j]·w[j]`
// is correct from the first sweep regardless of the init).
let mut h = gram.dot(&w);
for iter in 0..self.max_iter {
n_iter = iter + 1;
let mut w_max = F::zero();
let mut d_w_max = F::zero();
if self.selection == CoordSelection::Random {
order.shuffle(&mut rng);
}
for &j in &order {
let w_old = w[j];
// tmp ≡ direct `rho_j` + (XⱼᵀXⱼ/n)·w_old:
// (q[j] − H[j])/n + col_norms[j]·w[j], since
// Xⱼᵀr = q[j] − (Q·w)[j] and col_norms[j] = Q[j,j]/n.
let tmp = (q[j] - h[j]) / n_f + col_norms[j] * w_old;
// Soft-threshold for L1 (alpha_l1), then divide by the
// ElasticNet denominator (col_norm + alpha_l2). Identical
// to the direct path's per-coordinate update, just
// Gram-sourced.
let w_new = if denominators[j] > F::zero() {
let thresholded = if self.positive {
soft_threshold_positive(tmp, alpha_l1)
} else {
soft_threshold(tmp, alpha_l1)
};
thresholded / denominators[j]
} else {
F::zero()
};
if w_new != w_old {
// H += (w_new − w_old) · Q.column(j).
let delta = w_new - w_old;
let col = gram.column(j);
for i in 0..n_features {
h[i] = h[i] + delta * col[i];
}
}
let change = (w_new - w_old).abs();
if change > d_w_max {
d_w_max = change;
}
if w_new.abs() > w_max {
w_max = w_new.abs();
}
w[j] = w_new;
}
// SAME dual-gap stopping as the direct path: reuse the
// residual-based `enet_dual_gap` on (x_work, target) — equal to
// the Gram gap to fp precision, so `n_iter_` matches.
let last_iter = iter == self.max_iter - 1;
if w_max == F::zero() || d_w_max / w_max < d_w_tol || last_iter {
let dual_gap = enet_dual_gap(&x_work, &target, &w, self.alpha, self.l1_ratio);
let gap_raw = dual_gap * n_f;
if gap_raw < tol_scaled {
let intercept = compute_intercept(&x_mean, &y_mean, &w);
return Ok(FittedElasticNet {
coefficients: w,
intercept,
n_iter,
dual_gap,
});
}
}
}
// Did not converge within max_iter; return the current solution.
let intercept = compute_intercept(&x_mean, &y_mean, &w);
let dual_gap = enet_dual_gap(&x_work, &target, &w, self.alpha, self.l1_ratio);
return Ok(FittedElasticNet {
coefficients: w,
intercept,
n_iter,
dual_gap,
});
}
// Direct path: the CD loop maintains `residual = y_work − X_work·w`,
// adding back `X_j·w_old` per coordinate before recomputing `rho_j`. With
// a non-zero warm-start `w`, seed the residual with that running
// contribution removed. For the cold path (`w == 0`) `X_work·w` is the
// zero vector and the subtraction is a byte-identical no-op, so this is
// gated on warm start to leave the default path provably untouched.
if self.warm_start && self.coef_init.is_some() {
residual = &residual - &x_work.dot(&w);
}
for iter in 0..self.max_iter {
n_iter = iter + 1;
let mut w_max = F::zero();
let mut d_w_max = F::zero();
if self.selection == CoordSelection::Random {
order.shuffle(&mut rng);
}
for &j in &order {
let col_j = x_work.column(j);
let w_old = w[j];
// Add back contribution of current coefficient j to residual.
if w_old != F::zero() {
for i in 0..n_samples {
residual[i] = residual[i] + col_j[i] * w_old;
}
}
// Unpenalized correlation: X_j^T r / n.
let rho_j = col_j.dot(&residual) / n_f;
// Apply soft-threshold for L1, then divide by (col_norm + alpha_l2).
// For `positive=True`, use the non-negative soft-threshold so the
// coefficient is never negative (sklearn `_cd_fast.pyx:191-195`); the
// L2 term in the denominator is unchanged.
let w_new = if denominators[j] > F::zero() {
let thresholded = if self.positive {
soft_threshold_positive(rho_j, alpha_l1)
} else {
soft_threshold(rho_j, alpha_l1)
};
thresholded / denominators[j]
} else {
F::zero()
};
// Update residual with new coefficient.
if w_new != F::zero() {
for i in 0..n_samples {
residual[i] = residual[i] - col_j[i] * w_new;
}
}
// Track the largest coordinate update and the largest
// coefficient magnitude this sweep (`_cd_fast.pyx:201-205`).
let change = (w_new - w_old).abs();
if change > d_w_max {
d_w_max = change;
}
if w_new.abs() > w_max {
w_max = w_new.abs();
}
w[j] = w_new;
}
// sklearn's two-level convergence gate (`_cd_fast.pyx:207-251`):
// only when coordinates barely moved (relative gate) or on the
// last iteration do we compute the (expensive) dual gap, and we
// break only if the UN-normalized gap clears `tol · (target·target)`.
let last_iter = iter == self.max_iter - 1;
if w_max == F::zero() || d_w_max / w_max < d_w_tol || last_iter {
// `enet_dual_gap` returns the gap divided by `n` (the
// `dual_gap_` attribute scaling, REQ-12); multiply back to the
// un-normalized objective sklearn compares against
// `tol · (target·target)` (`:249`). The L2/beta term is already
// included in `enet_dual_gap`.
let dual_gap = enet_dual_gap(&x_work, &target, &w, self.alpha, self.l1_ratio);
let gap_raw = dual_gap * n_f;
if gap_raw < tol_scaled {
let intercept = compute_intercept(&x_mean, &y_mean, &w);
return Ok(FittedElasticNet {
coefficients: w,
intercept,
n_iter,
dual_gap,
});
}
}
}
// Return best solution found even without full convergence.
let intercept = compute_intercept(&x_mean, &y_mean, &w);
let dual_gap = enet_dual_gap(&x_work, &target, &w, self.alpha, self.l1_ratio);
Ok(FittedElasticNet {
coefficients: w,
intercept,
n_iter,
dual_gap,
})
}
}
/// Compute intercept from the centered means and fitted coefficients.
fn compute_intercept<F: Float + 'static>(
x_mean: &Option<Array1<F>>,
y_mean: &Option<F>,
w: &Array1<F>,
) -> F {
if let (Some(xm), Some(ym)) = (x_mean, y_mean) {
*ym - xm.dot(w)
} else {
F::zero()
}
}
impl<F: Float + Send + Sync + ScalarOperand + 'static> Predict<Array2<F>> for FittedElasticNet<F> {
type Output = Array1<F>;
type Error = FerroError;
/// Predict target values for the given feature matrix.
///
/// Computes `X @ coefficients + intercept`.
///
/// # Errors
///
/// Returns [`FerroError::ShapeMismatch`] if the number of features
/// does not match the fitted model.
fn predict(&self, x: &Array2<F>) -> Result<Array1<F>, FerroError> {
let n_features = x.ncols();
if n_features != self.coefficients.len() {
return Err(FerroError::ShapeMismatch {
expected: vec![self.coefficients.len()],
actual: vec![n_features],
context: "number of features must match fitted model".into(),
});
}
let preds = x.dot(&self.coefficients) + self.intercept;
Ok(preds)
}
}
impl<F: Float + Send + Sync + ScalarOperand + 'static> HasCoefficients<F> for FittedElasticNet<F> {
/// Returns the learned coefficient vector.
fn coefficients(&self) -> &Array1<F> {
&self.coefficients
}
/// Returns the learned intercept term.
fn intercept(&self) -> F {
self.intercept
}
}
// Pipeline integration.
impl<F> PipelineEstimator<F> for ElasticNet<F>
where
F: Float + FromPrimitive + ScalarOperand + Send + Sync + 'static,
{
/// Fit the model and return it as a boxed pipeline estimator.
///
/// # Errors
///
/// Propagates any [`FerroError`] from `fit`.
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 FittedElasticNet<F>
where
F: Float + ScalarOperand + Send + Sync + 'static,
{
/// Generate predictions via the pipeline interface.
///
/// # Errors
///
/// Propagates any [`FerroError`] from `predict`.
fn predict_pipeline(&self, x: &Array2<F>) -> Result<Array1<F>, FerroError> {
self.predict(x)
}
}
#[cfg(test)]
mod tests {
use super::*;
use approx::assert_relative_eq;
use ndarray::array;
// ---- soft_threshold helpers ----
#[test]
fn test_soft_threshold_positive() {
assert_relative_eq!(soft_threshold(5.0_f64, 1.0), 4.0);
}
#[test]
fn test_soft_threshold_negative() {
assert_relative_eq!(soft_threshold(-5.0_f64, 1.0), -4.0);
}
#[test]
fn test_soft_threshold_within_band() {
assert_relative_eq!(soft_threshold(0.5_f64, 1.0), 0.0);
assert_relative_eq!(soft_threshold(-0.5_f64, 1.0), 0.0);
assert_relative_eq!(soft_threshold(0.0_f64, 1.0), 0.0);
}
// ---- Builder ----
#[test]
fn test_default_builder() {
let m = ElasticNet::<f64>::new();
assert_relative_eq!(m.alpha, 1.0);
assert_relative_eq!(m.l1_ratio, 0.5);
assert_eq!(m.max_iter, 1000);
assert!(m.fit_intercept);
}
#[test]
fn test_builder_setters() {
let m = ElasticNet::<f64>::new()
.with_alpha(0.5)
.with_l1_ratio(0.2)
.with_max_iter(500)
.with_tol(1e-6)
.with_fit_intercept(false);
assert_relative_eq!(m.alpha, 0.5);
assert_relative_eq!(m.l1_ratio, 0.2);
assert_eq!(m.max_iter, 500);
assert!(!m.fit_intercept);
}
// ---- Validation errors ----
#[test]
fn test_negative_alpha_error() {
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 result = ElasticNet::<f64>::new().with_alpha(-1.0).fit(&x, &y);
assert!(result.is_err());
}
#[test]
fn test_l1_ratio_out_of_range_error() {
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 result = ElasticNet::<f64>::new().with_l1_ratio(1.5).fit(&x, &y);
assert!(result.is_err());
}
#[test]
fn test_shape_mismatch_error() {
let x = Array2::from_shape_vec((3, 1), vec![1.0, 2.0, 3.0]).unwrap();
let y = array![1.0, 2.0];
let result = ElasticNet::<f64>::new().fit(&x, &y);
assert!(result.is_err());
}
// ---- Correctness ----
#[test]
fn test_lasso_limit_l1_ratio_one() {
// With l1_ratio=1, ElasticNet should behave like Lasso.
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 = ElasticNet::<f64>::new().with_alpha(0.0).with_l1_ratio(1.0);
let fitted = model.fit(&x, &y).unwrap();
assert_relative_eq!(fitted.coefficients()[0], 2.0, epsilon = 1e-4);
assert_relative_eq!(fitted.intercept(), 1.0, epsilon = 1e-4);
}
#[test]
fn test_ridge_limit_l1_ratio_zero() {
// With l1_ratio=0 and alpha=0, should recover OLS.
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 = ElasticNet::<f64>::new().with_alpha(0.0).with_l1_ratio(0.0);
let fitted = model.fit(&x, &y).unwrap();
assert_relative_eq!(fitted.coefficients()[0], 2.0, epsilon = 1e-4);
assert_relative_eq!(fitted.intercept(), 1.0, epsilon = 1e-4);
}
#[test]
fn test_sparsity_with_high_l1_ratio() {
// High alpha with l1_ratio=1 should zero out irrelevant features.
let x = Array2::from_shape_vec(
(10, 3),
vec![
1.0, 0.0, 0.0, 2.0, 0.0, 0.0, 3.0, 0.0, 0.0, 4.0, 0.0, 0.0, 5.0, 0.0, 0.0, 6.0,
0.0, 0.0, 7.0, 0.0, 0.0, 8.0, 0.0, 0.0, 9.0, 0.0, 0.0, 10.0, 0.0, 0.0,
],
)
.unwrap();
let y = array![2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0];
let model = ElasticNet::<f64>::new().with_alpha(5.0).with_l1_ratio(1.0);
let fitted = model.fit(&x, &y).unwrap();
assert_relative_eq!(fitted.coefficients()[1], 0.0, epsilon = 1e-10);
assert_relative_eq!(fitted.coefficients()[2], 0.0, epsilon = 1e-10);
}
#[test]
fn test_higher_alpha_shrinks_more() {
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 low = ElasticNet::<f64>::new()
.with_alpha(0.01)
.with_l1_ratio(0.5)
.fit(&x, &y)
.unwrap();
let high = ElasticNet::<f64>::new()
.with_alpha(2.0)
.with_l1_ratio(0.5)
.fit(&x, &y)
.unwrap();
assert!(high.coefficients()[0].abs() <= low.coefficients()[0].abs());
}
#[test]
fn test_no_intercept() {
let x = Array2::from_shape_vec((4, 1), vec![1.0, 2.0, 3.0, 4.0]).unwrap();
let y = array![2.0, 4.0, 6.0, 8.0];
let fitted = ElasticNet::<f64>::new()
.with_alpha(0.0)
.with_l1_ratio(0.5)
.with_fit_intercept(false)
.fit(&x, &y)
.unwrap();
assert_relative_eq!(fitted.intercept(), 0.0, epsilon = 1e-10);
}
#[test]
fn test_predict_correct_length() {
let x = Array2::from_shape_vec((4, 1), vec![1.0, 2.0, 3.0, 4.0]).unwrap();
let y = array![2.0, 4.0, 6.0, 8.0];
let fitted = ElasticNet::<f64>::new()
.with_alpha(0.01)
.fit(&x, &y)
.unwrap();
let preds = fitted.predict(&x).unwrap();
assert_eq!(preds.len(), 4);
}
#[test]
fn test_predict_feature_mismatch() {
let x_train = Array2::from_shape_vec((3, 2), vec![1.0, 0.0, 2.0, 0.0, 3.0, 0.0]).unwrap();
let y = array![1.0, 2.0, 3.0];
let fitted = ElasticNet::<f64>::new()
.with_alpha(0.01)
.fit(&x_train, &y)
.unwrap();
let x_bad = Array2::from_shape_vec((3, 1), vec![1.0, 2.0, 3.0]).unwrap();
let result = fitted.predict(&x_bad);
assert!(result.is_err());
}
#[test]
fn test_has_coefficients_length() {
let x = Array2::from_shape_vec((3, 2), vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]).unwrap();
let y = array![1.0, 2.0, 3.0];
let fitted = ElasticNet::<f64>::new()
.with_alpha(0.1)
.fit(&x, &y)
.unwrap();
assert_eq!(fitted.coefficients().len(), 2);
}
#[test]
fn test_pipeline_integration() {
let x = Array2::from_shape_vec((4, 1), vec![1.0, 2.0, 3.0, 4.0]).unwrap();
let y = array![3.0, 5.0, 7.0, 9.0];
let model = ElasticNet::<f64>::new().with_alpha(0.01);
let fitted_pipe = model.fit_pipeline(&x, &y).unwrap();
let preds = fitted_pipe.predict_pipeline(&x).unwrap();
assert_eq!(preds.len(), 4);
}
// ---- positive=True (REQ-8) ----
#[test]
fn test_soft_threshold_positive_helper() {
// Non-negative branch: max(x - t, 0). Negative side clamps to 0.
assert_relative_eq!(soft_threshold_positive(5.0_f64, 1.0), 4.0);
assert_relative_eq!(soft_threshold_positive(-5.0_f64, 1.0), 0.0);
assert_relative_eq!(soft_threshold_positive(0.5_f64, 1.0), 0.0);
assert_relative_eq!(soft_threshold_positive(-0.5_f64, 1.0), 0.0);
assert_relative_eq!(soft_threshold_positive(0.0_f64, 1.0), 0.0);
}
/// Oracle fixture from live sklearn 1.5.2 (R-CHAR-3):
/// `X = [[1,3],[2,1],[3,4],[4,2],[5,5],[6,1],[2,4],[5,2]]`,
/// `y = X[:,0] - 2*X[:,1] + noise`.
fn positive_oracle_fixture() -> (Array2<f64>, Array1<f64>) {
let x: Array2<f64> = array![
[1.0, 3.0],
[2.0, 1.0],
[3.0, 4.0],
[4.0, 2.0],
[5.0, 5.0],
[6.0, 1.0],
[2.0, 4.0],
[5.0, 2.0],
];
let noise = array![0.1, -0.2, 0.15, 0.0, -0.1, 0.05, 0.2, -0.05];
let y: Array1<f64> = (0..8)
.map(|i| 1.0 * x[[i, 0]] - 2.0 * x[[i, 1]] + noise[i])
.collect();
(x, y)
}
#[test]
fn elasticnet_positive_matches_sklearn() {
// Live sklearn 1.5.2 oracle:
// ElasticNet(alpha=0.3, l1_ratio=0.5, positive=True)
// -> coef_ [1.13685345, 0.0], intercept_ -5.96023707
// (unconstrained ElasticNet(alpha=0.3, l1_ratio=0.5)
// -> coef_ [0.9081389, -1.7687475], intercept_ -0.29568051).
let (x, y) = positive_oracle_fixture();
let fit_res = ElasticNet::<f64>::new()
.with_alpha(0.3)
.with_l1_ratio(0.5)
.with_positive(true)
.fit(&x, &y);
assert!(fit_res.is_ok(), "positive fit should succeed");
let fitted = match fit_res {
Ok(f) => f,
Err(_) => return,
};
let coef = fitted.coefficients();
assert_relative_eq!(coef[0], 1.13685345, epsilon = 1e-5);
assert_relative_eq!(coef[1], 0.0, epsilon = 1e-5);
assert_relative_eq!(fitted.intercept(), -5.96023707, epsilon = 1e-4);
// All coefficients are non-negative.
for &c in coef.iter() {
assert!(c >= 0.0, "coefficient {c} should be non-negative");
}
// Differs materially from the unconstrained solution (~1.77 gap on
// feature 1), confirming the constraint is non-tautological.
let unc_res = ElasticNet::<f64>::new()
.with_alpha(0.3)
.with_l1_ratio(0.5)
.fit(&x, &y);
assert!(unc_res.is_ok(), "unconstrained fit should succeed");
let unconstrained = match unc_res {
Ok(f) => f,
Err(_) => return,
};
assert!((coef[1] - unconstrained.coefficients()[1]).abs() > 1.0);
}
#[test]
fn elasticnet_positive_false_unchanged() {
// positive=false (default) must be byte-identical to the plain fit.
let (x, y) = positive_oracle_fixture();
let default_res = ElasticNet::<f64>::new()
.with_alpha(0.3)
.with_l1_ratio(0.5)
.fit(&x, &y);
assert!(default_res.is_ok(), "default fit should succeed");
let default_fit = match default_res {
Ok(f) => f,
Err(_) => return,
};
let false_res = ElasticNet::<f64>::new()
.with_alpha(0.3)
.with_l1_ratio(0.5)
.with_positive(false)
.fit(&x, &y);
assert!(
false_res.is_ok(),
"explicit positive=false fit should succeed"
);
let explicit_false = match false_res {
Ok(f) => f,
Err(_) => return,
};
assert_eq!(
default_fit.coefficients(),
explicit_false.coefficients(),
"positive=false must be byte-identical to the default fit"
);
assert_eq!(default_fit.intercept(), explicit_false.intercept());
}
// ---- n_iter_ / dual_gap_ (REQ-12) ----
/// Centered fixture for the dual-gap oracle (R-CHAR-3):
/// `X = [[1,2],[2,1],[3,4],[4,3],[5,5]]`, `y = [3,2.5,7.1,6,11.2]`,
/// centered by column mean / target mean (the design the CD solves under
/// `fit_intercept`).
fn centered_dual_gap_fixture() -> Option<(Array2<f64>, Array1<f64>)> {
let x: Array2<f64> = array![[1.0, 2.0], [2.0, 1.0], [3.0, 4.0], [4.0, 3.0], [5.0, 5.0],];
let y: Array1<f64> = array![3.0, 2.5, 7.1, 6.0, 11.2];
let x_mean = x.mean_axis(Axis(0))?;
let y_mean = y.mean()?;
Some((&x - &x_mean, &y - y_mean))
}
fn raw_dual_gap_fixture() -> (Array2<f64>, Array1<f64>) {
let x: Array2<f64> = array![[1.0, 2.0], [2.0, 1.0], [3.0, 4.0], [4.0, 3.0], [5.0, 5.0],];
let y: Array1<f64> = array![3.0, 2.5, 7.1, 6.0, 11.2];
(x, y)
}
#[test]
fn enet_dual_gap_formula_matches_numpy() {
// numpy/sklearn-computed oracle points (NOT from ferrolearn),
// alpha=0.3, l1_ratio=0.5:
// gap(w=[0.5,1.0]) = 0.6369296296 (far-from-optimum)
// gap(w=[0.77323348,1.35480299]) = 0.0001057556 (the optimum)
let (xc, yc) = match centered_dual_gap_fixture() {
Some(f) => f,
None => return,
};
let far = enet_dual_gap(&xc, &yc, &array![0.5, 1.0], 0.3, 0.5);
assert_relative_eq!(far, 0.6369296296, epsilon = 1e-5);
let opt = enet_dual_gap(&xc, &yc, &array![0.77323348, 1.35480299], 0.3, 0.5);
assert_relative_eq!(opt, 0.0001057556, epsilon = 1e-7);
}
#[test]
fn enet_fitted_dual_gap_and_n_iter() {
// ElasticNet(alpha=0.3, l1_ratio=0.5) on the same fixture: dual_gap_
// converged near sklearn's 0.000106; n_iter_ within [1, max_iter].
let (x, y) = raw_dual_gap_fixture();
let fit_res = ElasticNet::<f64>::new()
.with_alpha(0.3)
.with_l1_ratio(0.5)
.fit(&x, &y);
assert!(fit_res.is_ok(), "fit should succeed");
let fitted = match fit_res {
Ok(f) => f,
Err(_) => return,
};
let gap = fitted.dual_gap();
assert!(gap >= 0.0, "dual_gap should be non-negative, got {gap}");
assert!(gap < 1e-3, "dual_gap should be converged-small, got {gap}");
let n_iter = fitted.n_iter();
assert!(n_iter >= 1, "n_iter should be at least 1, got {n_iter}");
assert!(n_iter <= 1000, "n_iter should be <= max_iter, got {n_iter}");
}
#[test]
fn enet_fields_dont_change_coef() {
// Regression guard: the additive n_iter_/dual_gap_ fields must not
// perturb coef_/intercept_. Compared against sklearn's converged
// coef_ = [0.77323348, 1.35480299]: with REQ-13's dual-gap stopping
// criterion the stop point is identical to sklearn, so the comparison
// is tight (1e-7) — matching sklearn BETTER, never loosened.
let (x, y) = raw_dual_gap_fixture();
let fit_res = ElasticNet::<f64>::new()
.with_alpha(0.3)
.with_l1_ratio(0.5)
.fit(&x, &y);
assert!(fit_res.is_ok(), "fit should succeed");
let fitted = match fit_res {
Ok(f) => f,
Err(_) => return,
};
assert_relative_eq!(fitted.coefficients()[0], 0.77323348, epsilon = 1e-7);
assert_relative_eq!(fitted.coefficients()[1], 1.35480299, epsilon = 1e-7);
}
#[test]
fn enet_dual_gap_stopping_matches_sklearn_coef_and_niter() {
// REQ-13: sklearn's relative-change + dual-gap stopping criterion.
// Live sklearn 1.5.2 oracle (R-CHAR-3):
// X=[[1,2],[2,1],[3,4],[4,3],[5,5]], y=[3,2.5,7.1,6,11.2]
// ElasticNet(alpha=0.3, l1_ratio=0.5).fit(X,y)
// -> coef_=[0.77323348, 1.35480299], n_iter_=16,
// dual_gap_=0.00010575563
let (x, y) = raw_dual_gap_fixture();
let fit_res = ElasticNet::<f64>::new()
.with_alpha(0.3)
.with_l1_ratio(0.5)
.fit(&x, &y);
assert!(fit_res.is_ok(), "fit should succeed");
let fitted = match fit_res {
Ok(f) => f,
Err(_) => return,
};
// Coef matches sklearn TIGHTLY now that the stopping point is identical.
assert_relative_eq!(fitted.coefficients()[0], 0.77323348, epsilon = 1e-7);
assert_relative_eq!(fitted.coefficients()[1], 1.35480299, epsilon = 1e-7);
// n_iter_ matches sklearn's 1-based dual-gap iteration count exactly.
assert_eq!(fitted.n_iter(), 16, "n_iter_ must match sklearn's 16");
// dual_gap_ (the /n attribute) stays the REQ-12 value.
assert_relative_eq!(fitted.dual_gap(), 0.00010575563, epsilon = 1e-7);
}
#[test]
fn enet_dual_gap_stopping_second_alpha() {
// Generalization check at alpha=0.1 (live sklearn 1.5.2 oracle):
// ElasticNet(alpha=0.1, l1_ratio=0.5).fit(X,y)
// -> coef_=[0.76514609, 1.47598354], n_iter_=19,
// dual_gap_=9.422349e-05
let (x, y) = raw_dual_gap_fixture();
let fit_res = ElasticNet::<f64>::new()
.with_alpha(0.1)
.with_l1_ratio(0.5)
.fit(&x, &y);
assert!(fit_res.is_ok(), "fit should succeed");
let fitted = match fit_res {
Ok(f) => f,
Err(_) => return,
};
assert_relative_eq!(fitted.coefficients()[0], 0.76514609, epsilon = 1e-7);
assert_relative_eq!(fitted.coefficients()[1], 1.47598354, epsilon = 1e-7);
assert_eq!(fitted.n_iter(), 19, "n_iter_ must match sklearn's 19");
assert_relative_eq!(fitted.dual_gap(), 9.422349e-05, epsilon = 1e-7);
}
// ---- selection='random' + random_state (REQ-10) ----
/// Oracle fixture for the selection tests (R-CHAR-3, live sklearn 1.5.2):
/// `X = [[1,2],[2,1],[3,4],[4,3],[5,5]]`, `y = [3,2.5,7.1,6,11.2]`,
/// `alpha=0.3`, `l1_ratio=0.5`.
fn selection_fixture() -> (Array2<f64>, Array1<f64>) {
let x: Array2<f64> = array![[1.0, 2.0], [2.0, 1.0], [3.0, 4.0], [4.0, 3.0], [5.0, 5.0],];
let y: Array1<f64> = array![3.0, 2.5, 7.1, 6.0, 11.2];
(x, y)
}
#[test]
fn enet_selection_cyclic_default_unchanged() {
// Default ElasticNet selection is Cyclic; coef must stay byte-identical
// to the prior cyclic path. Live sklearn 1.5.2 oracle (R-CHAR-3):
// ElasticNet(alpha=0.3, l1_ratio=0.5, selection='cyclic')
// -> coef_ [0.77323348, 1.35480299].
let (x, y) = selection_fixture();
// Default selection is Cyclic.
assert_eq!(ElasticNet::<f64>::new().selection, CoordSelection::Cyclic);
let default_res = ElasticNet::<f64>::new()
.with_alpha(0.3)
.with_l1_ratio(0.5)
.fit(&x, &y);
assert!(default_res.is_ok(), "default fit should succeed");
let default_fit = match default_res {
Ok(f) => f,
Err(_) => return,
};
// Matches sklearn's cyclic oracle tightly.
assert_relative_eq!(default_fit.coefficients()[0], 0.77323348, epsilon = 1e-7);
assert_relative_eq!(default_fit.coefficients()[1], 1.35480299, epsilon = 1e-7);
// Explicitly-constructed Cyclic is byte-identical to the default.
let explicit_res = ElasticNet::<f64>::new()
.with_alpha(0.3)
.with_l1_ratio(0.5)
.with_selection(CoordSelection::Cyclic)
.fit(&x, &y);
assert!(explicit_res.is_ok(), "explicit cyclic fit should succeed");
let explicit_cyclic = match explicit_res {
Ok(f) => f,
Err(_) => return,
};
assert_eq!(
default_fit.coefficients(),
explicit_cyclic.coefficients(),
"explicit Cyclic must be byte-identical to the default"
);
assert_eq!(default_fit.intercept(), explicit_cyclic.intercept());
}
// HONEST CAVEAT: exact bit-match to sklearn's `selection='random'` is
// numpy-MT19937-RNG-blocked (Rust `StdRng` != numpy MT19937), so the random
// path below verifies convergence-to-the-unique-optimum, NOT bitwise sklearn
// parity. The cyclic default IS bit-exact to sklearn (test above).
#[test]
fn enet_selection_random_converges_to_optimum() {
// Live sklearn 1.5.2 oracle (R-CHAR-3):
// ElasticNet(alpha=0.3, l1_ratio=0.5, selection='random',
// random_state=0)
// -> coef_ [0.77289352, 1.35505598] (same unique optimum,
// ~3e-4 from cyclic [0.77323348, 1.35480299] due to
// stopping-within-tol; NOT bit-identical to cyclic).
let (x, y) = selection_fixture();
let fit_res = ElasticNet::<f64>::new()
.with_alpha(0.3)
.with_l1_ratio(0.5)
.with_selection(CoordSelection::Random)
.with_random_state(0)
.fit(&x, &y);
assert!(fit_res.is_ok(), "random-selection fit should succeed");
let fitted = match fit_res {
Ok(f) => f,
Err(_) => return,
};
let coef = fitted.coefficients();
// Every coefficient finite.
for &c in coef.iter() {
assert!(c.is_finite(), "coefficient {c} must be finite");
}
// Converges to the unique cyclic optimum within tol.
let cyclic = [0.77323348_f64, 1.35480299_f64];
assert!(
(coef[0] - cyclic[0]).abs() < 1e-2,
"coef[0]={} should be within 1e-2 of cyclic {}",
coef[0],
cyclic[0]
);
assert!(
(coef[1] - cyclic[1]).abs() < 1e-2,
"coef[1]={} should be within 1e-2 of cyclic {}",
coef[1],
cyclic[1]
);
// Support set matches: both coefficients strictly positive.
assert!(coef[0] > 0.0, "coef[0] should be in the support");
assert!(coef[1] > 0.0, "coef[1] should be in the support");
}
// ---- precompute / Gram path (REQ-11) ----
#[test]
fn enet_precompute_matches_sklearn() -> Result<(), FerroError> {
// REQ-11: Gram (precompute=True) coordinate-descent path.
// Live sklearn 1.5.2 oracle (R-CHAR-3):
// X=[[1,2],[2,1],[3,4],[4,3],[5,5]], y=[3,2.5,7.1,6,11.2]
// ElasticNet(alpha=0.3, l1_ratio=0.5, precompute=True).fit(X,y)
// -> coef_=[0.7732334821, 1.3548029901], n_iter_=16
// (same optimum as precompute=False to ~1e-10).
let (x, y) = raw_dual_gap_fixture();
let fitted = ElasticNet::<f64>::new()
.with_alpha(0.3)
.with_l1_ratio(0.5)
.with_precompute(true)
.fit(&x, &y)?;
assert_relative_eq!(fitted.coefficients()[0], 0.7732334821, epsilon = 1e-7);
assert_relative_eq!(fitted.coefficients()[1], 1.3548029901, epsilon = 1e-7);
assert_eq!(fitted.n_iter(), 16, "n_iter_ must match sklearn's 16");
Ok(())
}
#[test]
fn enet_precompute_default_false_unchanged() -> Result<(), FerroError> {
// Default `precompute` is `false`; the default fit must be byte-identical
// to an explicitly-direct (precompute=false) fit (no perturbation).
assert!(
!ElasticNet::<f64>::new().precompute,
"default precompute is false"
);
let (x, y) = raw_dual_gap_fixture();
let default_fit = ElasticNet::<f64>::new()
.with_alpha(0.3)
.with_l1_ratio(0.5)
.fit(&x, &y)?;
let explicit_direct = ElasticNet::<f64>::new()
.with_alpha(0.3)
.with_l1_ratio(0.5)
.with_precompute(false)
.fit(&x, &y)?;
assert_eq!(
default_fit.coefficients(),
explicit_direct.coefficients(),
"explicit precompute=false must be byte-identical to the default"
);
assert_eq!(default_fit.intercept(), explicit_direct.intercept());
Ok(())
}
#[test]
fn enet_precompute_equals_direct() -> Result<(), FerroError> {
// The Gram path reaches the SAME unique optimum as the direct path,
// via different (reassociated) arithmetic — coef within 1e-6.
let (x, y) = raw_dual_gap_fixture();
let direct = ElasticNet::<f64>::new()
.with_alpha(0.3)
.with_l1_ratio(0.5)
.with_precompute(false)
.fit(&x, &y)?;
let gram = ElasticNet::<f64>::new()
.with_alpha(0.3)
.with_l1_ratio(0.5)
.with_precompute(true)
.fit(&x, &y)?;
assert_relative_eq!(
gram.coefficients()[0],
direct.coefficients()[0],
epsilon = 1e-6
);
assert_relative_eq!(
gram.coefficients()[1],
direct.coefficients()[1],
epsilon = 1e-6
);
Ok(())
}
// ---- warm_start (REQ-9) ----
/// Oracle fixture for the warm-start tests (R-CHAR-3, live sklearn 1.5.2):
/// `X = [[1,2],[2,1],[3,4],[4,3],[5,5]]`, `y = [3,2.5,7.1,6,11.2]`,
/// `alpha=0.5`, `l1_ratio=0.5`.
fn warm_start_fixture() -> (Array2<f64>, Array1<f64>) {
let x: Array2<f64> = array![[1.0, 2.0], [2.0, 1.0], [3.0, 4.0], [4.0, 3.0], [5.0, 5.0],];
let y: Array1<f64> = array![3.0, 2.5, 7.1, 6.0, 11.2];
(x, y)
}
#[test]
fn enet_warm_start_from_converged_matches_sklearn() -> Result<(), FerroError> {
// REQ-9: warm-start coordinate descent from the prior converged coef.
// Live sklearn 1.5.2 oracle (R-CHAR-3):
// ElasticNet(alpha=0.5, l1_ratio=0.5).fit(X,y)
// -> coef_=[0.7643620892, 1.2564536255], n_iter_=14 (cold)
// refit ElasticNet(alpha=0.5, l1_ratio=0.5, warm_start=True) from the
// converged coef_
// -> coef_=[0.7642996441, 1.2564980309], n_iter_=1 (warm)
let (x, y) = warm_start_fixture();
let cold = ElasticNet::<f64>::new()
.with_alpha(0.5)
.with_l1_ratio(0.5)
.fit(&x, &y)?;
assert_relative_eq!(cold.coefficients()[0], 0.7643620892, epsilon = 1e-6);
assert_relative_eq!(cold.coefficients()[1], 1.2564536255, epsilon = 1e-6);
assert_eq!(cold.n_iter(), 14, "cold n_iter_ must match sklearn's 14");
let warm = ElasticNet::<f64>::new()
.with_alpha(0.5)
.with_l1_ratio(0.5)
.with_warm_start(true)
.with_coef_init(cold.coefficients().to_owned())
.fit(&x, &y)?;
assert_relative_eq!(warm.coefficients()[0], 0.7642996441, epsilon = 1e-6);
assert_relative_eq!(warm.coefficients()[1], 1.2564980309, epsilon = 1e-6);
assert_eq!(warm.n_iter(), 1, "warm n_iter_ must match sklearn's 1");
Ok(())
}
#[test]
fn enet_warm_start_default_unchanged() -> Result<(), FerroError> {
// Default `warm_start=false`/`coef_init=None`; the default fit must be
// byte-identical to before (the cold zeros-init path is untouched).
assert!(
!ElasticNet::<f64>::new().warm_start,
"default warm_start is false"
);
assert!(
ElasticNet::<f64>::new().coef_init.is_none(),
"default coef_init is None"
);
let (x, y) = warm_start_fixture();
let default_fit = ElasticNet::<f64>::new()
.with_alpha(0.5)
.with_l1_ratio(0.5)
.fit(&x, &y)?;
let explicit_cold = ElasticNet::<f64>::new()
.with_alpha(0.5)
.with_l1_ratio(0.5)
.with_warm_start(false)
.fit(&x, &y)?;
// Bit-identical: same coordinate-descent start point (zeros).
assert_eq!(
default_fit.coefficients()[0].to_bits(),
explicit_cold.coefficients()[0].to_bits()
);
assert_eq!(
default_fit.coefficients()[1].to_bits(),
explicit_cold.coefficients()[1].to_bits()
);
assert_eq!(
default_fit.intercept().to_bits(),
explicit_cold.intercept().to_bits()
);
assert_eq!(default_fit.n_iter(), explicit_cold.n_iter());
Ok(())
}
#[test]
fn enet_warm_start_none_coef_init_equals_cold() -> Result<(), FerroError> {
// `warm_start=true` but no `coef_init` falls back to the zeros init,
// byte-identical to a plain cold fit (warm_start gates only whether
// `coef_init`, when present, is used).
let (x, y) = warm_start_fixture();
let cold = ElasticNet::<f64>::new()
.with_alpha(0.5)
.with_l1_ratio(0.5)
.fit(&x, &y)?;
let warm_no_init = ElasticNet::<f64>::new()
.with_alpha(0.5)
.with_l1_ratio(0.5)
.with_warm_start(true)
.fit(&x, &y)?;
assert_eq!(
cold.coefficients()[0].to_bits(),
warm_no_init.coefficients()[0].to_bits()
);
assert_eq!(
cold.coefficients()[1].to_bits(),
warm_no_init.coefficients()[1].to_bits()
);
assert_eq!(cold.n_iter(), warm_no_init.n_iter());
Ok(())
}
#[test]
fn enet_warm_start_coef_init_wrong_len_errors() {
// `coef_init` length (1) != n_features (2) must raise ShapeMismatch.
let (x, y) = warm_start_fixture();
let result = ElasticNet::<f64>::new()
.with_alpha(0.5)
.with_l1_ratio(0.5)
.with_warm_start(true)
.with_coef_init(array![0.0])
.fit(&x, &y);
assert!(
matches!(result, Err(FerroError::ShapeMismatch { .. })),
"wrong-length coef_init must return ShapeMismatch, got {result:?}"
);
}
}