ferrolearn_linear/

glm.rs

//! Generalized Linear Models (GLM).
//!
//! This module provides IRLS-based GLM regressors for count and positive
//! continuous data:
//!
//! - **[`GLMRegressor`]** — Generic GLM with selectable [`GLMFamily`]
//! - **[`PoissonRegressor`]** — Convenience wrapper with Poisson family
//! - **[`GammaRegressor`]** — Convenience wrapper with Gamma family
//! - **[`TweedieRegressor`]** — Convenience wrapper with Tweedie family
//!
//! All models use Iteratively Reweighted Least Squares (IRLS) and L2
//! regularization. The link function is fixed to **log** for Poisson and Gamma
//! (their sklearn losses are log-link only); [`TweedieRegressor`] selects its
//! [`Link`] via a `link` configuration (`auto`/`identity`/`log`), matching
//! `sklearn/linear_model/_glm/glm.py:889-903`.
//!
//! # Examples
//!
//! ```
//! use ferrolearn_linear::PoissonRegressor;
//! use ferrolearn_core::{Fit, Predict};
//! use ndarray::{array, Array1, Array2};
//!
//! let x = Array2::from_shape_vec((4, 1), vec![1.0, 2.0, 3.0, 4.0]).unwrap();
//! let y = array![2.0, 5.0, 10.0, 20.0];
//!
//! let model = PoissonRegressor::<f64>::new().with_alpha(0.0);
//! let fitted = model.fit(&x, &y).unwrap();
//! let preds = fitted.predict(&x).unwrap();
//! assert_eq!(preds.len(), 4);
//! ```
//!
//! ## REQ status (per `.design/linear/glm.md`, mirrors `sklearn/linear_model/_glm/glm.py` @ 1.5.2, commit 156ef14)
//!
//! Binary classification (R-DEFER-2): SHIPPED = impl + tests + green oracle
//! verification; NOT-STARTED = concrete open blocker referenced by `#`-number.
//! The public estimator types re-exported at the crate root are the consumer
//! surface (R-DEFER-1; no `ferrolearn-python` GLM binding yet).
//!
//! | REQ | Status | Evidence |
//! |---|---|---|
//! | REQ-4 (penalized objective: mean half-deviance + ½·alpha, intercept-free) | SHIPPED | `fn weighted_ridge_solve` adds the L2 penalty `weight_sum * alpha` to feature columns only, skipping the intercept column (`intercept_col`), matching sklearn's mean-deviance objective + unpenalized intercept (`glm.py:229-258`: `obj = average(½·deviance) + ½·alpha·‖coef‖²`, `l2_reg_strength = self.alpha`). Oracle parity tests `glm_poisson_intercept_unpenalized` (alpha=1e6 → `intercept_ = log(mean y)`, coef → 0) and `glm_poisson_penalty_scaling` (alpha=1.0 → `coef_=[0.34151720,0.18859745]`, `intercept_=-0.37680132`) green in `tests/divergence_glm_fit.rs`. |
//! | REQ-1 (Poisson family + log link) | SHIPPED | #548. `fn fit_glm_irls` with `GLMFamily::Poisson` (`variance => mu`, log-link Fisher scoring) under the REQ-4 mean-deviance/unpenalized-intercept objective fits `PoissonRegressor` to sklearn's `PoissonRegressor` (`HalfPoissonLoss`, log link, `glm.py:589-590`) at BOTH alpha=0 and alpha>0. Consumer: `PoissonRegressor::fit` (crate-root export). Oracle parity tests `glm_poisson_alpha_half_parity` (alpha=0.5 → live coef `[0.38388476754733647,0.2024000617918683]`, int `-0.519356533563308`), `glm_solver_param_invariant`, `glm_poisson_sample_weight`, `glm_poisson_penalty_scaling` green in `tests/divergence_glm_fit.rs` (the alpha=0 MLE matches to <1e-9, the module-header note). |
//! | REQ-2 (Gamma family + log link) | SHIPPED | #549. `GLMFamily::Gamma` (`variance => mu²`) drives `w = mu²/V(mu)` log-link IRLS; `GammaRegressor` matches sklearn's `GammaRegressor` (`HalfGammaLoss`, log link, y-domain `0 < y`, `glm.py:721-722`) at alpha=0 and alpha>0. The `y == 0` rejection (`HalfGammaLoss` open at 0) is enforced under REQ-14 (`YDomain::Positive`). Consumer: `GammaRegressor::fit` (crate-root export). Oracle tests `glm_gamma_alpha_half_parity` (alpha=0.5 → live coef `[0.24773782526507374,0.11636425618936652]`, int `0.3599464049766692`), `glm_gamma_sample_weight`, `glm_gamma_rejects_zero_y` green. |
//! | REQ-3 (Tweedie family + power) | SHIPPED | #550. `GLMFamily::Tweedie(p)` (`variance => mu^p`) with the link resolved from `power`/`link` (REQ-8); `TweedieRegressor(power=p)` matches sklearn's `TweedieRegressor` for the log-link powers `p>0` AND the identity-link `p<=0` (`HalfTweedieLoss`/`HalfTweedieLossIdentity`, `glm.py:889-903`) at alpha=0 and alpha>0. Verified live against the oracle for `p ∈ {0,1,1.5,2,3}` to <1e-8. Consumer: `TweedieRegressor::fit` (crate-root export). Oracle tests `glm_tweedie_alpha_half_parity` (`power=1.5, alpha=0.5` → live coef `[0.25606046404981164,0.11657692670900446]`, int `0.3563978246931595`), `glm_tweedie_power0_identity_link`, `glm_tweedie_power0_predict_identity_inverse`, `glm_tweedie_power2_rejects_zero_y` green. |
//! | REQ-5 (intercept init = link(weighted_mean(y))) | SHIPPED | #552. In `fn fit_glm_irls`, when `fit_intercept` AND NOT (warm_start with `coef_init`), the intercept entry is seeded at `coef[0] = link.link(weighted_mean(y))` (feature coefs stay 0) and `eta`/`mu` are recomputed from that seed, mirroring sklearn's `coef[-1] = link.link(np.average(y, weights=sample_weight))` (`glm.py:251-256`); `weighted_mean(y) = Σ(sᵢ·yᵢ)/Σ(sᵢ)`. warm_start with an explicit `coef_init` (REQ-11) takes precedence (the warm seed overrides the init), and a non-finite seed (e.g. `log(0)` for all-zero Poisson `y`) falls back to the previous cold start (intercept 0) with NO panic/NaN (R-CODE-2). The penalized GLM objective is convex, so the converged `coef_`/`intercept_` are init-invariant — all 22 pre-existing oracle tests stay byte-identical at convergence. Consumer: each estimator's `Fit::fit` (crate-root export). Oracle tests `glm_intercept_init_matches_sklearn_first_iterate` (constant `y=7` → `max_iter=1` intercept = `log(7) = 1.9459101490553132`, == live sklearn's first iterate; feature coef 0), `glm_intercept_init_converged_optimum_unchanged` (alpha=0.5 optimum unchanged vs oracle), `glm_intercept_init_all_zero_y_no_nan` (non-finite-seed fallback, finite result) green in `tests/divergence_glm_fit.rs`. |
//! | REQ-13 (score = D², deviance-explained) | SHIPPED | #559. `#[must_use] pub fn score(&self, x, y) -> Result<F, FerroError>` on `FittedGLMRegressor` computes `D² = 1 − (deviance + constant)/(deviance_null + constant)` (`glm.py:365-438`): `μ = predict(x)`, the null model predicts the (unweighted) mean `ȳ` for every sample, and the per-family unit deviance comes from `GLMFamily::unit_deviance` (Poisson `2·(y·ln(y/μ) − y + μ)`, y=0→`2μ`; Gamma `2·(−ln(y/μ) + (y−μ)/μ)`; Tweedie p=0 `(y−μ)²`; general-p `2·(y^(2−p)/((1−p)(2−p)) − y·μ^(1−p)/(1−p) + μ^(2−p)/(2−p))`), verified term-for-term against `sklearn/_loss/loss.py` (`HalfPoissonLoss:728-742`, `HalfGammaLoss:754-773`, `HalfTweedieLoss:789-837`). `GLMFamily::constant_to_optimal_zero` restores sklearn's `+ constant` so the degenerate constant-`y` boundary matches the oracle. `score` re-validates the y-domain (`YDomain::for_power`), mirroring `glm.py:413-417`. Consumer: the crate-root-exported `FittedGLMRegressor::score` (a public method on the boundary fitted type). Oracle tests `glm_poisson_d2_score` (D²=0.7979479374534378), `glm_gamma_d2_score` (0.8987486959882107), `glm_tweedie_power0_d2_score` (0.9319946452476573, == R²), `glm_tweedie_d2_score` (0.9277805586816806), `glm_score_rejects_out_of_domain_y` green in `tests/divergence_glm_fit.rs`; all 14 pre-existing glm divergence tests stay green. |
//! | REQ-7 (predict = link.inverse) | SHIPPED | `fn predict` applies `self.link.inverse(eta)` (`Link::Log => exp`, `Link::Identity => eta`), mirroring `glm.py:362` (`y_pred = link.inverse(raw_prediction)`). Consumer: the crate-root-exported `FittedGLMRegressor::predict` used by every wrapper; oracle test `glm_tweedie_power0_predict_identity_inverse` (identity link → raw linear predictor `[0.4,6.3,12.2,18.1]`) green in `tests/divergence_glm_fit.rs`. |
//! | REQ-8 (Tweedie link='auto'/identity/log) | SHIPPED | `pub enum Link { Log, Identity }` + `pub enum LinkConfig { Auto, Log, Identity }` with `LinkConfig::resolve(power)`: Auto → identity for `power <= 0`, log otherwise (`glm.py:889-893`). `TweedieRegressor.link: LinkConfig` (default `Auto`) is resolved at fit time and threaded into `fit_glm_irls`'s link-parameterized IRLS (`w = dmu_deta^2/V(mu)`, `z = eta + (y-mu)/dmu_deta`) and the fitted struct. Consumer: `TweedieRegressor::fit` (crate-root export); oracle test `glm_tweedie_power0_identity_link` (`coef_=[5.9]`, `intercept_=-5.5`, OLS) green. Poisson/Gamma wire `Link::Log` explicitly. |
//! | REQ-10 (solver param: lbfgs/newton-cholesky) | SHIPPED | #556. **R-DEV-2 (API parity):** `pub enum Solver { Lbfgs, NewtonCholesky }` + a `pub solver: Solver` field (default `Solver::Lbfgs`) on `GLMRegressor`/`PoissonRegressor`/`GammaRegressor`/`TweedieRegressor`, plus `#[must_use] fn with_solver`, mirroring sklearn's validated `solver` constructor parameter `StrOptions({"lbfgs","newton-cholesky"})` default `"lbfgs"` (`glm.py:140-145, :155`); the two-variant enum mirrors the `StrOptions` constraint. **R-DEV-7 (implementation differs, observable contract preserved):** ferrolearn fits all GLMs via IRLS/Fisher-scoring (`fn fit_glm_irls`) regardless of `solver` — the penalized GLM objective is convex, so IRLS reaches the SAME minimizer as both sklearn solvers (verified live: `PoissonRegressor(alpha=0.5)` gives coef `[0.38388523,0.20239975]`, int `-0.51935749` for `lbfgs` AND `newton-cholesky`, identical to ~1e-9). Consumer: each estimator's `Fit::fit` (crate-root export) — the `solver` field is part of the boundary estimator ABI. Oracle test `glm_solver_param_invariant` (fits with `Solver::Lbfgs` and `Solver::NewtonCholesky`, both coef/intercept match the solver-invariant live sklearn 1.5.2 oracle to 1e-4) green in `tests/divergence_glm_fit.rs`; the 19 pre-existing glm divergence tests stay green. |
//! | REQ-9 (Tweedie default power=0.0) | SHIPPED | `TweedieRegressor::new` sets `power: 0.0` (sklearn default, `glm.py:867`). Consumer: `TweedieRegressor::default`/`new` (crate-root export); oracle test `glm_tweedie_default_power` (`new().power == 0.0`) green. |
//! | REQ-11 (warm_start) | SHIPPED | #557. **R-DEV-2 (API parity):** `pub warm_start: bool` (default `false`) + `#[must_use] fn with_warm_start` on `GLMRegressor`/`PoissonRegressor`/`GammaRegressor`/`TweedieRegressor`, mirroring sklearn's `warm_start` parameter (`"boolean"`, default `False`, `glm.py:146, :158, :576, :708, :874`). **R-DEV-7 (Rust analog — immutable-estimator design, observable contract preserved):** sklearn's `warm_start=True` reuses the stateful `self.coef_`/`self.intercept_` mutated across `fit` calls as the optimizer's start (`glm.py:243-254`); ferrolearn's estimators are immutable (`fit(&self, ...)` never mutates `self`, no `self.coef_` to reuse), so the warm-start point is supplied EXPLICITLY via `pub coef_init: Option<(Array1<F>, F)>` + `#[must_use] fn with_coef_init(coef, intercept)`. `fn fit_glm_irls` seeds the IRLS coefficient vector (and derived `eta`/`mu`) from `coef_init` when `warm_start && coef_init.is_some()` (validating `feature_coef.len() == n_features`, else `ShapeMismatch`); otherwise the cold start (`coef = 0`) is byte-for-byte preserved. The penalized GLM objective is convex, so the converged `coef_`/`intercept_` are warm-start-INVARIANT — the init only changes the starting point (and iteration count), never the optimum — so the warm fit matches the cold fit AND the sklearn oracle (`glm.py:244-256`). Consumer: each estimator's `Fit::fit` (crate-root export) — the `warm_start`/`coef_init` fields are part of the boundary estimator ABI. Oracle tests `glm_warm_start_observable_contract` (warm fit from a perturbed init == cold fit == live sklearn 1.5.2 oracle `coef_=[0.38388477,0.20240006]`, `intercept_=-0.51935653` to 1e-6/1e-4) and `glm_warm_start_init_used` (seeding the exact optimum with `max_iter=1` lands at the solution, a cold `max_iter=1` fit does not — proves the init is genuinely used) green in `tests/divergence_glm_fit.rs`; the 20 pre-existing glm divergence tests stay green (all cold-start, byte-identical). |
//! | REQ-12 (sample_weight) | SHIPPED | `fn fit_with_sample_weight` on `GLMRegressor`/`PoissonRegressor`/`GammaRegressor`/`TweedieRegressor` threads an `Array1<F>` `sample_weight` into `fn fit_glm_irls`, where the IRLS `W` diagonal becomes `s_i * w_irls,i` (`weights[i] = weights[i] * sample_weight[i]`) and the L2-penalty scale is `weight_sum = S = sum_i s_i` (`sample_weight.iter().fold(..)`), matching sklearn's `sample_weight`-averaged deviance objective normalized by `sum(sample_weight)` (`glm.py:229-242`; `_check_sample_weight`, `glm.py:208-211`). Consumer: each estimator's `Fit::fit` (crate-root export) delegates with an all-ones weight vector, so the unweighted path is byte-identical (`weight_sum = n_samples`). Oracle tests `glm_poisson_sample_weight` (coef `[0.35738828,0.19717462]`, int `-0.43719203`) and `glm_gamma_sample_weight` (coef `[0.23049054,0.11350454]`, int `0.41955357`) green in `tests/divergence_glm_fit.rs`; the 8 pre-existing unweighted oracle tests stay green. |
//! | REQ-15 (non-finite input rejected) | SHIPPED | The shared IRLS entry `fn fit_glm_irls` — which every estimator (`PoissonRegressor`/`GammaRegressor`/`TweedieRegressor`/`GLMRegressor`) routes through — rejects any NaN/+/-inf in X, y, or `sample_weight` BEFORE the y-domain check and the IRLS loop with `FerroError::InvalidParameter`, mirroring sklearn's `_validate_data(force_all_finite=True)` (`glm.py:189-196`) + `_check_sample_weight` (default `force_all_finite=True`, `glm.py:211`) → `ValueError("Input X contains NaN.")` / `"... contains infinity ..."`. Placed ONCE at the shared entry (R-DEFER-8 single instance). `.iter().any(|v| !v.is_finite())` catches both NaN and Inf; the finite path is byte-identical (the unweighted `Fit::fit` delegates an all-ones weight). Verified vs the live sklearn 1.5.2 oracle (R-CHAR-3): `PoissonRegressor`/`GammaRegressor`/`TweedieRegressor`(`.fit`) raise `ValueError` for NaN/+inf/-inf in X, NaN/inf in y, and NaN/inf in sample_weight (`tests/divergence_linear_nonfinite_batch3.rs::{poisson,gamma,tweedie}_*`). Non-test consumer: each estimator's crate-root-exported `Fit::fit`. (#2261) |
//! | REQ-14 (n_iter_ + per-family y-domain validation) | SHIPPED | #560. Per-family y-domain guard in `fn fit_glm_irls`: `YDomain::for_power(family.domain_power())` then `y.iter().any(|&yi| !y_domain.contains(yi))` → `FerroError::InvalidParameter{name:"y", reason:"Some value(s) of y are out of the valid range of the loss '<loss>'."}`, mirroring sklearn's `if not base_loss.in_y_true_range(y): raise ValueError(...)` (`glm.py:221-225`). The valid range is keyed on the family's Tweedie `power` (NOT the link — verified vs the live oracle that `HalfTweedieLoss(p).interval_y_true == HalfTweedieLossIdentity(p).interval_y_true`): `power <= 0` unconstrained (Normal), `0 < power < 2` → `y >= 0` (Poisson `power=1`), `power >= 2` → `y > 0` (Gamma `power=2`, open at 0). `FittedGLMRegressor` gains `n_iter: usize` (the IRLS iteration count captured in the convergence loop) with `#[must_use] pub fn n_iter(&self) -> usize` — sklearn's `n_iter_` is the lbfgs count (`glm.py:110-114, :283`); ferrolearn's is the IRLS count (solvers differ, both report iterations-to-convergence). Consumer: `FittedGLMRegressor::n_iter` accessor on the crate-root-exported fitted type. Oracle tests `glm_gamma_rejects_zero_y` (Gamma rejects `y==0`, accepts `y>0`), `glm_tweedie_power2_rejects_zero_y` (`power=2.0` rejects `y==0`; `power=1.5` accepts it), `glm_poisson_rejects_negative_y` (rejects `y<0`, accepts `y==0`), `glm_n_iter_exposed` (`1 <= n_iter() <= max_iter`) green in `tests/divergence_glm_fit.rs`; the 10 pre-existing glm divergence tests stay green (all their `y` are in-domain). |

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, ScalarOperand};
use num_traits::{Float, FromPrimitive};

// ---------------------------------------------------------------------------
// Link
// ---------------------------------------------------------------------------

/// The link function `g` of a Generalized Linear Model, mapping the mean `mu`
/// to the linear predictor `eta = g(mu)` (and back via the inverse link `h`,
/// `mu = h(eta)`).
///
/// Mirrors the link carried by sklearn's loss classes
/// (`sklearn/linear_model/_glm/glm.py:119-131`): `HalfPoissonLoss`,
/// `HalfGammaLoss` and `HalfTweedieLoss` use the **log** link
/// (`y_pred = exp(X @ coef + intercept)`); `HalfSquaredError` and
/// `HalfTweedieLossIdentity` use the **identity** link
/// (`y_pred = X @ coef + intercept`).
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum Link {
    /// Log link: `g(mu) = ln(mu)`, inverse `h(eta) = exp(eta)`.
    ///
    /// Used by Poisson, Gamma and Tweedie-with-`power > 0` losses.
    Log,
    /// Identity link: `g(mu) = mu`, inverse `h(eta) = eta`.
    ///
    /// Used by the Normal/least-squares loss and Tweedie-with-`power <= 0`.
    Identity,
}

impl Link {
    /// Inverse link `h(eta) = mu`: maps the linear predictor to the mean.
    ///
    /// - [`Link::Log`] → `exp(eta)`
    /// - [`Link::Identity`] → `eta`
    ///
    /// Mirrors `link.inverse(raw_prediction)` in `glm.py:362`.
    #[must_use]
    fn inverse<F: Float>(self, eta: F) -> F {
        match self {
            Link::Log => eta.exp(),
            Link::Identity => eta,
        }
    }

    /// Forward link `g(mu) = eta`: maps the mean to the linear predictor.
    ///
    /// - [`Link::Log`] → `ln(mu)`
    /// - [`Link::Identity`] → `mu`
    ///
    /// Mirrors `link.link(...)` used by sklearn to seed the intercept at
    /// `link.link(average(y))` (`glm.py:254-256`).
    #[must_use]
    fn link<F: Float>(self, mu: F) -> F {
        match self {
            Link::Log => mu.ln(),
            Link::Identity => mu,
        }
    }

    /// Link derivative of the mean w.r.t. the linear predictor, `dmu/deta`,
    /// used to form the IRLS working weight and response.
    ///
    /// - [`Link::Log`] (`mu = exp(eta)`) → `dmu/deta = mu`
    /// - [`Link::Identity`] (`mu = eta`) → `dmu/deta = 1`
    #[must_use]
    fn dmu_deta<F: Float>(self, mu: F) -> F {
        match self {
            Link::Log => mu,
            Link::Identity => F::one(),
        }
    }
}

/// Configuration of the GLM link function, resolved to a concrete [`Link`] at
/// fit time.
///
/// Mirrors sklearn's `TweedieRegressor(link={'auto','identity','log'})`
/// (`glm.py:861, :889-903`). `Auto` selects the link from the Tweedie `power`:
/// identity for `power <= 0` (Normal), log otherwise (Poisson/Gamma/etc.).
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum LinkConfig {
    /// Resolve the link from the Tweedie `power` at fit time
    /// (`power <= 0` → identity, `power > 0` → log). The default.
    Auto,
    /// Force the log link regardless of `power`.
    Log,
    /// Force the identity link regardless of `power`.
    Identity,
}

impl LinkConfig {
    /// Resolve to a concrete [`Link`] given the Tweedie `power`.
    ///
    /// Mirrors `TweedieRegressor._get_loss` (`glm.py:889-903`):
    /// - `Auto` → identity for `power <= 0`, log for `power > 0`
    /// - `Log` → log; `Identity` → identity.
    #[must_use]
    fn resolve(self, power: f64) -> Link {
        match self {
            LinkConfig::Auto => {
                if power <= 0.0 {
                    Link::Identity
                } else {
                    Link::Log
                }
            }
            LinkConfig::Log => Link::Log,
            LinkConfig::Identity => Link::Identity,
        }
    }
}

// ---------------------------------------------------------------------------
// Solver
// ---------------------------------------------------------------------------

/// The optimization algorithm requested for fitting a GLM, mirroring sklearn's
/// `solver` constructor parameter (`sklearn/linear_model/_glm/glm.py:140-145`,
/// `StrOptions({"lbfgs", "newton-cholesky"})`, default `"lbfgs"`).
///
/// **Implementation note (R-DEV-7 — Rust analog).** ferrolearn fits all GLMs via
/// IRLS / Fisher-scoring (`fn fit_glm_irls`); the `solver` parameter is accepted
/// for scikit-learn API parity (R-DEV-2) and selects the requested optimizer's
/// *contract*. The penalized GLM objective is convex, so IRLS converges to the
/// **same** minimizer as both sklearn's `lbfgs` (scipy L-BFGS-B) and
/// `newton-cholesky` (Newton-Raphson with an inner Cholesky solve) — all three
/// are descent methods on one convex objective. Therefore the observable
/// contract (`coef_` / `intercept_`) matches sklearn for **either** `solver`
/// value, and ferrolearn does not vary the numerical path between them
/// (verified live: `PoissonRegressor(alpha=0.5)` gives the same fitted
/// attributes to ~1e-9 for `lbfgs` and `newton-cholesky`).
///
/// The type system constrains valid values to the two sklearn options, mirroring
/// the role of sklearn's `StrOptions` parameter constraint.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum Solver {
    /// L-BFGS-B (sklearn's default `"lbfgs"`): a quasi-Newton optimizer on the
    /// penalized loss + gradient (`glm.py:263-284`). In ferrolearn the fit is
    /// performed by IRLS, which reaches the same convex optimum (R-DEV-7).
    Lbfgs,
    /// Newton-Cholesky (sklearn's `"newton-cholesky"`): Newton-Raphson steps
    /// with an inner Cholesky solve, equivalent in exact arithmetic to iterated
    /// reweighted least squares (`glm.py:72-78, :285-296`). In ferrolearn the
    /// fit is performed by IRLS, which reaches the same convex optimum (R-DEV-7).
    NewtonCholesky,
}

// ---------------------------------------------------------------------------
// GLMFamily
// ---------------------------------------------------------------------------

/// The distributional family for a Generalized Linear Model.
///
/// Determines the variance function V(mu):
/// - **Poisson**: V(mu) = mu
/// - **Gamma**: V(mu) = mu^2
/// - **Tweedie(p)**: V(mu) = mu^p
#[derive(Debug, Clone, Copy)]
pub enum GLMFamily {
    /// Poisson family — variance proportional to the mean.
    Poisson,
    /// Gamma family — variance proportional to the squared mean.
    Gamma,
    /// Tweedie family with power parameter `p`.
    ///
    /// - `p = 0` gives Normal (constant variance)
    /// - `p = 1` gives Poisson
    /// - `p = 2` gives Gamma
    /// - `1 < p < 2` gives compound Poisson-Gamma
    Tweedie(f64),
}

impl GLMFamily {
    /// Compute the variance function V(mu) for a given mean `mu`.
    fn variance<F: Float + FromPrimitive>(&self, mu: F) -> F {
        match self {
            GLMFamily::Poisson => mu,
            GLMFamily::Gamma => mu * mu,
            GLMFamily::Tweedie(p) => {
                let power = F::from(*p).unwrap_or_else(F::zero);
                mu.powf(power)
            }
        }
    }

    /// The Tweedie power that determines this family's target (`y`) domain.
    ///
    /// The sklearn EDM losses derive their `interval_y_true` from the Tweedie
    /// `power` alone (`HalfPoissonLoss` is `power = 1`, `HalfGammaLoss` is
    /// `power = 2`); the link function does NOT change the valid `y` range
    /// (verified against the live oracle: `HalfTweedieLoss(power).interval_y_true
    /// == HalfTweedieLossIdentity(power).interval_y_true`). See
    /// [`YDomain::for_power`].
    #[must_use]
    fn domain_power(&self) -> f64 {
        match self {
            GLMFamily::Poisson => 1.0,
            GLMFamily::Gamma => 2.0,
            GLMFamily::Tweedie(p) => *p,
        }
    }

    /// The per-sample **unit deviance** `dev(y, mu)` of this family, the
    /// quantity whose `sample_weight`-weighted sum forms the D² numerator and
    /// denominator (`fn score`).
    ///
    /// This equals `2 * (loss(y, mu) + constant_to_optimal_zero(y))` of the
    /// matching sklearn EDM loss — i.e. twice the half-deviance with the
    /// `raw_prediction`-independent constant restored so a perfect prediction
    /// has zero deviance (`sklearn/_loss/loss.py`: `HalfPoissonLoss` `:728-742`,
    /// `HalfGammaLoss` `:754-773`, `HalfTweedieLoss` `:789-837`). Verified
    /// term-for-term against the live sklearn 1.5.2 oracle.
    ///
    /// Per family (`p` = Tweedie power):
    /// - **Poisson** (`p = 1`): `2·(y·ln(y/mu) − y + mu)`, with the convention
    ///   `y·ln(y/mu) → 0` at `y == 0` (so `dev = 2·mu`).
    /// - **Gamma** (`p = 2`): `2·(−ln(y/mu) + (y − mu)/mu)`.
    /// - **Tweedie `p == 0`** (Normal): `(y − mu)²`.
    /// - **Tweedie general** `p ∉ {0, 1, 2}`:
    ///   `2·( max(y,0)^(2−p)/((1−p)(2−p)) − y·mu^(1−p)/(1−p) + mu^(2−p)/(2−p) )`,
    ///   matching `HalfTweedieLoss.loss` (`loss.py:792-794`); `max(y, 0)` guards
    ///   `y == 0` (the `y^(2−p)` term is then 0).
    #[must_use]
    fn unit_deviance<F: Float + FromPrimitive>(&self, y: F, mu: F) -> F {
        let two = F::from(2.0).unwrap_or_else(|| F::one() + F::one());
        match self {
            GLMFamily::Poisson => {
                // 2·(y·ln(y/mu) − y + mu); y·ln(y/mu) → 0 at y == 0.
                let y_term = if y > F::zero() {
                    y * (y / mu).ln()
                } else {
                    F::zero()
                };
                two * (y_term - y + mu)
            }
            GLMFamily::Gamma => {
                // 2·(−ln(y/mu) + (y − mu)/mu).
                two * (F::zero() - (y / mu).ln() + (y - mu) / mu)
            }
            GLMFamily::Tweedie(p) => self.tweedie_unit_deviance(*p, y, mu),
        }
    }

    /// Tweedie unit deviance for an arbitrary power `p`, dispatching the
    /// `p ∈ {1, 2}` special cases to Poisson/Gamma and `p == 0` to the Normal
    /// squared error, matching `HalfTweedieLoss` taking the `p → 0, 1, 2` limits
    /// (`loss.py:796-797`).
    #[must_use]
    fn tweedie_unit_deviance<F: Float + FromPrimitive>(&self, p: f64, y: F, mu: F) -> F {
        let two = F::from(2.0).unwrap_or_else(|| F::one() + F::one());
        if p == 0.0 {
            // Normal: (y − mu)².
            let d = y - mu;
            return d * d;
        }
        if p == 1.0 {
            return GLMFamily::Poisson.unit_deviance(y, mu);
        }
        if p == 2.0 {
            return GLMFamily::Gamma.unit_deviance(y, mu);
        }
        // General p: 2·( max(y,0)^(2−p)/((1−p)(2−p)) − y·mu^(1−p)/(1−p)
        //                 + mu^(2−p)/(2−p) ).
        let pf = F::from(p).unwrap_or_else(F::zero);
        let one = F::one();
        let one_minus_p = one - pf;
        let two_minus_p = two - pf;
        let y_pos = if y > F::zero() { y } else { F::zero() };
        let t1 = y_pos.powf(two_minus_p) / (one_minus_p * two_minus_p);
        let t2 = y * mu.powf(one_minus_p) / one_minus_p;
        let t3 = mu.powf(two_minus_p) / two_minus_p;
        two * (t1 - t2 + t3)
    }

    /// The `sample_weight`-independent constant `constant_to_optimal_zero(y)` of
    /// the matching sklearn EDM loss — the term sklearn DROPS from its
    /// half-loss so that a perfect prediction scores zero deviance/2 (it is
    /// restored when forming the unit deviance).
    ///
    /// Mirrors `sklearn/_loss/loss.py`: `HalfPoissonLoss.constant_to_optimal_zero`
    /// `:738-742` = `xlogy(y, y) − y` (`xlogy(0, 0) = 0`); `HalfGammaLoss` `:769-773`
    /// = `−ln(y) − 1`; `HalfSquaredError`/Tweedie-identity `:453-458` = `0`;
    /// `HalfTweedieLoss.constant_to_optimal_zero` `:819-837` dispatches `p ∈ {0,1,2}`
    /// and is `max(y,0)^(2−p)/((1−p)(2−p))` otherwise.
    ///
    /// `score` adds this constant to both the model and the null half-deviance so
    /// the D² expression `1 − (deviance + constant)/(deviance_null + constant)`
    /// is byte-for-byte sklearn's (`glm.py:419-438`); it only affects the result
    /// at the degenerate constant-`y` boundary (`deviance_null == 0`), where
    /// sklearn returns `0` for Poisson (`constant ≠ 0`) and `NaN` for
    /// Gamma/Normal (`constant == 0`).
    #[must_use]
    fn constant_to_optimal_zero<F: Float + FromPrimitive>(&self, y: F) -> F {
        match self {
            GLMFamily::Poisson => {
                // xlogy(y, y) − y; xlogy(0, 0) = 0.
                let xlogy = if y > F::zero() { y * y.ln() } else { F::zero() };
                xlogy - y
            }
            GLMFamily::Gamma => {
                // −ln(y) − 1.
                F::zero() - y.ln() - F::one()
            }
            GLMFamily::Tweedie(p) => self.tweedie_constant_to_optimal_zero(*p, y),
        }
    }

    /// Tweedie `constant_to_optimal_zero` for an arbitrary power `p`, dispatching
    /// `p ∈ {0, 1, 2}` to Normal/Poisson/Gamma (`loss.py:819-831`) and using
    /// `max(y,0)^(2−p)/((1−p)(2−p))` otherwise (`loss.py:832-837`).
    #[must_use]
    fn tweedie_constant_to_optimal_zero<F: Float + FromPrimitive>(&self, p: f64, y: F) -> F {
        if p == 0.0 {
            // HalfSquaredError: 0.
            return F::zero();
        }
        if p == 1.0 {
            return GLMFamily::Poisson.constant_to_optimal_zero(y);
        }
        if p == 2.0 {
            return GLMFamily::Gamma.constant_to_optimal_zero(y);
        }
        let two = F::from(2.0).unwrap_or_else(|| F::one() + F::one());
        let pf = F::from(p).unwrap_or_else(F::zero);
        let one = F::one();
        let one_minus_p = one - pf;
        let two_minus_p = two - pf;
        let y_pos = if y > F::zero() { y } else { F::zero() };
        y_pos.powf(two_minus_p) / (one_minus_p * two_minus_p)
    }
}

/// The valid target (`y`) domain of a GLM loss, derived from its Tweedie
/// `power`, mirroring sklearn's `BaseLoss.interval_y_true` /
/// `in_y_true_range(y)` (`sklearn/linear_model/_glm/glm.py:221-225`).
///
/// The interval depends on `power` only (NOT the link); verified against the
/// live sklearn 1.5.2 oracle:
///
/// | power range          | valid `y`     | sklearn loss / interval                            |
/// |----------------------|---------------|----------------------------------------------------|
/// | `power <= 0`         | any real `y`  | `HalfSquaredError`/Tweedie identity, `(-inf, inf)` |
/// | `0 < power < 2`      | `y >= 0`      | Poisson (`power = 1`): `[0, inf)`                  |
/// | `power >= 2`         | `y > 0`       | Gamma (`power = 2`): `(0, inf)`                    |
///
/// (`power < 0` is the identity-link case, unconstrained; `0 < power < 1` is not
/// a standard Tweedie EDM but sklearn's `HalfTweedieLoss(power).interval_y_true`
/// is `[0, inf)` there, so we treat it as `y >= 0`.)
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum YDomain {
    /// `y` may be any real number (`power <= 0`, Normal/identity).
    Unconstrained,
    /// `y >= 0` (closed at 0; `0 < power < 2`, e.g. Poisson `power = 1`).
    NonNegative,
    /// `y > 0` (open at 0; `power >= 2`, e.g. Gamma `power = 2`).
    Positive,
}

impl YDomain {
    /// Resolve the valid `y` domain from a Tweedie `power`.
    ///
    /// Mirrors the boundaries of `HalfTweedieLoss(power).interval_y_true`
    /// (live sklearn 1.5.2 oracle): `power <= 0` → unconstrained, `0 < power < 2`
    /// → `y >= 0`, `power >= 2` → `y > 0`.
    #[must_use]
    fn for_power(power: f64) -> Self {
        if power <= 0.0 {
            YDomain::Unconstrained
        } else if power < 2.0 {
            YDomain::NonNegative
        } else {
            YDomain::Positive
        }
    }

    /// Human-readable description of the loss whose domain this is, for the
    /// error message (mirrors sklearn's `loss.__class__.__name__`,
    /// `glm.py:224`).
    #[must_use]
    fn loss_name(self) -> &'static str {
        match self {
            YDomain::Unconstrained => "HalfSquaredError",
            YDomain::NonNegative => "HalfPoissonLoss",
            YDomain::Positive => "HalfGammaLoss",
        }
    }

    /// Whether a single target value `yi` is inside this domain.
    #[must_use]
    fn contains<F: Float>(self, yi: F) -> bool {
        match self {
            YDomain::Unconstrained => true,
            YDomain::NonNegative => yi >= F::zero(),
            YDomain::Positive => yi > F::zero(),
        }
    }
}

// ---------------------------------------------------------------------------
// GLMRegressor
// ---------------------------------------------------------------------------

/// Generalized Linear Model regressor.
///
/// Fitted via IRLS with a log link function. The [`GLMFamily`] controls
/// the assumed variance-mean relationship.
///
/// # Type Parameters
///
/// - `F`: The floating-point type (`f32` or `f64`).
#[derive(Debug, Clone)]
pub struct GLMRegressor<F> {
    /// Distributional family (Poisson, Gamma, or Tweedie).
    pub family: GLMFamily,
    /// L2 regularization strength.
    pub alpha: F,
    /// Maximum number of IRLS iterations.
    pub max_iter: usize,
    /// Convergence tolerance on the maximum coefficient change.
    pub tol: F,
    /// Whether to fit an intercept (bias) term.
    pub fit_intercept: bool,
    /// Optimization algorithm requested, mirroring sklearn's `solver` parameter
    /// (`glm.py:140-145`, default `"lbfgs"`).
    ///
    /// ferrolearn fits via IRLS / Fisher-scoring regardless of this value
    /// (R-DEV-7); the parameter is accepted for sklearn API parity (R-DEV-2) and
    /// the observable `coef_` / `intercept_` match sklearn for either value
    /// because IRLS reaches the same convex optimum as both `lbfgs` and
    /// `newton-cholesky`. See [`Solver`].
    pub solver: Solver,
    /// Whether to warm-start the optimizer from an explicit initial point,
    /// mirroring sklearn's `warm_start` constructor parameter (default `false`,
    /// `glm.py:146, :158, :244`).
    ///
    /// **R-DEV-2 (API parity):** the name and default match sklearn. **R-DEV-7
    /// (Rust analog):** sklearn reuses the stateful `self.coef_` / `self.intercept_`
    /// across `fit` calls (`glm.py:244-250`); ferrolearn's estimators are immutable
    /// (`fit(&self, ...)` never mutates `self`), so the warm-start point is supplied
    /// EXPLICITLY via [`GLMRegressor::with_coef_init`]. When `warm_start == true`
    /// and an init is set, the IRLS seeds from it; otherwise it cold-starts at
    /// `coef = 0`. Because the GLM objective is convex, the converged
    /// `coef_` / `intercept_` are warm-start-invariant — identical to the cold
    /// fit and to sklearn regardless of the init.
    pub warm_start: bool,
    /// Explicit warm-start initial point `(feature_coefficients, intercept)` —
    /// the ferrolearn analog of sklearn reusing `self.coef_` / `self.intercept_`
    /// across `fit` calls (R-DEV-7, `glm.py:244-250`).
    ///
    /// Only consulted when [`GLMRegressor::warm_start`] is `true`. The
    /// feature-coefficient vector length must equal the number of features in `X`
    /// (else [`FerroError::ShapeMismatch`] at fit time). Set via
    /// [`GLMRegressor::with_coef_init`]; typically a previous fit's
    /// `coefficients()` / `intercept()`.
    pub coef_init: Option<(Array1<F>, F)>,
}

impl<F: Float + FromPrimitive> GLMRegressor<F> {
    /// Create a new `GLMRegressor` with the given family.
    ///
    /// Defaults: `alpha = 1.0`, `max_iter = 100`, `tol = 1e-4`,
    /// `fit_intercept = true`, `solver = Solver::Lbfgs` (sklearn default).
    #[must_use]
    pub fn new(family: GLMFamily) -> Self {
        Self {
            family,
            alpha: F::one(),
            max_iter: 100,
            tol: F::from(1e-4).unwrap_or_else(F::epsilon),
            fit_intercept: true,
            solver: Solver::Lbfgs,
            warm_start: false,
            coef_init: None,
        }
    }

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

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

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

    /// Set 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 the optimization [`Solver`], mirroring sklearn's `solver` parameter
    /// (`glm.py:140-145`, default `"lbfgs"`).
    ///
    /// ferrolearn fits via IRLS regardless of the value (R-DEV-7); the parameter
    /// is accepted for sklearn API parity (R-DEV-2) and both values produce the
    /// same observable `coef_` / `intercept_` (IRLS reaches the same convex
    /// optimum as both `lbfgs` and `newton-cholesky`).
    #[must_use]
    pub fn with_solver(mut self, solver: Solver) -> Self {
        self.solver = solver;
        self
    }

    /// Enable or disable warm-starting, mirroring sklearn's `warm_start`
    /// parameter (default `false`, `glm.py:146, :158`).
    ///
    /// **R-DEV-2 / R-DEV-7.** sklearn reuses the stateful `self.coef_` /
    /// `self.intercept_` across `fit` calls (`glm.py:244-250`); ferrolearn's
    /// immutable estimators take the warm-start point EXPLICITLY via
    /// [`GLMRegressor::with_coef_init`]. When `warm_start` is `true` and an init
    /// is set, the IRLS seeds from it; otherwise it cold-starts at `coef = 0`. The
    /// convex GLM objective makes the converged `coef_` / `intercept_`
    /// warm-start-invariant.
    #[must_use]
    pub fn with_warm_start(mut self, warm_start: bool) -> Self {
        self.warm_start = warm_start;
        self
    }

    /// Set the explicit warm-start initial point `(feature_coefficients,
    /// intercept)` — the ferrolearn analog of sklearn reusing `self.coef_` /
    /// `self.intercept_` (R-DEV-7, `glm.py:244-250`).
    ///
    /// Only consulted when [`GLMRegressor::warm_start`] is `true`. `coef`'s length
    /// must equal the number of features in `X` (validated at fit time). Pass a
    /// previous fit's `coefficients()` / `intercept()` to resume from it.
    #[must_use]
    pub fn with_coef_init(mut self, coef: Array1<F>, intercept: F) -> Self {
        self.coef_init = Some((coef, intercept));
        self
    }
}

impl<F: Float + Send + Sync + ScalarOperand + FromPrimitive + 'static> GLMRegressor<F> {
    /// Fit the GLM via IRLS with per-sample weights `sample_weight`.
    ///
    /// Mirrors sklearn's `fit(X, y, sample_weight)` (`glm.py:170`,
    /// `:229-242`): the deviance term is a `sample_weight`-weighted average
    /// (normalized by `S = sum_i s_i`), so the IRLS `W` diagonal becomes
    /// `s_i * w_irls,i` and the L2 penalty scales with `S`. The working
    /// response `z_i` is unchanged per sample. Passing an all-ones weight
    /// vector reproduces [`Fit::fit`] exactly.
    ///
    /// # Errors
    ///
    /// - [`FerroError::ShapeMismatch`] — `sample_weight.len() != n_samples`,
    ///   or `y` length mismatch.
    /// - [`FerroError::InvalidParameter`] — a negative sample weight, negative
    ///   alpha, or (log link) negative `y`.
    /// - [`FerroError::InsufficientSamples`] — zero samples.
    pub fn fit_with_sample_weight(
        &self,
        x: &Array2<F>,
        y: &Array1<F>,
        sample_weight: &Array1<F>,
    ) -> Result<FittedGLMRegressor<F>, FerroError> {
        fit_glm_irls(
            x,
            y,
            sample_weight,
            &self.family,
            Link::Log,
            self.alpha,
            self.max_iter,
            self.tol,
            self.fit_intercept,
            self.warm_start,
            self.coef_init.as_ref(),
        )
    }
}

/// Fitted GLM regressor.
///
/// Stores the learned coefficients and intercept on the link scale, together
/// with the [`Link`] used at fit time. Predictions are
/// `link.inverse(X @ coef + intercept)` — `exp(...)` for [`Link::Log`], the raw
/// linear predictor for [`Link::Identity`] (`glm.py:362`).
#[derive(Debug, Clone)]
pub struct FittedGLMRegressor<F> {
    /// Learned coefficient vector on the link scale.
    coefficients: Array1<F>,
    /// Learned intercept on the link scale.
    intercept: F,
    /// Link function applied by `predict` (inverse link maps `eta` to `mu`).
    link: Link,
    /// Distributional family, retained so [`FittedGLMRegressor::score`] can
    /// compute this family's unit deviance for the D² (deviance-explained)
    /// score (`glm.py:365-438`).
    family: GLMFamily,
    /// Number of IRLS iterations actually run (until convergence or `max_iter`).
    ///
    /// Mirrors sklearn's fitted `n_iter_` attribute (`glm.py:110-114, :283`),
    /// the solver's iteration count. ferrolearn's solver is IRLS (not lbfgs), so
    /// this is the **IRLS** iteration count; both report iterations-to-convergence
    /// but the exact value is solver-dependent.
    n_iter: usize,
}

impl<F> FittedGLMRegressor<F> {
    /// Number of IRLS iterations run during the fit
    /// (until convergence or `max_iter`).
    ///
    /// Mirrors scikit-learn's fitted `n_iter_` attribute
    /// (`sklearn/linear_model/_glm/glm.py:110-114, :283`), which reports the
    /// solver's iteration count. ferrolearn's solver is **IRLS** (sklearn's
    /// default is lbfgs), so this is the IRLS iteration count, not the lbfgs
    /// one — both report iterations-to-convergence, but the exact value differs
    /// because the solvers differ. The value is in `1..=max_iter`.
    #[must_use]
    pub fn n_iter(&self) -> usize {
        self.n_iter
    }
}

impl<F: Float + Send + Sync + ScalarOperand + FromPrimitive + 'static> FittedGLMRegressor<F> {
    /// Compute **D²**, the fraction of deviance explained — the GLM
    /// generalization of R² (which is the special case for the Normal family).
    ///
    /// Mirrors `sklearn`'s `_GeneralizedLinearRegressor.score`
    /// (`sklearn/linear_model/_glm/glm.py:365-438`):
    ///
    /// ```text
    /// D² = 1 − (deviance + constant) / (deviance_null + constant)
    /// ```
    ///
    /// where `deviance` is the `sample_weight`-mean half-deviance of the fitted
    /// model's predictions `mu = predict(X)`, `deviance_null` is the same for the
    /// constant null model that predicts the weighted mean of `y` for every
    /// sample (`y_pred = link.inverse(link(mean_w(y))) = mean_w(y)`), and
    /// `constant` is the family's `constant_to_optimal_zero(y)` mean
    /// (`glm.py:419-438`). Because the dropped factor-of-2 and the `constant`
    /// cancel whenever `deviance_null ≠ 0`, the result equals the unit-deviance
    /// ratio `1 − Σ sᵢ·dev(yᵢ, muᵢ) / Σ sᵢ·dev(yᵢ, ȳ_w)`. The best possible
    /// score is `1.0`; it can be negative (an arbitrarily worse-than-null model).
    ///
    /// This is the unweighted path (`sample_weight = None`); the weighted mean
    /// `ȳ_w` reduces to the plain mean of `y`.
    ///
    /// **Degenerate (constant-`y`) boundary.** When `y` is constant the null
    /// deviance is `0`, so sklearn evaluates `1 − (deviance + constant)/0`. The
    /// returned value is sklearn's exact algebraic form, but the result there
    /// depends on the fitted model reproducing `mu == ȳ` to full precision: with
    /// sklearn's lbfgs solver `deviance + constant == deviance_null + constant`
    /// (both equal the float roundoff of `loss + constant`), giving `0` for
    /// Poisson and `NaN` for Gamma/Normal (`constant == 0`). ferrolearn's IRLS
    /// converges `mu` to `ȳ` only to solver tolerance (not bit-exactly), so on a
    /// constant-`y` Poisson input the ratio is `~1 ± ε` (≈ `0`) rather than
    /// exactly `0`; this is a fit-precision artifact of the degenerate input, not
    /// of the D² formula (non-degenerate inputs match the oracle to `< 1e-6`).
    ///
    /// `y` is re-validated against the family's target domain, mirroring
    /// sklearn's `if not base_loss.in_y_true_range(y): raise ValueError(...)` in
    /// `score` (`glm.py:413-417`).
    ///
    /// # Errors
    ///
    /// - [`FerroError::ShapeMismatch`] — `X` feature count mismatch (via
    ///   [`Predict::predict`]) or `y` length mismatch.
    /// - [`FerroError::InvalidParameter`] — a value of `y` is out of the family's
    ///   valid target range.
    #[must_use = "the computed D² score should be used"]
    pub fn score(&self, x: &Array2<F>, y: &Array1<F>) -> Result<F, FerroError> {
        let mu = self.predict(x)?;

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

        // Re-validate y against the family's target domain, mirroring sklearn's
        // `score`: `if not base_loss.in_y_true_range(y): raise ValueError(...)`
        // (`glm.py:413-417`). Same per-family domain as `fit` (keyed on power).
        let y_domain = YDomain::for_power(self.family.domain_power());
        if y.iter().any(|&yi| !y_domain.contains(yi)) {
            return Err(FerroError::InvalidParameter {
                name: "y".into(),
                reason: format!(
                    "Some value(s) of y are out of the valid range of the loss '{}'.",
                    y_domain.loss_name()
                ),
            });
        }

        // Weighted mean of y (unweighted => plain mean). The null model predicts
        // `ȳ` for every sample (`glm.py:431`: `link.link(average(y))` mapped back
        // through `link.inverse` in the null deviance is just `ȳ`).
        let n = F::from(y.len()).unwrap_or_else(F::one);
        let y_bar = y.iter().fold(F::zero(), |acc, &yi| acc + yi) / n;

        // sklearn forms `deviance + constant` and `deviance_null + constant`,
        // where `deviance = mean(loss)`, `deviance_null = mean(loss_null)` and
        // `constant = mean(constant_to_optimal_zero(y))` (`glm.py:419-437`).
        // Because per sample `loss + constant_to_optimal_zero = ½·unit_deviance`
        // (the `loss` drops exactly the constant the unit deviance restores), the
        // two terms sklearn divides are simply the mean half unit deviances:
        //   deviance + constant      = mean(½·unit_deviance(y, mu_model))
        //   deviance_null + constant = mean(½·unit_deviance(y, ȳ)).
        // We accumulate `½·unit_deviance` on each side directly — adding `constant`
        // again would double-count it. The factor ½ cancels in the ratio for any
        // nonzero denominator; we keep it so the constant-`y` boundary
        // (`deviance_null == 0`) reproduces sklearn's exact `½·unit_dev_model / 0`
        // → it must, however, also restore the dropped `constant` at that boundary
        // to distinguish Poisson (0) from Gamma/Normal (NaN). We therefore compute
        // sklearn's `(deviance + constant)` / `(deviance_null + constant)` form
        // with `deviance = mean(loss)`, `loss = ½·unit_deviance − constant`.
        let half = F::from(0.5).unwrap_or_else(|| F::one() / (F::one() + F::one()));
        let mut sum_loss_model = F::zero();
        let mut sum_loss_null = F::zero();
        let mut sum_const = F::zero();
        for (&yi, &mui) in y.iter().zip(mu.iter()) {
            let c = self.family.constant_to_optimal_zero(yi);
            sum_loss_model = sum_loss_model + (half * self.family.unit_deviance(yi, mui) - c);
            sum_loss_null = sum_loss_null + (half * self.family.unit_deviance(yi, y_bar) - c);
            sum_const = sum_const + c;
        }
        let deviance = sum_loss_model / n;
        let deviance_null = sum_loss_null / n;
        let constant = sum_const / n;

        // `1 − (deviance + constant) / (deviance_null + constant)` (glm.py:438).
        // Away from the constant-`y` boundary this equals the unit-deviance ratio
        // `1 − Σ dev(y, mu) / Σ dev(y, ȳ)`; at `deviance_null == 0` the restored
        // `constant` reproduces sklearn's family-dependent result — 0 for Poisson
        // (`constant ≠ 0`), NaN for Gamma/Normal (`constant == 0`) — verified
        // against the live sklearn 1.5.2 oracle.
        Ok(F::one() - (deviance + constant) / (deviance_null + constant))
    }
}

// ---------------------------------------------------------------------------
// Convenience wrappers
// ---------------------------------------------------------------------------

/// Poisson regressor — GLM with Poisson family and log link.
///
/// Suitable for modelling count data (y >= 0, integer-valued).
///
/// # Type Parameters
///
/// - `F`: The floating-point type (`f32` or `f64`).
#[derive(Debug, Clone)]
pub struct PoissonRegressor<F> {
    /// L2 regularization strength.
    pub alpha: F,
    /// Maximum number of IRLS iterations.
    pub max_iter: usize,
    /// Convergence tolerance.
    pub tol: F,
    /// Whether to fit an intercept.
    pub fit_intercept: bool,
    /// Optimization algorithm requested, mirroring sklearn's `solver` parameter
    /// (`glm.py:140-145`, default `"lbfgs"`).
    ///
    /// ferrolearn fits via IRLS regardless of this value (R-DEV-7); the parameter
    /// is accepted for sklearn API parity (R-DEV-2) and the observable
    /// `coef_` / `intercept_` match sklearn for either value. See [`Solver`].
    pub solver: Solver,
    /// Whether to warm-start from an explicit initial point, mirroring sklearn's
    /// `warm_start` parameter (default `false`, `glm.py:146, :576`). See
    /// [`GLMRegressor::warm_start`] for the R-DEV-2 / R-DEV-7 rationale.
    pub warm_start: bool,
    /// Explicit warm-start initial point `(feature_coefficients, intercept)` — the
    /// ferrolearn analog of sklearn reusing `self.coef_` / `self.intercept_`
    /// (R-DEV-7, `glm.py:244-250`). Consulted only when `warm_start` is `true`. Set
    /// via [`PoissonRegressor::with_coef_init`].
    pub coef_init: Option<(Array1<F>, F)>,
}

impl<F: Float + FromPrimitive> PoissonRegressor<F> {
    /// Create a new `PoissonRegressor` with default settings.
    ///
    /// Defaults: `alpha = 1.0`, `max_iter = 100`, `tol = 1e-4`,
    /// `fit_intercept = true`, `solver = Solver::Lbfgs` (sklearn default).
    #[must_use]
    pub fn new() -> Self {
        Self {
            alpha: F::one(),
            max_iter: 100,
            tol: F::from(1e-4).unwrap_or_else(F::epsilon),
            fit_intercept: true,
            solver: Solver::Lbfgs,
            warm_start: false,
            coef_init: None,
        }
    }

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

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

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

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

    /// Set the optimization [`Solver`], mirroring sklearn's `solver` parameter
    /// (`glm.py:140-145`, default `"lbfgs"`).
    ///
    /// ferrolearn fits via IRLS regardless of the value (R-DEV-7); both values
    /// produce the same observable `coef_` / `intercept_` (sklearn API parity,
    /// R-DEV-2).
    #[must_use]
    pub fn with_solver(mut self, solver: Solver) -> Self {
        self.solver = solver;
        self
    }

    /// Enable or disable warm-starting, mirroring sklearn's `warm_start`
    /// parameter (default `false`, `glm.py:146, :576`). See
    /// [`GLMRegressor::with_warm_start`] for the R-DEV-2 / R-DEV-7 rationale; the
    /// warm-start point is supplied explicitly via
    /// [`PoissonRegressor::with_coef_init`].
    #[must_use]
    pub fn with_warm_start(mut self, warm_start: bool) -> Self {
        self.warm_start = warm_start;
        self
    }

    /// Set the explicit warm-start initial point `(feature_coefficients,
    /// intercept)` — the ferrolearn analog of sklearn reusing `self.coef_` /
    /// `self.intercept_` (R-DEV-7, `glm.py:244-250`). Only consulted when
    /// `warm_start` is `true`; `coef`'s length must equal the number of features.
    #[must_use]
    pub fn with_coef_init(mut self, coef: Array1<F>, intercept: F) -> Self {
        self.coef_init = Some((coef, intercept));
        self
    }
}

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

impl<F: Float + Send + Sync + ScalarOperand + FromPrimitive + 'static> PoissonRegressor<F> {
    /// Fit the Poisson GLM via IRLS with per-sample weights `sample_weight`.
    ///
    /// Mirrors sklearn's `PoissonRegressor.fit(X, y, sample_weight)`
    /// (`glm.py:170`, `:229-242`). See
    /// [`GLMRegressor::fit_with_sample_weight`] for the weighting semantics;
    /// an all-ones weight vector reproduces [`Fit::fit`] exactly.
    ///
    /// # Errors
    ///
    /// See [`GLMRegressor::fit_with_sample_weight`].
    pub fn fit_with_sample_weight(
        &self,
        x: &Array2<F>,
        y: &Array1<F>,
        sample_weight: &Array1<F>,
    ) -> Result<FittedGLMRegressor<F>, FerroError> {
        fit_glm_irls(
            x,
            y,
            sample_weight,
            &GLMFamily::Poisson,
            Link::Log,
            self.alpha,
            self.max_iter,
            self.tol,
            self.fit_intercept,
            self.warm_start,
            self.coef_init.as_ref(),
        )
    }
}

/// Gamma regressor — GLM with Gamma family and log link.
///
/// Suitable for modelling positive continuous data (y > 0).
///
/// # Type Parameters
///
/// - `F`: The floating-point type (`f32` or `f64`).
#[derive(Debug, Clone)]
pub struct GammaRegressor<F> {
    /// L2 regularization strength.
    pub alpha: F,
    /// Maximum number of IRLS iterations.
    pub max_iter: usize,
    /// Convergence tolerance.
    pub tol: F,
    /// Whether to fit an intercept.
    pub fit_intercept: bool,
    /// Optimization algorithm requested, mirroring sklearn's `solver` parameter
    /// (`glm.py:140-145`, default `"lbfgs"`).
    ///
    /// ferrolearn fits via IRLS regardless of this value (R-DEV-7); the parameter
    /// is accepted for sklearn API parity (R-DEV-2) and the observable
    /// `coef_` / `intercept_` match sklearn for either value. See [`Solver`].
    pub solver: Solver,
    /// Whether to warm-start from an explicit initial point, mirroring sklearn's
    /// `warm_start` parameter (default `false`, `glm.py:146, :708`). See
    /// [`GLMRegressor::warm_start`] for the R-DEV-2 / R-DEV-7 rationale.
    pub warm_start: bool,
    /// Explicit warm-start initial point `(feature_coefficients, intercept)` — the
    /// ferrolearn analog of sklearn reusing `self.coef_` / `self.intercept_`
    /// (R-DEV-7, `glm.py:244-250`). Consulted only when `warm_start` is `true`. Set
    /// via [`GammaRegressor::with_coef_init`].
    pub coef_init: Option<(Array1<F>, F)>,
}

impl<F: Float + FromPrimitive> GammaRegressor<F> {
    /// Create a new `GammaRegressor` with default settings.
    ///
    /// Defaults: `alpha = 1.0`, `max_iter = 100`, `tol = 1e-4`,
    /// `fit_intercept = true`, `solver = Solver::Lbfgs` (sklearn default).
    #[must_use]
    pub fn new() -> Self {
        Self {
            alpha: F::one(),
            max_iter: 100,
            tol: F::from(1e-4).unwrap_or_else(F::epsilon),
            fit_intercept: true,
            solver: Solver::Lbfgs,
            warm_start: false,
            coef_init: None,
        }
    }

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

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

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

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

    /// Set the optimization [`Solver`], mirroring sklearn's `solver` parameter
    /// (`glm.py:140-145`, default `"lbfgs"`).
    ///
    /// ferrolearn fits via IRLS regardless of the value (R-DEV-7); both values
    /// produce the same observable `coef_` / `intercept_` (sklearn API parity,
    /// R-DEV-2).
    #[must_use]
    pub fn with_solver(mut self, solver: Solver) -> Self {
        self.solver = solver;
        self
    }

    /// Enable or disable warm-starting, mirroring sklearn's `warm_start`
    /// parameter (default `false`, `glm.py:146, :708`). See
    /// [`GLMRegressor::with_warm_start`] for the R-DEV-2 / R-DEV-7 rationale; the
    /// warm-start point is supplied explicitly via
    /// [`GammaRegressor::with_coef_init`].
    #[must_use]
    pub fn with_warm_start(mut self, warm_start: bool) -> Self {
        self.warm_start = warm_start;
        self
    }

    /// Set the explicit warm-start initial point `(feature_coefficients,
    /// intercept)` — the ferrolearn analog of sklearn reusing `self.coef_` /
    /// `self.intercept_` (R-DEV-7, `glm.py:244-250`). Only consulted when
    /// `warm_start` is `true`; `coef`'s length must equal the number of features.
    #[must_use]
    pub fn with_coef_init(mut self, coef: Array1<F>, intercept: F) -> Self {
        self.coef_init = Some((coef, intercept));
        self
    }
}

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

impl<F: Float + Send + Sync + ScalarOperand + FromPrimitive + 'static> GammaRegressor<F> {
    /// Fit the Gamma GLM via IRLS with per-sample weights `sample_weight`.
    ///
    /// Mirrors sklearn's `GammaRegressor.fit(X, y, sample_weight)`
    /// (`glm.py:170`, `:229-242`). See
    /// [`GLMRegressor::fit_with_sample_weight`] for the weighting semantics;
    /// an all-ones weight vector reproduces [`Fit::fit`] exactly.
    ///
    /// # Errors
    ///
    /// See [`GLMRegressor::fit_with_sample_weight`].
    pub fn fit_with_sample_weight(
        &self,
        x: &Array2<F>,
        y: &Array1<F>,
        sample_weight: &Array1<F>,
    ) -> Result<FittedGLMRegressor<F>, FerroError> {
        fit_glm_irls(
            x,
            y,
            sample_weight,
            &GLMFamily::Gamma,
            Link::Log,
            self.alpha,
            self.max_iter,
            self.tol,
            self.fit_intercept,
            self.warm_start,
            self.coef_init.as_ref(),
        )
    }
}

/// Tweedie regressor — GLM with Tweedie family and log link.
///
/// The `power` parameter controls the variance-mean relationship:
/// V(mu) = mu^power.
///
/// # Type Parameters
///
/// - `F`: The floating-point type (`f32` or `f64`).
#[derive(Debug, Clone)]
pub struct TweedieRegressor<F> {
    /// Tweedie power parameter.
    pub power: f64,
    /// Link-function configuration (`Auto`/`Log`/`Identity`).
    ///
    /// `Auto` (the default) resolves to the identity link for `power <= 0`
    /// (Normal) and the log link for `power > 0`, matching sklearn's
    /// `link='auto'` (`glm.py:889-893`).
    pub link: LinkConfig,
    /// L2 regularization strength.
    pub alpha: F,
    /// Maximum number of IRLS iterations.
    pub max_iter: usize,
    /// Convergence tolerance.
    pub tol: F,
    /// Whether to fit an intercept.
    pub fit_intercept: bool,
    /// Optimization algorithm requested, mirroring sklearn's `solver` parameter
    /// (`glm.py:140-145`, default `"lbfgs"`).
    ///
    /// ferrolearn fits via IRLS regardless of this value (R-DEV-7); the parameter
    /// is accepted for sklearn API parity (R-DEV-2) and the observable
    /// `coef_` / `intercept_` match sklearn for either value. See [`Solver`].
    pub solver: Solver,
    /// Whether to warm-start from an explicit initial point, mirroring sklearn's
    /// `warm_start` parameter (default `false`, `glm.py:146, :874`). See
    /// [`GLMRegressor::warm_start`] for the R-DEV-2 / R-DEV-7 rationale.
    pub warm_start: bool,
    /// Explicit warm-start initial point `(feature_coefficients, intercept)` — the
    /// ferrolearn analog of sklearn reusing `self.coef_` / `self.intercept_`
    /// (R-DEV-7, `glm.py:244-250`). Consulted only when `warm_start` is `true`. Set
    /// via [`TweedieRegressor::with_coef_init`].
    pub coef_init: Option<(Array1<F>, F)>,
}

impl<F: Float + FromPrimitive> TweedieRegressor<F> {
    /// Create a new `TweedieRegressor` with default settings.
    ///
    /// Defaults match sklearn's `TweedieRegressor.__init__` (`glm.py:864-887`):
    /// `power = 0.0` (Normal), `link = LinkConfig::Auto`, `alpha = 1.0`,
    /// `max_iter = 100`, `tol = 1e-4`, `fit_intercept = true`,
    /// `solver = Solver::Lbfgs` (sklearn default). With the default
    /// `power = 0.0` and `Auto` link, the model is Normal/identity-link (OLS).
    #[must_use]
    pub fn new() -> Self {
        Self {
            power: 0.0,
            link: LinkConfig::Auto,
            alpha: F::one(),
            max_iter: 100,
            tol: F::from(1e-4).unwrap_or_else(F::epsilon),
            fit_intercept: true,
            solver: Solver::Lbfgs,
            warm_start: false,
            coef_init: None,
        }
    }

    /// Set the Tweedie power parameter.
    #[must_use]
    pub fn with_power(mut self, power: f64) -> Self {
        self.power = power;
        self
    }

    /// Set the link-function configuration (`Auto`/`Log`/`Identity`).
    ///
    /// Mirrors sklearn's `link={'auto','identity','log'}` (`glm.py:861`).
    #[must_use]
    pub fn with_link(mut self, link: LinkConfig) -> Self {
        self.link = link;
        self
    }

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

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

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

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

    /// Set the optimization [`Solver`], mirroring sklearn's `solver` parameter
    /// (`glm.py:140-145`, default `"lbfgs"`).
    ///
    /// ferrolearn fits via IRLS regardless of the value (R-DEV-7); both values
    /// produce the same observable `coef_` / `intercept_` (sklearn API parity,
    /// R-DEV-2).
    #[must_use]
    pub fn with_solver(mut self, solver: Solver) -> Self {
        self.solver = solver;
        self
    }

    /// Enable or disable warm-starting, mirroring sklearn's `warm_start`
    /// parameter (default `false`, `glm.py:146, :874`). See
    /// [`GLMRegressor::with_warm_start`] for the R-DEV-2 / R-DEV-7 rationale; the
    /// warm-start point is supplied explicitly via
    /// [`TweedieRegressor::with_coef_init`].
    #[must_use]
    pub fn with_warm_start(mut self, warm_start: bool) -> Self {
        self.warm_start = warm_start;
        self
    }

    /// Set the explicit warm-start initial point `(feature_coefficients,
    /// intercept)` — the ferrolearn analog of sklearn reusing `self.coef_` /
    /// `self.intercept_` (R-DEV-7, `glm.py:244-250`). Only consulted when
    /// `warm_start` is `true`; `coef`'s length must equal the number of features.
    #[must_use]
    pub fn with_coef_init(mut self, coef: Array1<F>, intercept: F) -> Self {
        self.coef_init = Some((coef, intercept));
        self
    }
}

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

impl<F: Float + Send + Sync + ScalarOperand + FromPrimitive + 'static> TweedieRegressor<F> {
    /// Fit the Tweedie GLM via IRLS with per-sample weights `sample_weight`.
    ///
    /// Mirrors sklearn's `TweedieRegressor.fit(X, y, sample_weight)`
    /// (`glm.py:170`, `:229-242`). The link is resolved from `link`/`power`
    /// as in [`Fit::fit`] (`glm.py:889-903`); see
    /// [`GLMRegressor::fit_with_sample_weight`] for the weighting semantics.
    /// An all-ones weight vector reproduces [`Fit::fit`] exactly.
    ///
    /// # Errors
    ///
    /// See [`GLMRegressor::fit_with_sample_weight`].
    pub fn fit_with_sample_weight(
        &self,
        x: &Array2<F>,
        y: &Array1<F>,
        sample_weight: &Array1<F>,
    ) -> Result<FittedGLMRegressor<F>, FerroError> {
        let link = self.link.resolve(self.power);
        fit_glm_irls(
            x,
            y,
            sample_weight,
            &GLMFamily::Tweedie(self.power),
            link,
            self.alpha,
            self.max_iter,
            self.tol,
            self.fit_intercept,
            self.warm_start,
            self.coef_init.as_ref(),
        )
    }
}

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

/// Cholesky solve for `A x = b`.
fn cholesky_solve<F: Float>(a: &Array2<F>, b: &Array1<F>) -> Result<Array1<F>, FerroError> {
    let n = a.nrows();
    let mut l = Array2::<F>::zeros((n, n));

    for i in 0..n {
        for j in 0..=i {
            let mut s = a[[i, j]];
            for k in 0..j {
                s = s - l[[i, k]] * l[[j, k]];
            }
            if i == j {
                if s <= F::zero() {
                    return Err(FerroError::NumericalInstability {
                        message: "Cholesky: matrix not positive definite".into(),
                    });
                }
                l[[i, j]] = s.sqrt();
            } else {
                l[[i, j]] = s / l[[j, j]];
            }
        }
    }

    let mut z = Array1::<F>::zeros(n);
    for i in 0..n {
        let mut s = b[i];
        for k in 0..i {
            s = s - l[[i, k]] * z[k];
        }
        z[i] = s / l[[i, i]];
    }

    let mut x_sol = Array1::<F>::zeros(n);
    for i in (0..n).rev() {
        let mut s = z[i];
        for k in (i + 1)..n {
            s = s - l[[k, i]] * x_sol[k];
        }
        x_sol[i] = s / l[[i, i]];
    }

    Ok(x_sol)
}

/// Gaussian elimination with partial pivoting.
fn gaussian_solve<F: Float>(
    n: usize,
    a: &Array2<F>,
    b: &Array1<F>,
) -> Result<Array1<F>, FerroError> {
    let mut aug = Array2::<F>::zeros((n, n + 1));
    for i in 0..n {
        for j in 0..n {
            aug[[i, j]] = a[[i, j]];
        }
        aug[[i, n]] = b[i];
    }

    for col in 0..n {
        let mut max_val = aug[[col, col]].abs();
        let mut max_row = col;
        for row in (col + 1)..n {
            let v = aug[[row, col]].abs();
            if v > max_val {
                max_val = v;
                max_row = row;
            }
        }

        if max_val < F::from(1e-12).unwrap_or_else(F::epsilon) {
            return Err(FerroError::NumericalInstability {
                message: "singular matrix in Gaussian elimination".into(),
            });
        }

        if max_row != col {
            for j in 0..=n {
                let tmp = aug[[col, j]];
                aug[[col, j]] = aug[[max_row, j]];
                aug[[max_row, j]] = tmp;
            }
        }

        let pivot = aug[[col, col]];
        for row in (col + 1)..n {
            let factor = aug[[row, col]] / pivot;
            for j in col..=n {
                let above = aug[[col, j]];
                aug[[row, j]] = aug[[row, j]] - factor * above;
            }
        }
    }

    let mut x_sol = Array1::<F>::zeros(n);
    for i in (0..n).rev() {
        let mut s = aug[[i, n]];
        for j in (i + 1)..n {
            s = s - aug[[i, j]] * x_sol[j];
        }
        if aug[[i, i]].abs() < F::from(1e-12).unwrap_or_else(F::epsilon) {
            return Err(FerroError::NumericalInstability {
                message: "near-zero pivot in back substitution".into(),
            });
        }
        x_sol[i] = s / aug[[i, i]];
    }

    Ok(x_sol)
}

/// Solve the weighted ridge system `(X^T W X + P) w = X^T W z`, where the
/// penalty matrix `P` adds the L2 regularization to the diagonal of the
/// feature columns only.
///
/// # Penalty scaling and the intercept (sklearn parity, `glm.py:229-258`)
///
/// scikit-learn minimizes the per-sample-MEAN half-deviance plus an L2 prior on
/// the feature coefficients (NOT the intercept):
///
/// ```text
/// J(w) = 1/(2*S) * sum_i s_i * deviance_i + 1/2 * alpha * ||coef||^2,
/// ```
///
/// with `S = sum_i s_i` (= `n_samples` for unweighted fits). Its stationarity
/// condition is `(1/S) * grad[sum 1/2 dev] + alpha * w_features = 0`.
///
/// The IRLS normal equations `X^T W X w = X^T W z` are the linearization of the
/// SUMMED half-deviance `sum_i 1/2 dev_i` (no `1/S` factor): `X^T W X` is the
/// summed-scale Hessian. To make those summed equations correspond to sklearn's
/// mean-scale objective we multiply the penalty by `S` (the sum of weights, =
/// `n_samples` unweighted) before adding it to the diagonal: the added penalty
/// is `weight_sum * alpha`, applied to feature columns only, leaving the
/// intercept (column 0 of the augmented design, when present) unpenalized.
///
/// `weight_sum` is `S = sum_i s_i`, the sum of the GLM `sample_weight`
/// (= `n_samples` for an all-ones / unweighted fit); `intercept_col` is
/// `Some(0)` when an intercept column was prepended to `x`, `None` otherwise.
fn weighted_ridge_solve<F: Float + FromPrimitive>(
    x: &Array2<F>,
    z: &Array1<F>,
    weights: &Array1<F>,
    alpha: F,
    weight_sum: F,
    intercept_col: Option<usize>,
) -> Result<Array1<F>, FerroError> {
    let (n_samples, n_features) = x.dim();

    let mut xtwx = Array2::<F>::zeros((n_features, n_features));
    let mut xtwz = Array1::<F>::zeros(n_features);

    for i in 0..n_samples {
        let wi = weights[i];
        let xi = x.row(i);
        for r in 0..n_features {
            xtwz[r] = xtwz[r] + wi * xi[r] * z[i];
            for c in 0..n_features {
                xtwx[[r, c]] = xtwx[[r, c]] + wi * xi[r] * xi[c];
            }
        }
    }

    // Add L2 regularization. The IRLS normal equations are at the SUMMED-deviance
    // scale, so to match sklearn's MEAN-deviance objective the diagonal penalty is
    // `weight_sum * alpha` (glm.py:229-242). The intercept column is excluded:
    // sklearn's `l2_reg_strength = self.alpha` weights only `||coef||^2`
    // (glm.py:258), never the intercept.
    let penalty = weight_sum * alpha;
    for i in 0..n_features {
        if Some(i) == intercept_col {
            continue;
        }
        xtwx[[i, i]] = xtwx[[i, i]] + penalty;
    }

    cholesky_solve(&xtwx, &xtwz).or_else(|_| gaussian_solve(n_features, &xtwx, &xtwz))
}

/// Core IRLS fitting logic shared by all GLM variants.
///
/// The IRLS update is parameterized by the [`Link`]: with linear predictor
/// `eta = X @ coef`, mean `mu = link.inverse(eta)` and link derivative
/// `dmu/deta`, the standard Fisher-scoring working weight and response are
/// `w = (dmu/deta)^2 / V(mu)` and `z = eta + (y - mu) / (dmu/deta)`
/// (`glm.py:362` for the inverse-link mapping). For [`Link::Log`]
/// (`dmu/deta = mu`) this is `w = mu^2 / V(mu)`, `z = eta + (y - mu)/mu`,
/// byte-identical to the previous log-only code. For [`Link::Identity`] with
/// `V(mu) = mu^0 = 1` (Normal/`power = 0`), `w = 1`, `z = y`, so IRLS reduces to
/// ordinary least squares.
///
/// # Warm start (R-DEV-2 + R-DEV-7)
///
/// sklearn's `warm_start=True` reuses the previous fit's `self.coef_` /
/// `self.intercept_` as the optimizer's starting point (`glm.py:243-254`).
/// ferrolearn's estimators are immutable (`fit(&self, ...)` returns a fresh
/// fitted object and never mutates `self`), so there is no `self.coef_` to reuse
/// across calls; the Rust-idiomatic analog (R-DEV-7) is an EXPLICIT initial point
/// supplied via `coef_init = Some((feature_coef, intercept))`. When
/// `warm_start == true` and `coef_init` is provided, the IRLS coefficient vector
/// (and the derived `eta` / `mu`) are seeded from `coef_init` instead of the
/// cold start (`coef = 0`, `eta = link(y)`). The feature-coefficient length must
/// equal `n_features` (else [`FerroError::ShapeMismatch`]). Otherwise the cold
/// start is kept byte-for-byte.
///
/// Because the penalized GLM objective is convex, the converged `coef_` /
/// `intercept_` are warm-start-INVARIANT: the init only changes the starting
/// point (and so the iteration count), never the optimum — so the warm-started
/// fit matches both the cold-start fit and the sklearn oracle (`glm.py:244-256`
/// reaches the same minimizer regardless of the seed).
#[allow(
    clippy::too_many_arguments,
    reason = "shared IRLS core threads the link and the warm-start init \
    alongside the family/penalty/convergence parameters; splitting into a config \
    struct would obscure the 1:1 mapping to sklearn's fit signature"
)]
fn fit_glm_irls<F: Float + Send + Sync + ScalarOperand + FromPrimitive + 'static>(
    x: &Array2<F>,
    y: &Array1<F>,
    sample_weight: &Array1<F>,
    family: &GLMFamily,
    link: Link,
    alpha: F,
    max_iter: usize,
    tol: F,
    fit_intercept: bool,
    warm_start: bool,
    coef_init: Option<&(Array1<F>, F)>,
) -> Result<FittedGLMRegressor<F>, FerroError> {
    let (n_samples, n_features_orig) = 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 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, mirroring sklearn's
    // `self._validate_data(X, y, ..., y_numeric=True)` (`glm.py:189-196`) which
    // keeps the default `force_all_finite=True`, so `check_array` rejects any
    // NaN or +/-inf in X OR y with a `ValueError` BEFORE the IRLS loop. sklearn
    // also validates `sample_weight` via `_check_sample_weight` (default
    // `force_all_finite=True`, `glm.py:211`). This is the SHARED IRLS entry that
    // every estimator (PoissonRegressor/GammaRegressor/TweedieRegressor/
    // GLMRegressor) routes through, so the check lands ONCE here.
    // `.iter().any(|v| !v.is_finite())` rejects both NaN and Inf (bounds-safe,
    // no panic, R-CODE-2), matching the crate idiom (`ridge.rs`). The finite
    // path is byte-identical (the guard never fires on finite input). The
    // sample_weight finiteness check MUST precede the non-negative check below:
    // sklearn's `_check_sample_weight` runs `check_array(force_all_finite=True)`
    // first, so a `-inf` weight raises the infinity `ValueError` (verified live),
    // NOT a non-negative error. This guard also precedes the per-family y-domain
    // validation so a non-finite y is reported as a finiteness failure
    // (`is_finite()` is `false` for NaN/Inf), matching sklearn's `check_array`
    // running before `in_y_true_range`.
    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(),
        });
    }

    // sklearn's GLM family does NOT reject finite negative `sample_weight`: it
    // calls `_check_sample_weight(sample_weight, X, dtype=loss_dtype)` WITHOUT
    // `only_non_negative=True` (`glm.py:211`), so a finite negative weight passes
    // validation and reaches the optimizer (verified live: `PoissonRegressor()
    // .fit(X, y, sample_weight=[..,-1.0])` fits, coef `[-0.04268902, 0.24080026]`).
    // Negative weights are mathematically valid for the GLM objective — the
    // deviance is averaged `sum_i s_i * deviance_i` (`glm.py:229-242`), so a
    // negative `s_i` simply NEGATES that sample's contribution. ferrolearn's IRLS
    // uses `sample_weight` LINEARLY (`weights[i] *= sample_weight[i]` feeding the
    // `X^T W X` / `X^T W z` accumulation in `weighted_ridge_solve`, NO `sqrt`), so
    // a negative weight flows through unchanged and reaches the same optimum as
    // sklearn's lbfgs. The non-finite guard above (#2261) still rejects NaN/Inf
    // weights (sklearn's `_check_sample_weight` runs `check_array(force_all_finite
    // =True)` first), matching exception parity for non-finite input.

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

    if alpha < F::zero() {
        return Err(FerroError::InvalidParameter {
            name: "alpha".into(),
            reason: "must be non-negative".into(),
        });
    }

    // Per-family target-domain validation, mirroring sklearn's
    //   `if not linear_loss.base_loss.in_y_true_range(y): raise ValueError(...)`
    // (`glm.py:221-225`). The valid `y` range is determined by the family's
    // Tweedie `power` (NOT the link — verified against the live oracle):
    //   * `power <= 0`  (Normal/identity)        -> y unconstrained
    //   * `0 < power < 2` (Poisson `power = 1`)  -> y >= 0   (closed at 0)
    //   * `power >= 2`    (Gamma `power = 2`)     -> y > 0    (open at 0)
    // Confirmed live-oracle behaviors: `GammaRegressor().fit(X, [0,1,2])` and
    // `TweedieRegressor(power=2.0).fit(X, [0,1,2])` raise ValueError;
    // `PoissonRegressor().fit(X, [-1,1,2])` raises; `TweedieRegressor(power=1.5)`
    // accepts y == 0.
    let min_y = F::from(1e-10).unwrap_or_else(F::epsilon);
    let y_domain = YDomain::for_power(family.domain_power());
    if y.iter().any(|&yi| !y_domain.contains(yi)) {
        return Err(FerroError::InvalidParameter {
            name: "y".into(),
            reason: format!(
                "Some value(s) of y are out of the valid range of the loss '{}'.",
                y_domain.loss_name()
            ),
        });
    }

    // Build design matrix (optionally prepend intercept column).
    let n_cols = if fit_intercept {
        n_features_orig + 1
    } else {
        n_features_orig
    };

    let mut x_design = Array2::<F>::zeros((n_samples, n_cols));
    if fit_intercept {
        for i in 0..n_samples {
            x_design[[i, 0]] = F::one();
            for j in 0..n_features_orig {
                x_design[[i, j + 1]] = x[[i, j]];
            }
        }
    } else {
        x_design.assign(x);
    }

    // The intercept (column 0 of the augmented design when fitting one) is
    // excluded from the L2 penalty, matching sklearn (glm.py:258).
    let intercept_col = if fit_intercept { Some(0) } else { None };

    // Sum of sample weights `S = sum_i s_i`. For an unweighted GLM every weight
    // is 1, so this is `n_samples` (byte-identical to the previous unweighted
    // path). It scales the L2 penalty so the summed-deviance IRLS normal
    // equations minimize sklearn's mean-deviance objective normalized by
    // `sum(sample_weight)` (glm.py:229-242).
    let weight_sum = sample_weight.iter().fold(F::zero(), |acc, &si| acc + si);

    // For the log link, clamp y away from 0 so `ln(y)` and `mu` stay finite.
    // For the identity link y is used as-is (mu = eta = y, no positivity
    // constraint).
    let y_safe: Array1<F> = match link {
        Link::Log => y.mapv(|v| if v < min_y { min_y } else { v }),
        Link::Identity => y.clone(),
    };

    // Initialise mu = y_safe, eta = g(mu) = link(mu): log link → ln(y),
    // identity link → y.
    let mut mu: Array1<F> = y_safe.clone();
    let mut eta: Array1<F> = match link {
        Link::Log => y_safe.mapv(|v| v.ln()),
        Link::Identity => y_safe.clone(),
    };

    // Coefficient initialization.
    //
    // COLD START (default, `warm_start == false` or no `coef_init`): feature
    // coefficients are 0; when `fit_intercept` the intercept entry is seeded at
    // `link.link(weighted_mean(y))` (REQ-5, `glm.py:251-256`), and `eta`/`mu`
    // are recomputed from that seed (so `eta == link(mean_w(y))` constant). When
    // `fit_intercept` is false the intercept seed is skipped and `eta`/`mu` keep
    // the per-sample `link(y)` seed above. The objective is convex, so this
    // changes only the starting point (and the iteration count), never the
    // converged optimum — the cold fit still matches the sklearn oracle.
    //
    // WARM START (R-DEV-7 analog of sklearn's `warm_start=True`,
    // `glm.py:243-254`): seed the IRLS coefficient vector from the explicit
    // `coef_init = (feature_coef, intercept)` instead, then recompute
    // `eta = X_design @ coef` and `mu = link.inverse(eta)` so the first IRLS
    // weights/working-response are formed at the supplied point (sklearn likewise
    // hands `coef` straight to the optimizer). When `warm_start && coef_init`,
    // the explicit init takes precedence over the REQ-5 intercept seed.
    let mut coef = Array1::<F>::zeros(n_cols);
    if warm_start && let Some((feature_coef, intercept_init)) = coef_init {
        if feature_coef.len() != n_features_orig {
            return Err(FerroError::ShapeMismatch {
                expected: vec![n_features_orig],
                actual: vec![feature_coef.len()],
                context: "coef_init feature-coefficient length must match \
                          number of features in X"
                    .into(),
            });
        }
        // Place the supplied coefficients into the (optionally
        // intercept-augmented) design layout: column 0 is the intercept when
        // `fit_intercept`, features follow. When `fit_intercept == false` the
        // supplied intercept is ignored (the model has no intercept term),
        // matching sklearn's `coef = self.coef_` (no intercept) branch
        // (`glm.py:248-249`).
        if fit_intercept {
            coef[0] = *intercept_init;
            for (j, &cj) in feature_coef.iter().enumerate() {
                coef[j + 1] = cj;
            }
        } else {
            for (j, &cj) in feature_coef.iter().enumerate() {
                coef[j] = cj;
            }
        }
        // Recompute eta/mu at the warm-start point so the first IRLS step is
        // formed there (link.inverse with the same clamps as the loop body).
        eta = x_design.dot(&coef);
        match link {
            Link::Log => {
                let hi = F::from(20.0).unwrap_or_else(F::max_value);
                let lo = F::zero() - hi;
                let mmu = F::from(1e-10).unwrap_or_else(F::epsilon);
                let xmu = F::from(1e10).unwrap_or_else(F::max_value);
                for i in 0..n_samples {
                    let eta_i = eta[i].max(lo).min(hi);
                    eta[i] = eta_i;
                    mu[i] = link.inverse(eta_i).max(mmu).min(xmu);
                }
            }
            Link::Identity => {
                for i in 0..n_samples {
                    mu[i] = link.inverse(eta[i]);
                }
            }
        }
    } else if fit_intercept {
        // COLD-START INTERCEPT INITIALIZATION (REQ-5, `glm.py:251-256`).
        //
        // sklearn cold-starts with `coef = init_zero_coef(X)` and, when
        // `fit_intercept`, seeds the intercept entry at
        //   `coef[-1] = link.link(np.average(y, weights=sample_weight))`
        // (the feature coefficients stay 0). For the log link this is
        // `intercept_init = log(weighted_mean(y))`; for the identity link it is
        // `weighted_mean(y)`. We mirror this exactly: the intercept is design
        // column 0, so `coef[0] = link.link(weighted_mean(y))`.
        //
        // `weighted_mean(y) = Σ(s_i·y_i) / Σ(s_i)` (= plain mean for an all-ones
        // sample_weight). `weight_sum = Σ s_i` was computed above.
        //
        // Edge case (R-CODE-2, R-DEV-1): the seed must lie in the link's domain.
        // For the log link `weighted_mean(y)` can be 0 (e.g. all-zero Poisson
        // `y`, which is IN the Poisson domain `y >= 0`), where `log(0) = -inf` —
        // sklearn likewise computes `link.link(0) = -inf` there (`glm.py:254`)
        // and its lbfgs optimum STAYS at `-inf`: the penalized data deviance is
        // minimized as `mu -> 0` (`eta -> -inf`), and the L2 penalty is on the
        // feature coefficients only (`l2_reg_strength = self.alpha`,
        // `glm.py:258`), which are 0. So sklearn returns `intercept_ = -inf`,
        // `coef_ = 0`, and `predict = inverse_link(X @ 0 - inf) = exp(-inf) =
        // 0.0` EXACTLY for every sample. We short-circuit to that same degenerate
        // optimum: IRLS cannot iterate from `mu = 0` (the Fisher weights blow up),
        // so we skip the loop entirely and return `coef_ = 0`, `intercept_ =
        // link(weighted_mean(y))` (the non-finite seed). This is alpha-INVARIANT
        // (the penalty on coef is 0 since coef = 0; the deviance drives
        // `eta -> -inf` regardless of alpha) and fires ONLY for a non-finite seed
        // (log link + `mean(y) == 0`) — the exact all-zero case. Gamma requires
        // `y > 0` so `mean(y) > 0` always (domain guard rejects `y == 0`); the
        // identity link gives `link(0) = 0` (finite), so neither hits this path.
        let weighted_y_sum = y
            .iter()
            .zip(sample_weight.iter())
            .fold(F::zero(), |acc, (&yi, &si)| acc + si * yi);
        let denom = if weight_sum > F::zero() {
            weight_sum
        } else {
            F::from(n_samples).unwrap_or_else(F::one)
        };
        let weighted_mean_y = weighted_y_sum / denom;
        let intercept_init = link.link(weighted_mean_y);

        if !intercept_init.is_finite() {
            // Degenerate optimum (all-zero-y log-link case): short-circuit to
            // sklearn's lbfgs landing point `coef_ = 0`, `intercept_ = -inf`
            // (`glm.py:254`), skipping IRLS (which cannot iterate from `mu = 0`).
            // `predict = exp(-inf) = 0.0` then matches sklearn exactly. The
            // feature coefficients are already 0 (cold start); set the intercept
            // column to the non-finite seed and return immediately. `n_iter_ = 0`
            // records that the optimum was reached without any IRLS iteration.
            coef[0] = intercept_init;
            let intercept = coef[0];
            let coefficients = Array1::from_iter(coef.iter().skip(1).copied());
            return Ok(FittedGLMRegressor {
                coefficients,
                intercept,
                link,
                family: *family,
                n_iter: 0,
            });
        }

        {
            coef[0] = intercept_init;
            // Recompute eta/mu from the seeded coef (feature coefs are 0, so
            // `eta = intercept_init` for every sample, i.e. `eta = link(mean y)`).
            eta = x_design.dot(&coef);
            match link {
                Link::Log => {
                    let hi = F::from(20.0).unwrap_or_else(F::max_value);
                    let lo = F::zero() - hi;
                    let mmu = F::from(1e-10).unwrap_or_else(F::epsilon);
                    let xmu = F::from(1e10).unwrap_or_else(F::max_value);
                    for i in 0..n_samples {
                        let eta_i = eta[i].max(lo).min(hi);
                        eta[i] = eta_i;
                        mu[i] = link.inverse(eta_i).max(mmu).min(xmu);
                    }
                }
                Link::Identity => {
                    for i in 0..n_samples {
                        mu[i] = link.inverse(eta[i]);
                    }
                }
            }
        }
        // (The non-finite-seed branch above returns early, so there is no
        // finite-init fall-through here.)
    }

    let min_mu = F::from(1e-10).unwrap_or_else(F::epsilon);
    let max_mu = F::from(1e10).unwrap_or_else(F::max_value);

    // Count the IRLS iterations actually run (sklearn's `n_iter_`, `glm.py:283`).
    // At least one iteration always runs (`max_iter >= 1`); on convergence we
    // break after the iteration that satisfied the tolerance.
    let mut n_iter = 0usize;
    for _iter in 0..max_iter {
        n_iter += 1;
        let coef_old = coef.clone();

        // Compute IRLS weights and working response.
        let mut weights = Array1::<F>::zeros(n_samples);
        let mut z = Array1::<F>::zeros(n_samples);

        for i in 0..n_samples {
            // IRLS (Fisher scoring) with the configured link:
            //   dmu/deta  : Log => mu, Identity => 1
            //   weight w  = (dmu/deta)^2 / V(mu)
            //   response z = eta + (y - mu) / (dmu/deta)
            // For Log this is `w = mu^2/V(mu)`, `z = eta + (y - mu)/mu`,
            // byte-identical to the previous log-only code (clamped `mu_i`
            // throughout). For Identity + power=0 (V=1): w=1, z=y => OLS.
            match link {
                Link::Log => {
                    let mu_i = mu[i].max(min_mu).min(max_mu);
                    let var_i = family.variance(mu_i).max(min_mu);
                    let g_prime = F::one() / mu_i; // derivative of log link
                    z[i] = eta[i] + (y_safe[i] - mu_i) * g_prime;
                    weights[i] = F::one() / (g_prime * g_prime * var_i);
                }
                Link::Identity => {
                    let mu_i = mu[i];
                    // V(mu) for the identity link can see mu <= 0 (eta is
                    // unbounded); for power=0, V(mu)=mu^0=1 always. Clamp the
                    // magnitude for non-zero powers so V stays finite/positive.
                    let var_i = family.variance(mu_i.abs().max(min_mu)).max(min_mu);
                    let dmu_deta = link.dmu_deta(mu_i); // = 1
                    z[i] = eta[i] + (y_safe[i] - mu_i) / dmu_deta;
                    weights[i] = dmu_deta * dmu_deta / var_i;
                }
            }
            // Clamp the Fisher (GLM) working weight.
            if weights[i] < min_mu {
                weights[i] = min_mu;
            }
            // Apply the per-sample weight `s_i`: the `W` diagonal entry for
            // sample i becomes `s_i * w_irls,i` (standard weighted IRLS).
            // sklearn weights the deviance average by `sample_weight`
            // (glm.py:229-242); the working response `z_i` is unchanged. For
            // all-ones weights this is a no-op (byte-identical unweighted path).
            weights[i] = weights[i] * sample_weight[i];
        }

        // Solve weighted ridge. `weights` now carry `s_i * w_irls,i`;
        // `weight_sum = S = sum_i s_i` (= n_samples for an all-ones fit). The
        // penalty is scaled by `S` so the summed-deviance normal equations
        // minimize sklearn's sample-weight-averaged deviance objective
        // (glm.py:229-242). The intercept column (column 0 of the
        // augmented design, present iff `fit_intercept`) is left unpenalized
        // (glm.py:258, `l2_reg_strength = self.alpha` weighs only `||coef||^2`).
        coef = weighted_ridge_solve(&x_design, &z, &weights, alpha, weight_sum, intercept_col)?;

        // Update eta = X @ coef and mu = link.inverse(eta).
        eta = x_design.dot(&coef);
        match link {
            Link::Log => {
                let hi = F::from(20.0).unwrap_or_else(F::max_value);
                let lo = F::zero() - hi;
                for i in 0..n_samples {
                    // Clamp eta to prevent overflow in exp.
                    let eta_i = eta[i].max(lo).min(hi);
                    eta[i] = eta_i;
                    mu[i] = link.inverse(eta_i).max(min_mu).min(max_mu);
                }
            }
            Link::Identity => {
                // Identity link: eta is unbounded; mu = eta (no exp clamp).
                for i in 0..n_samples {
                    mu[i] = link.inverse(eta[i]);
                }
            }
        }

        // Check convergence.
        let max_change = coef
            .iter()
            .zip(coef_old.iter())
            .map(|(&c, &co)| (c - co).abs())
            .fold(F::zero(), |a, b| if b > a { b } else { a });

        if max_change < tol {
            break;
        }
    }

    // Extract intercept and feature coefficients.
    let (intercept, coefficients) = if fit_intercept {
        let intercept = coef[0];
        let coefficients = Array1::from_iter(coef.iter().skip(1).copied());
        (intercept, coefficients)
    } else {
        (F::zero(), coef)
    };

    Ok(FittedGLMRegressor {
        coefficients,
        intercept,
        link,
        family: *family,
        n_iter,
    })
}

// ---------------------------------------------------------------------------
// Fit — GLMRegressor
// ---------------------------------------------------------------------------

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

    /// Fit the GLM via IRLS.
    ///
    /// # Errors
    ///
    /// - [`FerroError::ShapeMismatch`] — sample count mismatch.
    /// - [`FerroError::InsufficientSamples`] — zero samples.
    /// - [`FerroError::InvalidParameter`] — negative alpha or negative y.
    fn fit(&self, x: &Array2<F>, y: &Array1<F>) -> Result<FittedGLMRegressor<F>, FerroError> {
        // GLMRegressor's families are all log-link (Poisson/Gamma/Tweedie>0).
        // The Tweedie identity link is exposed only through TweedieRegressor's
        // `link` configuration (sklearn similarly only exposes `link` on
        // `TweedieRegressor`, `glm.py:861`). The default (no sample_weight) path
        // delegates with an all-ones weight vector, matching sklearn's
        // `_check_sample_weight` default (`glm.py:208-211`); `weight_sum` then
        // equals `n_samples`, so this is byte-identical to the unweighted fit.
        self.fit_with_sample_weight(x, y, &Array1::ones(x.nrows()))
    }
}

// ---------------------------------------------------------------------------
// Fit — PoissonRegressor
// ---------------------------------------------------------------------------

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

    /// Fit the Poisson GLM via IRLS.
    ///
    /// # Errors
    ///
    /// See [`GLMRegressor::fit`].
    fn fit(&self, x: &Array2<F>, y: &Array1<F>) -> Result<FittedGLMRegressor<F>, FerroError> {
        // Poisson uses the log link only (`HalfPoissonLoss`, `glm.py:589-590`).
        // Default (no sample_weight) path: all-ones weights => byte-identical to
        // the unweighted fit (`weight_sum = n_samples`).
        self.fit_with_sample_weight(x, y, &Array1::ones(x.nrows()))
    }
}

// ---------------------------------------------------------------------------
// Fit — GammaRegressor
// ---------------------------------------------------------------------------

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

    /// Fit the Gamma GLM via IRLS.
    ///
    /// # Errors
    ///
    /// See [`GLMRegressor::fit`].
    fn fit(&self, x: &Array2<F>, y: &Array1<F>) -> Result<FittedGLMRegressor<F>, FerroError> {
        // Gamma uses the log link only (`HalfGammaLoss`, `glm.py:721-722`).
        // Default (no sample_weight) path: all-ones weights => byte-identical to
        // the unweighted fit (`weight_sum = n_samples`).
        self.fit_with_sample_weight(x, y, &Array1::ones(x.nrows()))
    }
}

// ---------------------------------------------------------------------------
// Fit — TweedieRegressor
// ---------------------------------------------------------------------------

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

    /// Fit the Tweedie GLM via IRLS.
    ///
    /// # Errors
    ///
    /// See [`GLMRegressor::fit`].
    fn fit(&self, x: &Array2<F>, y: &Array1<F>) -> Result<FittedGLMRegressor<F>, FerroError> {
        // Resolve the link from the configuration and Tweedie power, mirroring
        // `TweedieRegressor._get_loss` (`glm.py:889-903`): `auto` selects
        // identity for `power <= 0` (Normal/OLS) and log for `power > 0`.
        // Default (no sample_weight) path: all-ones weights => byte-identical to
        // the unweighted fit (`weight_sum = n_samples`).
        self.fit_with_sample_weight(x, y, &Array1::ones(x.nrows()))
    }
}

// ---------------------------------------------------------------------------
// Predict / HasCoefficients / Pipeline — FittedGLMRegressor
// ---------------------------------------------------------------------------

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

    /// Predict using the fitted GLM.
    ///
    /// Computes `link.inverse(X @ coefficients + intercept)` (`glm.py:362`):
    /// `exp(...)` for a [`Link::Log`] model (Poisson/Gamma/Tweedie with
    /// `power > 0`), and the raw linear predictor `X @ coef + intercept` for a
    /// [`Link::Identity`] model (Tweedie with `power <= 0`, Normal).
    ///
    /// # 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> {
        if x.ncols() != self.coefficients.len() {
            return Err(FerroError::ShapeMismatch {
                expected: vec![self.coefficients.len()],
                actual: vec![x.ncols()],
                context: "number of features must match fitted model".into(),
            });
        }
        let eta = x.dot(&self.coefficients) + self.intercept;
        let link = self.link;
        Ok(eta.mapv(|v| link.inverse(v)))
    }
}

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

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

// Pipeline integration for GLMRegressor.
impl<F> PipelineEstimator<F> for GLMRegressor<F>
where
    F: Float + FromPrimitive + 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 FittedGLMRegressor<F>
where
    F: Float + ScalarOperand + Send + Sync + 'static,
{
    fn predict_pipeline(&self, x: &Array2<F>) -> Result<Array1<F>, FerroError> {
        self.predict(x)
    }
}

// Pipeline integration for PoissonRegressor.
impl<F> PipelineEstimator<F> for PoissonRegressor<F>
where
    F: Float + FromPrimitive + 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))
    }
}

// Pipeline integration for GammaRegressor.
impl<F> PipelineEstimator<F> for GammaRegressor<F>
where
    F: Float + FromPrimitive + 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))
    }
}

// Pipeline integration for TweedieRegressor.
impl<F> PipelineEstimator<F> for TweedieRegressor<F>
where
    F: Float + FromPrimitive + 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))
    }
}

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

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

    // ---- GLMRegressor ----

    #[test]
    fn test_glm_poisson_defaults() {
        let m = GLMRegressor::<f64>::new(GLMFamily::Poisson);
        assert_relative_eq!(m.alpha, 1.0);
        assert_eq!(m.max_iter, 100);
        assert!(m.fit_intercept);
    }

    #[test]
    fn test_glm_builder() {
        let m = GLMRegressor::<f64>::new(GLMFamily::Gamma)
            .with_alpha(0.5)
            .with_max_iter(200)
            .with_tol(1e-6)
            .with_fit_intercept(false);
        assert_relative_eq!(m.alpha, 0.5);
        assert_eq!(m.max_iter, 200);
        assert!(!m.fit_intercept);
    }

    #[test]
    fn test_glm_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];
        assert!(
            GLMRegressor::<f64>::new(GLMFamily::Poisson)
                .fit(&x, &y)
                .is_err()
        );
    }

    #[test]
    fn test_glm_negative_alpha() {
        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];
        assert!(
            GLMRegressor::<f64>::new(GLMFamily::Poisson)
                .with_alpha(-1.0)
                .fit(&x, &y)
                .is_err()
        );
    }

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

        let fitted = GLMRegressor::<f64>::new(GLMFamily::Poisson)
            .with_alpha(0.0)
            .with_max_iter(200)
            .fit(&x, &y)
            .unwrap();
        let preds = fitted.predict(&x).unwrap();
        assert_eq!(preds.len(), 4);
        // Predictions should be positive.
        for &p in preds.iter() {
            assert!(p > 0.0);
        }
    }

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

        let fitted = GLMRegressor::<f64>::new(GLMFamily::Gamma)
            .with_alpha(0.0)
            .with_max_iter(200)
            .fit(&x, &y)
            .unwrap();
        let preds = fitted.predict(&x).unwrap();
        assert_eq!(preds.len(), 4);
        for &p in preds.iter() {
            assert!(p > 0.0);
        }
    }

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

        let fitted = GLMRegressor::<f64>::new(GLMFamily::Tweedie(1.5))
            .with_alpha(0.0)
            .with_max_iter(200)
            .fit(&x, &y)
            .unwrap();
        let preds = fitted.predict(&x).unwrap();
        assert_eq!(preds.len(), 4);
        for &p in preds.iter() {
            assert!(p > 0.0);
        }
    }

    #[test]
    fn test_glm_predict_feature_mismatch() {
        let x = 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 = GLMRegressor::<f64>::new(GLMFamily::Poisson)
            .fit(&x, &y)
            .unwrap();
        let x_bad = Array2::from_shape_vec((3, 1), vec![1.0, 2.0, 3.0]).unwrap();
        assert!(fitted.predict(&x_bad).is_err());
    }

    #[test]
    fn test_glm_has_coefficients() {
        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 = GLMRegressor::<f64>::new(GLMFamily::Poisson)
            .fit(&x, &y)
            .unwrap();
        assert_eq!(fitted.coefficients().len(), 2);
    }

    #[test]
    fn test_glm_pipeline() {
        let x = Array2::from_shape_vec((4, 1), vec![1.0, 2.0, 3.0, 4.0]).unwrap();
        let y = array![2.0, 5.0, 10.0, 20.0];
        let model = GLMRegressor::<f64>::new(GLMFamily::Poisson).with_alpha(0.0);
        let fitted = model.fit_pipeline(&x, &y).unwrap();
        let preds = fitted.predict_pipeline(&x).unwrap();
        assert_eq!(preds.len(), 4);
    }

    // ---- PoissonRegressor ----

    #[test]
    fn test_poisson_defaults() {
        let m = PoissonRegressor::<f64>::new();
        assert_relative_eq!(m.alpha, 1.0);
        assert_eq!(m.max_iter, 100);
        assert!(m.fit_intercept);
    }

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

        let fitted = PoissonRegressor::<f64>::new()
            .with_alpha(0.0)
            .with_max_iter(200)
            .fit(&x, &y)
            .unwrap();
        let preds = fitted.predict(&x).unwrap();
        assert_eq!(preds.len(), 4);
        for &p in preds.iter() {
            assert!(p > 0.0);
        }
    }

    #[test]
    fn test_poisson_pipeline() {
        let x = Array2::from_shape_vec((4, 1), vec![1.0, 2.0, 3.0, 4.0]).unwrap();
        let y = array![2.0, 5.0, 10.0, 20.0];
        let fitted = PoissonRegressor::<f64>::new()
            .with_alpha(0.0)
            .fit_pipeline(&x, &y)
            .unwrap();
        let preds = fitted.predict_pipeline(&x).unwrap();
        assert_eq!(preds.len(), 4);
    }

    // ---- GammaRegressor ----

    #[test]
    fn test_gamma_defaults() {
        let m = GammaRegressor::<f64>::new();
        assert_relative_eq!(m.alpha, 1.0);
        assert_eq!(m.max_iter, 100);
    }

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

        let fitted = GammaRegressor::<f64>::new()
            .with_alpha(0.0)
            .with_max_iter(200)
            .fit(&x, &y)
            .unwrap();
        let preds = fitted.predict(&x).unwrap();
        assert_eq!(preds.len(), 4);
        for &p in preds.iter() {
            assert!(p > 0.0);
        }
    }

    #[test]
    fn test_gamma_pipeline() {
        let x = Array2::from_shape_vec((4, 1), vec![1.0, 2.0, 3.0, 4.0]).unwrap();
        let y = array![2.0, 5.0, 10.0, 20.0];
        let fitted = GammaRegressor::<f64>::new()
            .with_alpha(0.0)
            .fit_pipeline(&x, &y)
            .unwrap();
        let preds = fitted.predict_pipeline(&x).unwrap();
        assert_eq!(preds.len(), 4);
    }

    // ---- TweedieRegressor ----

    #[test]
    fn test_tweedie_defaults() {
        let m = TweedieRegressor::<f64>::new();
        // sklearn TweedieRegressor default power=0.0 (Normal), link='auto'
        // (glm.py:867, :870).
        assert_relative_eq!(m.power, 0.0);
        assert_eq!(m.link, LinkConfig::Auto);
        assert_relative_eq!(m.alpha, 1.0);
    }

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

        let fitted = TweedieRegressor::<f64>::new()
            .with_power(1.5)
            .with_alpha(0.0)
            .with_max_iter(200)
            .fit(&x, &y)
            .unwrap();
        let preds = fitted.predict(&x).unwrap();
        assert_eq!(preds.len(), 4);
        for &p in preds.iter() {
            assert!(p > 0.0);
        }
    }

    #[test]
    fn test_tweedie_pipeline() {
        let x = Array2::from_shape_vec((4, 1), vec![1.0, 2.0, 3.0, 4.0]).unwrap();
        let y = array![2.0, 5.0, 10.0, 20.0];
        let fitted = TweedieRegressor::<f64>::new()
            .with_alpha(0.0)
            .fit_pipeline(&x, &y)
            .unwrap();
        let preds = fitted.predict_pipeline(&x).unwrap();
        assert_eq!(preds.len(), 4);
    }

    // ---- Link ----

    #[test]
    fn test_link_inverse() {
        assert_relative_eq!(Link::Log.inverse(0.0_f64), 1.0);
        assert_relative_eq!(Link::Identity.inverse(3.5_f64), 3.5);
    }

    #[test]
    fn test_link_config_resolve_auto() {
        // glm.py:889-893: auto -> identity for power<=0, log for power>0.
        assert_eq!(LinkConfig::Auto.resolve(0.0), Link::Identity);
        assert_eq!(LinkConfig::Auto.resolve(-1.0), Link::Identity);
        assert_eq!(LinkConfig::Auto.resolve(1.5), Link::Log);
        assert_eq!(LinkConfig::Log.resolve(0.0), Link::Log);
        assert_eq!(LinkConfig::Identity.resolve(2.0), Link::Identity);
    }

    #[test]
    fn test_tweedie_with_link_builder() {
        let m = TweedieRegressor::<f64>::new().with_link(LinkConfig::Log);
        assert_eq!(m.link, LinkConfig::Log);
    }

    // ---- Solver (sklearn API parity, glm.py:140-145) ----

    #[test]
    fn test_solver_default_lbfgs() {
        // sklearn default solver='lbfgs' (glm.py:155).
        assert_eq!(
            GLMRegressor::<f64>::new(GLMFamily::Poisson).solver,
            Solver::Lbfgs
        );
        assert_eq!(PoissonRegressor::<f64>::new().solver, Solver::Lbfgs);
        assert_eq!(GammaRegressor::<f64>::new().solver, Solver::Lbfgs);
        assert_eq!(TweedieRegressor::<f64>::new().solver, Solver::Lbfgs);
    }

    #[test]
    fn test_with_solver_builder() {
        assert_eq!(
            PoissonRegressor::<f64>::new()
                .with_solver(Solver::NewtonCholesky)
                .solver,
            Solver::NewtonCholesky
        );
        assert_eq!(
            GLMRegressor::<f64>::new(GLMFamily::Gamma)
                .with_solver(Solver::NewtonCholesky)
                .solver,
            Solver::NewtonCholesky
        );
    }

    // ---- Variance function ----

    #[test]
    fn test_variance_poisson() {
        let v = GLMFamily::Poisson.variance(3.0_f64);
        assert_relative_eq!(v, 3.0);
    }

    #[test]
    fn test_variance_gamma() {
        let v = GLMFamily::Gamma.variance(3.0_f64);
        assert_relative_eq!(v, 9.0);
    }

    #[test]
    fn test_variance_tweedie() {
        let v = GLMFamily::Tweedie(1.5).variance(4.0_f64);
        assert_relative_eq!(v, 4.0_f64.powf(1.5), epsilon = 1e-10);
    }

    #[test]
    fn test_glm_negative_y() {
        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];
        assert!(
            GLMRegressor::<f64>::new(GLMFamily::Poisson)
                .fit(&x, &y)
                .is_err()
        );
    }
}