limma-rust 0.1.0

//! Empirical Bayes moderation. Port of limma's `eBayes`/`.ebayes`
//! (`ebayes.R`), `squeezeVar` (`squeezeVar.R`), `fitFDist` (`fitFDist.R`),
//! `fitFDistRobustly` (`fitFDistRobustly.R`) and the `fstat.only` path of
//! `classifyTestsF` (`decidetests.R`). `trend` and `robust` are independently
//! selectable, including their combination (trended-robust via `loessFit`).
//!
//! `fitFDist` is ported in full, including the trended (covariate) fit and its
//! `notallok` branch, which drops unusable genes from the spline fit and then
//! extrapolates the trend back to them (`predict.ns`). The exact-likelihood
//! `fitFDistUnequalDF1` is ported in full as well — the constant and trended
//! (covariate) non-robust fits plus the `robust` sub-path (FDR-based outlier
//! detection, the `Recall` refit weighted by FDR, and the per-gene
//! `df2.shrunk` machinery) — so robust empirical Bayes works even when the
//! residual d.f. differ across genes.

use anyhow::{bail, Result};
use ndarray::{Array1, Array2};

use crate::fit::MArrayLM;
use crate::linalg::{cov2cor, eigh, qr_econ, solve_upper};
use crate::lowess::{loess_fit, loess_fit_unweighted};
use crate::special::{
    betai, chi2_sf, chisq_pdf, chisq_qf, digamma, f_lsf, f_pdf, f_qf, f_sf, gauss_legendre_01,
    ln_gamma, logmdigamma, t_isf, t_sf, t_two_sided_pvalue, trigamma, trigamma_inverse,
};
use crate::splines::NsBasis;
use crate::toptable::p_adjust_bh;
use crate::voom::choose_lowess_span;

/// Median of a slice (R `median` semantics: mean of the two central order
/// statistics for even length). Ignores nothing — caller pre-filters.
fn median(xs: &[f64]) -> f64 {
    let mut v: Vec<f64> = xs.to_vec();
    v.sort_by(|a, b| a.partial_cmp(b).unwrap());
    let n = v.len();
    if n == 0 {
        return f64::NAN;
    }
    if n % 2 == 1 {
        v[n / 2]
    } else {
        0.5 * (v[n / 2 - 1] + v[n / 2])
    }
}

/// Moment estimation of `(scale, df2)` of a scaled F-distribution from a sample
/// `x` with first d.f. `df1`. Port of `fitFDist` with `covariate=None`.
pub(crate) fn fit_fdist(x: &Array1<f64>, df1: &Array1<f64>) -> (f64, f64) {
    let n = x.len();
    if n == 0 {
        return (f64::NAN, f64::NAN);
    }
    if n == 1 {
        return (x[0], 0.0);
    }

    // Keep finite, non-negative variances with usable degrees of freedom.
    let mut xs: Vec<f64> = Vec::new();
    let mut d1: Vec<f64> = Vec::new();
    for i in 0..n {
        let ok = df1[i].is_finite() && df1[i] > 1e-15 && x[i].is_finite() && x[i] > -1e-15;
        if ok {
            xs.push(x[i]);
            d1.push(df1[i]);
        }
    }
    let nok = xs.len();
    if nok == 0 {
        return (f64::NAN, f64::NAN);
    }
    if nok == 1 {
        return (xs[0], 0.0);
    }

    // Offset exact/near zeros away from zero.
    for v in xs.iter_mut() {
        if *v < 0.0 {
            *v = 0.0;
        }
    }
    let mut m = median(&xs);
    if m == 0.0 {
        m = 1.0;
    }
    let floor = 1e-5 * m;
    for v in xs.iter_mut() {
        if *v < floor {
            *v = floor;
        }
    }

    // Work on log(F).
    let e: Vec<f64> = xs
        .iter()
        .zip(d1.iter())
        .map(|(&xv, &dv)| xv.ln() - digamma(dv / 2.0) + (dv / 2.0).ln())
        .collect();
    let emean: f64 = e.iter().sum::<f64>() / nok as f64;
    let mut evar: f64 = e.iter().map(|&ev| (ev - emean).powi(2)).sum::<f64>() / (nok as f64 - 1.0);
    let tri_mean: f64 = d1.iter().map(|&dv| trigamma(dv / 2.0)).sum::<f64>() / nok as f64;
    evar -= tri_mean;

    if evar > 0.0 {
        let df2 = 2.0 * trigamma_inverse(evar);
        let s20 = (emean + digamma(df2 / 2.0) - (df2 / 2.0).ln()).exp();
        (s20, df2)
    } else {
        // Pooled variance is the MLE of the scale in this limiting case.
        let s20 = xs.iter().sum::<f64>() / nok as f64;
        (s20, f64::INFINITY)
    }
}

/// Trended moment estimation: `fitFDist(x, df1, covariate)`. The log-variance
/// trend is modelled by a natural cubic spline in the covariate (R uses
/// `splines::ns(covariate, df, intercept=TRUE)`), so the scale `s2.0`
/// (`fit$scale`) becomes per-gene; the second d.f. (`fit$df2`) stays scalar.
/// Returns `(scale[n], df2)`.
///
/// Genes with a non-usable variance or `df1` (non-finite, non-positive, or zero
/// d.f.) are dropped from the spline fit (`ok` mask), exactly as in `fitFDist`;
/// the trend is then *extrapolated* back to those genes (`predict.ns`), so every
/// gene still receives a scale. The covariate may not contain `NA`, but `±Inf`
/// covariate values are clamped to the finite range `±1` (R's behaviour).
fn fit_fdist_trend(
    x: &Array1<f64>,
    df1: &Array1<f64>,
    covariate: &Array1<f64>,
) -> Result<(Array1<f64>, f64)> {
    let n = x.len();
    if n == 0 {
        return Ok((Array1::from(Vec::<f64>::new()), f64::NAN));
    }
    if n == 1 {
        return Ok((Array1::from_elem(1, x[0]), 0.0));
    }

    let df1_scalar = df1.len() == 1;
    if !df1_scalar && df1.len() != n {
        bail!("x and df1 have different lengths");
    }
    let df1_at = |i: usize| if df1_scalar { df1[0] } else { df1[i] };

    // `ok` mask: start from df1. A single non-usable scalar df1 means no fit.
    let mut ok = vec![false; n];
    if df1_scalar {
        if !(df1[0].is_finite() && df1[0] > 1e-15) {
            return Ok((Array1::from_elem(n, f64::NAN), f64::NAN));
        }
        ok.iter_mut().for_each(|b| *b = true);
    } else {
        for (i, b) in ok.iter_mut().enumerate() {
            *b = df1[i].is_finite() && df1[i] > 1e-15;
        }
    }

    // Covariate: NA disallowed; clamp ±Inf to the finite range ±1 (or to its
    // sign if every value is infinite).
    if covariate.iter().any(|v| v.is_nan()) {
        bail!("NA covariate values not allowed");
    }
    let mut cov = covariate.to_vec();
    let finite: Vec<f64> = cov.iter().copied().filter(|v| v.is_finite()).collect();
    if finite.len() < n {
        if finite.is_empty() {
            for v in cov.iter_mut() {
                *v = v.signum();
            }
        } else {
            let rmin = finite.iter().copied().fold(f64::INFINITY, f64::min);
            let rmax = finite.iter().copied().fold(f64::NEG_INFINITY, f64::max);
            for v in cov.iter_mut() {
                if *v == f64::NEG_INFINITY {
                    *v = rmin - 1.0;
                } else if *v == f64::INFINITY {
                    *v = rmax + 1.0;
                }
            }
        }
    }

    // Finalise `ok` with the variance checks.
    for (i, b) in ok.iter_mut().enumerate() {
        *b = *b && x[i].is_finite() && x[i] > -1e-15;
    }
    let ok_idx: Vec<usize> = (0..n).filter(|&i| ok[i]).collect();
    let nok = ok_idx.len();
    if nok == 0 {
        return Ok((Array1::from_elem(n, f64::NAN), f64::NAN));
    }
    if nok == 1 {
        return Ok((Array1::from_elem(n, x[ok_idx[0]]), 0.0));
    }
    let notallok = nok < n;
    let notok_idx: Vec<usize> = (0..n).filter(|&i| !ok[i]).collect();
    let cov_ok: Vec<f64> = ok_idx.iter().map(|&i| cov[i]).collect();

    // Spline d.f. from the *usable* count, capped at the number of distinct
    // usable covariate values.
    let n_unique = {
        let mut c = cov_ok.clone();
        c.sort_by(|a, b| a.partial_cmp(b).unwrap());
        c.dedup();
        c.len()
    };
    let mut splinedf = 1usize + (nok >= 3) as usize + (nok >= 6) as usize + (nok >= 30) as usize;
    splinedf = splinedf.min(n_unique);
    if splinedf < 2 {
        // Degenerate covariate: fall back to the constant-mean fit on the usable
        // subset, broadcast to all genes.
        let x_ok: Array1<f64> = ok_idx.iter().map(|&i| x[i]).collect();
        let df1_ok: Array1<f64> = if df1_scalar {
            Array1::from_elem(1, df1[0])
        } else {
            ok_idx.iter().map(|&i| df1[i]).collect()
        };
        let (scale, df2) = fit_fdist(&x_ok, &df1_ok);
        return Ok((Array1::from_elem(n, scale), df2));
    }

    // Offset exact/near-zero variances away from zero (as in fit_fdist), on the
    // usable subset.
    let mut xs: Vec<f64> = ok_idx.iter().map(|&i| x[i].max(0.0)).collect();
    let mut m = median(&xs);
    if m == 0.0 {
        m = 1.0;
    }
    let floor = 1e-5 * m;
    for v in xs.iter_mut() {
        if *v < floor {
            *v = floor;
        }
    }

    // e = log(x) + logmdigamma(df1/2),  logmdigamma(a) = log(a) - digamma(a).
    let e: Array1<f64> = ok_idx
        .iter()
        .enumerate()
        .map(|(k, &i)| {
            let d = df1_at(i);
            xs[k].ln() - digamma(d / 2.0) + (d / 2.0).ln()
        })
        .collect();

    // Fit the natural spline on the usable subset; keep the basis so the trend
    // can be evaluated (extrapolated) at the dropped genes.
    let basis = NsBasis::fit(&cov_ok, splinedf);
    let design = basis.eval(&cov_ok); // nok x splinedf
    let (q, r) = qr_econ(&design);
    let qte = q.t().dot(&e);
    let coef = solve_upper(&r, &qte);
    let fitted_ok = design.dot(&coef); // length nok

    // evar = RSS/(nok - splinedf) = mean of the residual `effects` (R uses the
    // full-QR residual effects; for a full-rank fit these agree exactly).
    let rss: f64 = e
        .iter()
        .zip(fitted_ok.iter())
        .map(|(a, b)| (a - b) * (a - b))
        .sum();
    let mut evar = rss / (nok as f64 - splinedf as f64);
    let tri_mean: f64 = ok_idx
        .iter()
        .map(|&i| trigamma(df1_at(i) / 2.0))
        .sum::<f64>()
        / nok as f64;
    evar -= tri_mean;

    // emean over all genes: fitted on the usable subset, predicted (extrapolated)
    // on the dropped genes.
    let mut emean = Array1::<f64>::zeros(n);
    for (k, &i) in ok_idx.iter().enumerate() {
        emean[i] = fitted_ok[k];
    }
    if notallok {
        let cov_notok: Vec<f64> = notok_idx.iter().map(|&i| cov[i]).collect();
        let pred = basis.eval(&cov_notok).dot(&coef);
        for (k, &i) in notok_idx.iter().enumerate() {
            emean[i] = pred[k];
        }
    }

    if evar > 0.0 {
        let df2 = 2.0 * trigamma_inverse(evar);
        let adj = digamma(df2 / 2.0) - (df2 / 2.0).ln(); // -logmdigamma(df2/2)
        let scale = emean.mapv(|em| (em + adj).exp());
        Ok((scale, df2))
    } else {
        // df2 = Inf with a covariate present: s2.0 = exp(emean).
        let scale = emean.mapv(f64::exp);
        Ok((scale, f64::INFINITY))
    }
}

/// R `quantile(x, probs, type = 7)` (the default). `sorted` must be ascending.
pub(crate) fn quantile_type7(sorted: &[f64], probs: &[f64]) -> Vec<f64> {
    let n = sorted.len();
    probs
        .iter()
        .map(|&p| {
            if n == 1 {
                return sorted[0];
            }
            let h = (n as f64 - 1.0) * p;
            let lo = h.floor() as usize;
            let hi = (lo + 1).min(n - 1);
            sorted[lo] + (h - lo as f64) * (sorted[hi] - sorted[lo])
        })
        .collect()
}

/// R `mean(x, trim)`: drop `floor(trim*n)` order statistics from each tail.
fn trimmed_mean(xs: &[f64], trim: f64) -> f64 {
    let n = xs.len();
    let mut v = xs.to_vec();
    v.sort_by(|a, b| a.partial_cmp(b).unwrap());
    let lo = (trim * n as f64).floor() as usize;
    let hi = n - lo;
    let s = &v[lo..hi];
    s.iter().sum::<f64>() / s.len() as f64
}

/// R `rank(x)` with the default `ties.method = "average"`.
fn rank_average(xs: &[f64]) -> Vec<f64> {
    let n = xs.len();
    let mut idx: Vec<usize> = (0..n).collect();
    idx.sort_by(|&a, &b| xs[a].partial_cmp(&xs[b]).unwrap());
    let mut ranks = vec![0.0_f64; n];
    let mut i = 0;
    while i < n {
        let mut j = i;
        while j + 1 < n && xs[idx[j + 1]] == xs[idx[i]] {
            j += 1;
        }
        let avg = (i + j) as f64 / 2.0 + 1.0;
        for &k in &idx[i..=j] {
            ranks[k] = avg;
        }
        i = j + 1;
    }
    ranks
}

/// R `order(x)`: ascending, ties broken by original position (stable).
fn order_indices(xs: &[f64]) -> Vec<usize> {
    let mut idx: Vec<usize> = (0..xs.len()).collect();
    idx.sort_by(|&a, &b| xs[a].partial_cmp(&xs[b]).unwrap().then(a.cmp(&b)));
    idx
}

/// F quantile allowing `df2 = +Inf` (limit `qf(p,df1,Inf) = qchisq(p,df1)/df1`).
fn qf_maybe_inf(p: f64, df1: f64, df2: f64) -> f64 {
    if df2.is_infinite() {
        chisq_qf(p, df1) / df1
    } else {
        f_qf(p, df1, df2)
    }
}

/// F density allowing `df2 = +Inf` (limit `df(x,df1,Inf) = dchisq(x*df1,df1)*df1`).
fn df_maybe_inf(x: f64, df1: f64, df2: f64) -> f64 {
    if df2.is_infinite() {
        chisq_pdf(x * df1, df1) * df1
    } else {
        f_pdf(x, df1, df2)
    }
}

/// Mean and variance of Winsorized `log F(df1, df2)` values, computed by
/// Gauss-Legendre quadrature on the link scale `q = f/(1+f)`. Port of the
/// `winsorizedMoments` closure in `fitFDistRobustly`. `wtp = (lower, upper)` is
/// `winsor.tail.p`; the quadrature interval spans the `(wtp.0, 1-wtp.1)`
/// F-quantiles.
fn winsorized_moments(
    df1: f64,
    df2: f64,
    wtp: (f64, f64),
    gnodes: &[f64],
    gweights: &[f64],
) -> (f64, f64) {
    let fq0 = qf_maybe_inf(wtp.0, df1, df2);
    let fq1 = qf_maybe_inf(1.0 - wtp.1, df1, df2);
    let zq0 = fq0.ln();
    let zq1 = fq1.ln();
    let q0 = fq0 / (1.0 + fq0);
    let q1 = fq1 / (1.0 + fq1);
    let q21 = q1 - q0;

    let nq = gnodes.len();
    let mut znodes = vec![0.0_f64; nq];
    let mut fvals = vec![0.0_f64; nq];
    let mut sum_fz = 0.0;
    for k in 0..nq {
        let node = q0 + q21 * gnodes[k];
        let fnode = node / (1.0 - node);
        let znode = fnode.ln();
        let fval = df_maybe_inf(fnode, df1, df2) / ((1.0 - node) * (1.0 - node));
        znodes[k] = znode;
        fvals[k] = fval;
        sum_fz += gweights[k] * fval * znode;
    }
    let m = q21 * sum_fz + (zq0 * wtp.0 + zq1 * wtp.1);

    let mut sum_fv = 0.0;
    for k in 0..nq {
        let d = znodes[k] - m;
        sum_fv += gweights[k] * fvals[k] * d * d;
    }
    let v = q21 * sum_fv + ((zq0 - m).powi(2) * wtp.0 + (zq1 - m).powi(2) * wtp.1);
    (m, v)
}

/// Brent's method root finder, a faithful port of R's `R_zeroin2` (`zeroin.c`).
/// Caller supplies the bracketing endpoints and their function values.
fn brent_root<F: Fn(f64) -> f64>(f: &F, ax: f64, bx: f64, fa0: f64, fb0: f64, tol: f64) -> f64 {
    let eps = f64::EPSILON;
    let (mut a, mut b) = (ax, bx);
    let (mut fa, mut fb) = (fa0, fb0);
    let mut c = a;
    let mut fc = fa;
    if fa == 0.0 {
        return a;
    }
    if fb == 0.0 {
        return b;
    }
    for _ in 0..1000 {
        let prev_step = b - a;
        if fc.abs() < fb.abs() {
            a = b;
            b = c;
            c = a;
            fa = fb;
            fb = fc;
            fc = fa;
        }
        let tol_act = 2.0 * eps * b.abs() + tol / 2.0;
        let mut new_step = (c - b) / 2.0;
        if new_step.abs() <= tol_act || fb == 0.0 {
            return b;
        }
        if prev_step.abs() >= tol_act && fa.abs() > fb.abs() {
            let cb = c - b;
            let (mut p, mut q);
            if a == c {
                let t1 = fb / fa;
                p = cb * t1;
                q = 1.0 - t1;
            } else {
                let mut qq = fa / fc;
                let t1 = fb / fc;
                let t2 = fb / fa;
                p = t2 * (cb * qq * (qq - t1) - (b - a) * (t1 - 1.0));
                qq = (qq - 1.0) * (t1 - 1.0) * (t2 - 1.0);
                q = qq;
            }
            if p > 0.0 {
                q = -q;
            } else {
                p = -p;
            }
            if p < 0.75 * cb * q - (tol_act * q).abs() / 2.0 && p < (prev_step * q / 2.0).abs() {
                new_step = p / q;
            }
        }
        if new_step.abs() < tol_act {
            new_step = if new_step > 0.0 { tol_act } else { -tol_act };
        }
        a = b;
        fa = fb;
        b += new_step;
        fb = f(b);
        if (fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0) {
            c = a;
            fc = fa;
        }
    }
    b
}

/// One-dimensional function minimizer: a port of R's `Brent_fmin`, the engine
/// behind `optimize()`. Combines golden-section search with successive
/// parabolic interpolation and returns the minimizing argument in `[ax, bx]`.
/// `tol` is the absolute tolerance on the argument; R's `optimize` default is
/// `.Machine$double.eps^0.25`.
fn brent_fmin<F: Fn(f64) -> f64>(ax: f64, bx: f64, f: F, tol: f64) -> f64 {
    // c is the squared inverse of the golden ratio.
    let c = 0.5 * (3.0 - 5.0_f64.sqrt());
    // eps is approximately the square root of the relative machine precision.
    let eps = f64::EPSILON.sqrt();

    let (mut a, mut b) = (ax, bx);
    let mut v = a + c * (b - a);
    let (mut w, mut x) = (v, v);
    let mut d = 0.0_f64;
    let mut e = 0.0_f64;
    let mut fx = f(x);
    let (mut fv, mut fw) = (fx, fx);
    let tol3 = tol / 3.0;

    loop {
        let xm = 0.5 * (a + b);
        let tol1 = eps * x.abs() + tol3;
        let t2 = 2.0 * tol1;
        // Stopping criterion.
        if (x - xm).abs() <= t2 - 0.5 * (b - a) {
            break;
        }
        let (mut p, mut q, mut r) = (0.0, 0.0, 0.0);
        if e.abs() > tol1 {
            // Fit a parabola.
            r = (x - w) * (fx - fv);
            q = (x - v) * (fx - fw);
            p = (x - v) * q - (x - w) * r;
            q = 2.0 * (q - r);
            if q > 0.0 {
                p = -p;
            } else {
                q = -q;
            }
            r = e;
            e = d;
        }

        if p.abs() >= (0.5 * q * r).abs() || p <= q * (a - x) || p >= q * (b - x) {
            // A golden-section step.
            e = if x < xm { b - x } else { a - x };
            d = c * e;
        } else {
            // A parabolic-interpolation step.
            d = p / q;
            let u = x + d;
            // f must not be evaluated too close to ax or bx.
            if u - a < t2 || b - u < t2 {
                d = if x < xm { tol1 } else { -tol1 };
            }
        }

        // f must not be evaluated too close to x.
        let u = if d.abs() >= tol1 {
            x + d
        } else if d > 0.0 {
            x + tol1
        } else {
            x - tol1
        };
        let fu = f(u);

        // Update a, b, v, w, and x.
        if fu <= fx {
            if u < x {
                b = x;
            } else {
                a = x;
            }
            v = w;
            fv = fw;
            w = x;
            fw = fx;
            x = u;
            fx = fu;
        } else {
            if u < x {
                a = u;
            } else {
                b = u;
            }
            if fu <= fw || w == x {
                v = w;
                fv = fw;
                w = u;
                fw = fu;
            } else if fu <= fv || v == x || v == w {
                v = u;
                fv = fu;
            }
        }
    }
    x
}

/// Non-robust scaled-F fit, dispatching on the covariate: `fitFDist(x, df1)`
/// (constant scale, broadcast to length `n`) or `fitFDist(x, df1, covariate)`
/// (a per-gene spline trend). Used inside `fitFDistRobustly` for its non-robust
/// reference estimate.
fn non_robust_fit(x: &[f64], df1: f64, covariate: Option<&[f64]>) -> Result<(Vec<f64>, f64)> {
    let n = x.len();
    match covariate {
        None => {
            let (s, d) = fit_fdist(&Array1::from(x.to_vec()), &Array1::from(vec![df1; n]));
            Ok((vec![s; n], d))
        }
        Some(cov) => {
            let (s_vec, d) = fit_fdist_trend(
                &Array1::from(x.to_vec()),
                &Array1::from(vec![df1; n]),
                &Array1::from(cov.to_vec()),
            )?;
            Ok((s_vec.to_vec(), d))
        }
    }
}

/// Hyperparameters of a scaled F-distribution from [`fit_fdist_unequal_df1`].
pub(crate) struct FDistUnequalDF1 {
    /// `fit$scale`: the prior variance scale `s0^2` (a scalar in the no-covariate
    /// core path).
    pub(crate) scale: f64,
    /// `fit$df2`: the prior degrees of freedom `2*d2`.
    pub(crate) df2: f64,
}

/// Shared preprocessing for the [`fit_fdist_unequal_df1`] family: applies
/// limma's NA-`x`, small-`df1`, and informative-count rules and returns the
/// per-observation quantities (`d1`, `xpos`, `e`, `w`) consumed by the
/// optimisation, together with the (possibly zeroed) prior weights.
enum UnequalPrep {
    /// Fewer than two informative observations: the scaled-F fit is undefined.
    Insufficient,
    Ready {
        n: usize,
        d1: Vec<f64>,
        xpos: Vec<f64>,
        e: Vec<f64>,
        w: Vec<f64>,
        pw: Option<Vec<f64>>,
        /// `n.informative == 2` forces limma's scalar, unweighted fit (covariate
        /// and prior weights dropped).
        drop_covariate: bool,
    },
}

/// Port of the shared head of `fitFDistUnequalDF1` (lines up to the `emean`
/// branch): validate `df1`/`prior.weights`, zero out the weight of NA `x` and
/// near-zero `df1`, count informative observations, and build `d1`, `xpos`,
/// `e = log(xpos) + logmdigamma(d1)` and `w = prior.weights / trigamma(d1)`.
fn unequal_df1_prep(x: &[f64], df1: &[f64], prior_weights: Option<&[f64]>) -> Result<UnequalPrep> {
    let n = x.len();
    if df1.len() != 1 && df1.len() != n {
        bail!("x and df1 are different lengths");
    }
    if df1.iter().any(|v| v.is_nan()) {
        bail!("NA df1 values");
    }
    let df1_at = |i: usize| if df1.len() == 1 { df1[0] } else { df1[i] };

    if let Some(pw) = prior_weights {
        if pw.len() != n {
            bail!("x and prior.weights are different lengths");
        }
        let m = pw.iter().copied().fold(f64::INFINITY, f64::min);
        if m.is_nan() {
            bail!("prior.weights contain NA values");
        }
        if m < 0.0 {
            bail!("prior.weights are negative");
        }
    }

    // Local mutable copies: NA x values and tiny df1 are made un-informative.
    let mut xv: Vec<f64> = x.to_vec();
    let mut df1f: Vec<f64> = (0..n).map(df1_at).collect();
    let mut pw: Option<Vec<f64>> = prior_weights.map(<[f64]>::to_vec);

    // NAs in x -> zero prior weight, x set to 0.
    if xv.iter().any(|v| v.is_nan()) {
        let na: Vec<bool> = xv.iter().map(|v| v.is_nan()).collect();
        match pw {
            None => pw = Some(na.iter().map(|&m| f64::from(!m)).collect()),
            Some(ref mut p) => {
                for (pi, &m) in p.iter_mut().zip(&na) {
                    if m {
                        *pi = 0.0;
                    }
                }
            }
        }
        for (xi, &m) in xv.iter_mut().zip(&na) {
            if m {
                *xi = 0.0;
            }
        }
    }

    // Treat very small df1 (< 0.01) as un-informative.
    if df1f.iter().copied().fold(f64::INFINITY, f64::min) < 0.01 {
        let small: Vec<bool> = df1f.iter().map(|&d| d < 0.01).collect();
        match pw {
            None => pw = Some(small.iter().map(|&m| f64::from(!m)).collect()),
            Some(ref mut p) => {
                for (pi, &m) in p.iter_mut().zip(&small) {
                    if m {
                        *pi = 0.0;
                    }
                }
            }
        }
        for (di, &m) in df1f.iter_mut().zip(&small) {
            if m {
                *di = 1.0;
            }
        }
    }

    // Need at least two informative observations.
    let informative: Vec<bool> = (0..n)
        .map(|j| xv[j] > 0.0 && pw.as_ref().is_none_or(|p| p[j] != 0.0))
        .collect();
    let n_informative = informative.iter().filter(|&&b| b).count();
    if n_informative < 2 {
        return Ok(UnequalPrep::Insufficient);
    }
    // With exactly two informative values, limma drops the covariate and prior
    // weights and reverts to the scalar fit.
    let drop_covariate = n_informative == 2;
    if drop_covariate {
        pw = None;
    }

    // Moment estimation on log(F), adjusted for d1. Avoid exact zeros.
    let xi: Vec<f64> = (0..n).filter(|&j| informative[j]).map(|j| xv[j]).collect();
    let med = median(&xi);
    let xpos: Vec<f64> = xv.iter().map(|&v| v.max(1e-12 * med)).collect();
    let d1: Vec<f64> = df1f.iter().map(|&d| d / 2.0).collect();
    let e: Vec<f64> = (0..n).map(|j| xpos[j].ln() + logmdigamma(d1[j])).collect();
    let mut w: Vec<f64> = d1.iter().map(|&d| 1.0 / trigamma(d)).collect();
    if let Some(ref p) = pw {
        for (wj, &pj) in w.iter_mut().zip(p) {
            *wj *= pj;
        }
    }

    Ok(UnequalPrep::Ready {
        n,
        d1,
        xpos,
        e,
        w,
        pw,
        drop_covariate,
    })
}

/// Maximise the scaled-F log-likelihood over `par = d2/(1+d2)` in `[1/2, 0.9998]`
/// and return the optimal `d2`. `emean` is the location of `log(F)`: a constant
/// vector in the core path, the per-gene `loessFit` trend in the covariate path.
/// `d2s20`/`log(d2s20)` are recomputed per observation (necessary when `emean`
/// varies); with a constant `emean` this is bit-identical to hoisting them out.
fn unequal_df1_optimize(d1: &[f64], xpos: &[f64], emean: &[f64], pw: Option<&[f64]>) -> f64 {
    let n = d1.len();
    let d1x: Vec<f64> = (0..n).map(|j| d1[j] * xpos[j]).collect();
    let minus_twice_loglik = |par: f64| -> f64 {
        let d2 = par / (1.0 - par);
        let lmd_d2 = logmdigamma(d2);
        let lgamma_d2 = ln_gamma(d2);
        let mut acc = 0.0;
        for j in 0..n {
            let d2s20 = d2 * (emean[j] - lmd_d2).exp();
            let term = -(d1[j] + d2) * (d1x[j] / d2s20).ln_1p() - d1[j] * d2s20.ln()
                + ln_gamma(d1[j] + d2)
                - lgamma_d2;
            acc += match pw {
                Some(p) => p[j] * term,
                None => term,
            };
        }
        -2.0 * acc
    };
    let tol = f64::EPSILON.powf(0.25);
    let par = brent_fmin(0.5, 0.9998, minus_twice_loglik, tol);
    par / (1.0 - par)
}

/// Maximum-likelihood scaled-F fit that allows the first degrees of freedom
/// `df1` to vary across observations. Port of limma's `fitFDistUnequalDF1`
/// (core path: `covariate = NULL`, `robust = FALSE`), selected by `squeezeVar`
/// whenever residual df differ between genes (e.g. data with NAs). Unlike the
/// moment-based [`fit_fdist`], it maximises the exact scaled-F log-likelihood
/// over `df2` with a one-dimensional optimisation, matching R's `optimize`.
///
/// `df1` may be a unit slice (broadcast to all observations) or length `n`.
/// `prior_weights`, when given, must be non-negative and length `n`. The trended
/// (covariate) variant is [`fit_fdist_unequal_df1_trend`]; the `robust` sub-path
/// is [`fit_fdist_unequal_df1_robust`].
pub(crate) fn fit_fdist_unequal_df1(
    x: &[f64],
    df1: &[f64],
    prior_weights: Option<&[f64]>,
) -> Result<FDistUnequalDF1> {
    match unequal_df1_prep(x, df1, prior_weights)? {
        UnequalPrep::Insufficient => Ok(FDistUnequalDF1 {
            scale: f64::NAN,
            df2: f64::NAN,
        }),
        UnequalPrep::Ready {
            n,
            d1,
            xpos,
            e,
            w,
            pw,
            ..
        } => {
            let sw: f64 = w.iter().sum();
            let emean: f64 = (0..n).map(|j| w[j] * e[j]).sum::<f64>() / sw;
            let emean_vec = vec![emean; n];
            let d2 = unequal_df1_optimize(&d1, &xpos, &emean_vec, pw.as_deref());
            let s20 = (emean - logmdigamma(d2)).exp();
            Ok(FDistUnequalDF1 {
                scale: s20,
                df2: 2.0 * d2,
            })
        }
    }
}

/// Trended (covariate) maximum-likelihood scaled-F fit allowing unequal `df1`.
/// Port of the `covariate != NULL`, `robust = FALSE` path of
/// `fitFDistUnequalDF1`: the location of `log(F)` is a weighted `loessFit` trend
/// in `covariate` rather than a single mean, so the returned `scale` is per-gene
/// (length `n`). The prior d.f. `df2 = 2*d2` stays scalar.
///
/// `span` defaults to `chooseLowessSpan(n, small.n = 500)` when `None`. With
/// only two informative observations limma drops the covariate, so the fit
/// reduces to the constant mean broadcast across genes.
pub(crate) fn fit_fdist_unequal_df1_trend(
    x: &[f64],
    df1: &[f64],
    covariate: &[f64],
    span: Option<f64>,
    prior_weights: Option<&[f64]>,
) -> Result<(Vec<f64>, f64)> {
    let n = x.len();
    if covariate.len() != n {
        bail!("x and covariate are different lengths");
    }
    if covariate.iter().any(|v| v.is_nan()) {
        bail!("NA covariate values not allowed");
    }
    match unequal_df1_prep(x, df1, prior_weights)? {
        UnequalPrep::Insufficient => Ok((vec![f64::NAN; n], f64::NAN)),
        UnequalPrep::Ready {
            d1,
            xpos,
            e,
            w,
            pw,
            drop_covariate,
            ..
        } => {
            // emean is the loess trend, or the scalar mean when limma dropped
            // the covariate (n.informative == 2).
            let emean: Vec<f64> = if drop_covariate {
                let sw: f64 = w.iter().sum();
                let m: f64 = (0..n).map(|j| w[j] * e[j]).sum::<f64>() / sw;
                vec![m; n]
            } else {
                let span = span.unwrap_or_else(|| choose_lowess_span(n, 500.0, 0.3, 1.0 / 3.0));
                let mut ws = w.clone();
                ws.sort_by(|a, b| a.partial_cmp(b).unwrap());
                let q75 = quantile_type7(&ws, &[0.75])[0];
                let wq: Vec<f64> = w.iter().map(|&wj| wj / q75).collect();
                loess_fit(&e, covariate, Some(&wq), span, 1, 1e-8, 1e2, true).fitted
            };
            let d2 = unequal_df1_optimize(&d1, &xpos, &emean, pw.as_deref());
            let scale: Vec<f64> = emean.iter().map(|&m| (m - logmdigamma(d2)).exp()).collect();
            Ok((scale, 2.0 * d2))
        }
    }
}

/// Dispatch a non-robust `fitFDistUnequalDF1` fit and return the per-gene scale
/// (the scalar `scale` broadcast to length `n` in the no-covariate path) with the
/// scalar `df2`. Used by the robust path for both its initial fit and the `Recall`
/// with FDR prior weights. `span` is ignored without a covariate, as in R.
fn unequal_nonrobust(
    x: &[f64],
    df1: &[f64],
    covariate: Option<&[f64]>,
    span: Option<f64>,
    prior_weights: Option<&[f64]>,
) -> Result<(Vec<f64>, f64)> {
    match covariate {
        None => {
            let r = fit_fdist_unequal_df1(x, df1, prior_weights)?;
            Ok((vec![r.scale; x.len()], r.df2))
        }
        Some(cov) => fit_fdist_unequal_df1_trend(x, df1, cov, span, prior_weights),
    }
}

/// Robust maximum-likelihood scaled-F fit allowing unequal `df1`. Port of the
/// `robust = TRUE` path of `fitFDistUnequalDF1`: an initial non-robust fit
/// identifies two-sided outliers by FDR, the distribution is refit with the FDR
/// as prior weights (R's `Recall`), and a qq-plot argument shrinks `df2` per gene
/// for the right-tail outliers (`df2.shrunk`). The returned `scale` is broadcast
/// to length `n` (no covariate) or per-gene (trended); `df2_shrunk` is the
/// per-gene prior d.f. that `squeezeVar` forwards as `df.prior`.
pub(crate) fn fit_fdist_unequal_df1_robust(
    x: &[f64],
    df1: &[f64],
    covariate: Option<&[f64]>,
    span: Option<f64>,
) -> Result<RobustFDist> {
    let n = x.len();
    if df1.len() != 1 && df1.len() != n {
        bail!("x and df1 are different lengths");
    }
    let df1_at = |i: usize| if df1.len() == 1 { df1[0] } else { df1[i] };

    // Clean copies matching the in-place edits at the head of fitFDistUnequalDF1:
    // NA x -> 0, tiny df1 (< 0.01) -> 1. The robust block and the Recall both
    // operate on these cleaned values.
    let xc: Vec<f64> = x
        .iter()
        .map(|&v| if v.is_nan() { 0.0 } else { v })
        .collect();
    let df1m: Vec<f64> = (0..n)
        .map(|i| {
            let d = df1_at(i);
            if d < 0.01 {
                1.0
            } else {
                d
            }
        })
        .collect();

    // Initial non-robust fit. Passing the ORIGINAL x/df1 is fine: the shared prep
    // applies the same NA/tiny-df1 cleaning internally, so s20/df2 match R's
    // values computed after the in-place edits.
    let (s20, df2) = unequal_nonrobust(x, df1, covariate, span, None)?;

    // R forces robust = FALSE (skipping this whole block) once fewer than three
    // observations are informative: < 2 returns NA, == 2 reverts to the scalar
    // non-robust fit. Either way the answer is exactly the initial fit.
    let n_informative = (0..n).filter(|&j| xc[j] > 0.0 && df1_at(j) >= 0.01).count();
    if n_informative < 3 {
        return Ok(RobustFDist {
            scale: s20,
            df2_shrunk: vec![df2; n],
        });
    }

    // Identify two-sided outliers via FDR. FStat uses the cleaned x and the
    // initial scale; the F tail probabilities use the modified df1.
    let fstat: Vec<f64> = (0..n).map(|j| xc[j] / s20[j]).collect();
    let right_p: Vec<f64> = (0..n).map(|j| f_sf(fstat[j], df1m[j], df2)).collect();
    let mut left_p: Vec<f64> = right_p.iter().map(|&r| 1.0 - r).collect();
    if left_p.iter().copied().fold(f64::INFINITY, f64::min) < 0.001 {
        for j in 0..n {
            if left_p[j] < 0.001 {
                // pf(FStat, df1, df2, lower.tail = TRUE): F lower-tail CDF,
                // evaluated through the incomplete beta to stay accurate.
                let y = df1m[j] * fstat[j] / (df2 + df1m[j] * fstat[j]);
                left_p[j] = betai(y, df1m[j] / 2.0, df2 / 2.0);
            }
        }
    }
    let two_sided_p: Vec<f64> = (0..n).map(|j| 2.0 * left_p[j].min(right_p[j])).collect();
    let mut fdr = p_adjust_bh(&two_sided_p);
    for f in &mut fdr {
        if *f > 0.3 {
            *f = 1.0;
        }
    }

    // No outliers -> non-robust estimates.
    if fdr.iter().copied().fold(f64::INFINITY, f64::min) == 1.0 {
        return Ok(RobustFDist {
            scale: s20,
            df2_shrunk: vec![df2; n],
        });
    }

    // Refit with FDR as prior weights (R's Recall). Pass the CLEANED xc and
    // MODIFIED df1m so the recursive prep does not re-zero genes that survived as
    // informative; span is dropped (R omits it), reverting to chooseLowessSpan.
    let (s20_pw, df2_pw) = unequal_nonrobust(&xc, &df1m, covariate, None, Some(&fdr))?;

    // qq-plot method: probability of not being a right outlier.
    let r = rank_average(&fstat);
    let prob_not_outlier: Vec<f64> = (0..n)
        .map(|j| {
            let uniform_p = (n as f64 - r[j] + 0.5) / n as f64;
            (right_p[j] / uniform_p).min(1.0)
        })
        .collect();

    // No right outliers -> robust estimates without df2 shrinkage.
    if prob_not_outlier
        .iter()
        .copied()
        .fold(f64::INFINITY, f64::min)
        == 1.0
    {
        return Ok(RobustFDist {
            scale: s20_pw,
            df2_shrunk: vec![df2_pw; n],
        });
    }

    // Posterior df for the right outliers. df2.outlier is derived from the Recall
    // df2 but evaluated at the INITIAL FStat / RightP.
    let o = order_indices(&right_p);
    let i = o[0]; // which.min(RightP)
    let min_right_p = right_p[i];
    let mut df2_shrunk: Vec<f64> = if min_right_p == 0.0 {
        prob_not_outlier.iter().map(|&p| p * df2_pw).collect()
    } else {
        let mut df2_outlier = (0.5_f64.ln() / min_right_p.ln()) * df2_pw;
        // Iterate once for accuracy using the log right-tail probability.
        let new_log_right_p = f_lsf(fstat[i], df1m[i], df2_outlier);
        df2_outlier *= 0.5_f64.ln() / new_log_right_p;
        (0..n)
            .map(|j| prob_not_outlier[j] * df2_pw + (1.0 - prob_not_outlier[j]) * df2_outlier)
            .collect()
    };

    // Force df2.shrunk to be monotonic in RightP.
    let mut df2_ordered: Vec<f64> = o.iter().map(|&k| df2_shrunk[k]).collect();
    let mut cs = 0.0;
    let mut m = vec![0.0; n];
    for k in 0..n {
        cs += df2_ordered[k];
        m[k] = cs / (k as f64 + 1.0);
    }
    let imin = order_indices(&m)[0]; // which.min(m)
    for slot in df2_ordered.iter_mut().take(imin + 1) {
        *slot = m[imin];
    }
    let mut run = f64::NEG_INFINITY;
    for slot in &mut df2_ordered {
        run = run.max(*slot);
        *slot = run;
    }
    for (pos, &k) in o.iter().enumerate() {
        df2_shrunk[k] = df2_ordered[pos];
    }

    Ok(RobustFDist {
        scale: s20_pw,
        df2_shrunk,
    })
}

/// Robust moment estimation of a scaled F-distribution. Port of
/// `fitFDistRobustly`. `scale` is the prior variance scale (`fit$scale`):
/// per-gene with a covariate (the trended path) and otherwise constant.
/// `df2_shrunk` is the per-observation posterior second d.f. (`fit$df2.shrunk`)
/// consumed by robust empirical Bayes as the prior d.f. The scalar global `df2`
/// is an intermediate that `squeezeVar` does not forward, so it is not retained.
pub(crate) struct RobustFDist {
    pub(crate) scale: Vec<f64>,
    pub(crate) df2_shrunk: Vec<f64>,
}

/// Untrended robust fit (`covariate = NULL`). The no-covariate path cannot
/// fail, so the [`fit_fdist_robustly_covariate`] `Result` is unwrapped.
pub(crate) fn fit_fdist_robustly(x: &[f64], df1: f64) -> RobustFDist {
    fit_fdist_robustly_covariate(x, df1, None)
        .expect("fitFDistRobustly without a covariate cannot fail")
}

/// Robust scaled-F moment estimation with optional intensity trend. Port of
/// `fitFDistRobustly` for scalar `df1` (the eBayes case: residual d.f. are
/// equal, so the per-gene `df1` of the R source collapses to a scalar and its
/// "transform x to constant df1" block is a no-op). With `covariate = Some(c)`
/// the log-variance trend is fitted by `loessFit(z, c, span = 0.4)`, making the
/// scale per-gene (`trend=TRUE, robust=TRUE`); without it the trend is a
/// Winsorized mean and the scale is constant. The non-usable-gene (`notallok`)
/// branch is not ported, so the covariate path errors unless every gene is
/// usable (this surfaces through the non-robust trended fit it calls first).
fn fit_fdist_robustly_covariate(
    x: &[f64],
    df1: f64,
    covariate: Option<&[f64]>,
) -> Result<RobustFDist> {
    let n = x.len();
    if n < 2 {
        let s = if n == 1 { x[0] } else { f64::NAN };
        let d = if n == 1 { 0.0 } else { f64::NAN };
        return Ok(RobustFDist {
            scale: vec![s; n],
            df2_shrunk: vec![d; n],
        });
    }
    if n == 2 {
        let (scale, df2) = non_robust_fit(x, df1, covariate)?;
        return Ok(RobustFDist {
            scale,
            df2_shrunk: vec![df2; n],
        });
    }

    let wtp = (0.05_f64, 0.1_f64);
    let prob = (0.05_f64, 0.9_f64); // (winsor.tail.p[1], 1 - winsor.tail.p[2])

    // Offset tiny/zero variances away from zero (R uses m*1e-12).
    let mut xs = x.to_vec();
    let m = median(&xs);
    let floor = m * 1e-12;
    for v in xs.iter_mut() {
        if *v < floor {
            *v = floor;
        }
    }

    // Non-robust estimate (a lower bound on df2; its scale is returned directly
    // in a couple of degenerate branches). With a covariate this is the spline
    // trend, so `nr_scale` is per-gene.
    let (nr_scale, nr_df2) = non_robust_fit(&xs, df1, covariate)?;

    // Too few observations to Winsorize: fall back to the non-robust fit.
    if wtp.0 < 1.0 / n as f64 && wtp.1 < 1.0 / n as f64 {
        return Ok(RobustFDist {
            scale: nr_scale,
            df2_shrunk: vec![nr_df2; n],
        });
    }

    let z: Vec<f64> = xs.iter().map(|v| v.ln()).collect();
    // Demean (Winsorized mean) or detrend (loess on the covariate). `ztrend` is
    // per-gene in both cases; constant when there is no covariate.
    let ztrend: Vec<f64> = match covariate {
        None => vec![trimmed_mean(&z, wtp.1); n],
        Some(cov) => loess_fit_unweighted(&z, cov, 0.4, 4).0,
    };
    let zresid: Vec<f64> = z.iter().zip(ztrend.iter()).map(|(zv, t)| zv - t).collect();

    // Moments of the Winsorized residuals.
    let mut zr_sorted = zresid.clone();
    zr_sorted.sort_by(|a, b| a.partial_cmp(b).unwrap());
    let zrq = quantile_type7(&zr_sorted, &[prob.0, prob.1]);
    let zwins: Vec<f64> = zresid.iter().map(|&v| v.clamp(zrq[0], zrq[1])).collect();
    let zwmean = zwins.iter().sum::<f64>() / n as f64;
    let zwvar = zwins.iter().map(|&v| (v - zwmean).powi(2)).sum::<f64>() / (n as f64 - 1.0);

    let (gnodes, gweights) = gauss_legendre_01(128);

    // Try df2 = Inf.
    let (mom_inf_mean, mom_inf_var) =
        winsorized_moments(df1, f64::INFINITY, wtp, &gnodes, &gweights);
    let funval_inf = (zwvar / mom_inf_var).ln();

    if funval_inf <= 0.0 {
        let ztrendcorrected: Vec<f64> = ztrend.iter().map(|t| t + zwmean - mom_inf_mean).collect();
        let s20: Vec<f64> = ztrendcorrected.iter().map(|v| v.exp()).collect();
        let fstat: Vec<f64> = z
            .iter()
            .zip(ztrendcorrected.iter())
            .map(|(zv, tc)| (zv - tc).exp())
            .collect();
        let tailp: Vec<f64> = fstat.iter().map(|&fv| chi2_sf(fv * df1, df1)).collect();
        let r = rank_average(&fstat);
        let df_pooled = n as f64 * df1;
        let mut df2_shrunk = vec![f64::INFINITY; n];
        let prob_not_outlier: Vec<f64> = (0..n)
            .map(|i| {
                let etp = (n as f64 - r[i] + 0.5) / n as f64;
                (tailp[i] / etp).min(1.0)
            })
            .collect();
        if prob_not_outlier.iter().any(|&p| p < 1.0) {
            for i in 0..n {
                if prob_not_outlier[i] < 1.0 {
                    df2_shrunk[i] = prob_not_outlier[i] * df_pooled;
                }
            }
            let o = order_indices(&tailp);
            let mut run = f64::NEG_INFINITY;
            for &i in &o {
                run = run.max(df2_shrunk[i]);
                df2_shrunk[i] = run;
            }
        }
        return Ok(RobustFDist {
            scale: s20,
            df2_shrunk,
        });
    }

    // Otherwise estimate a finite df2 by matching the Winsorized variance.
    if nr_df2.is_infinite() {
        return Ok(RobustFDist {
            scale: nr_scale,
            df2_shrunk: vec![nr_df2; n],
        });
    }
    let fun = |xx: f64| -> f64 {
        let df2 = xx / (1.0 - xx);
        let (_m, v) = winsorized_moments(df1, df2, wtp, &gnodes, &gweights);
        (zwvar / v).ln()
    };
    let rbx = nr_df2 / (1.0 + nr_df2);
    let funval_low = fun(rbx);
    let df2 = if funval_low >= 0.0 {
        nr_df2
    } else {
        let root = brent_root(&fun, rbx, 1.0, funval_low, funval_inf, 1e-8);
        root / (1.0 - root)
    };

    // Bias-correct the trend at the estimated df2.
    let (mom_mean, _mom_var) = winsorized_moments(df1, df2, wtp, &gnodes, &gweights);
    let ztrendcorrected: Vec<f64> = ztrend.iter().map(|t| t + zwmean - mom_mean).collect();
    let s20: Vec<f64> = ztrendcorrected.iter().map(|v| v.exp()).collect();

    // Posterior df2 per observation, shrinking outliers toward a smaller df2.
    let fstat: Vec<f64> = z
        .iter()
        .zip(ztrendcorrected.iter())
        .map(|(zv, tc)| (zv - tc).exp())
        .collect();
    let log_tailp: Vec<f64> = fstat.iter().map(|&fv| f_lsf(fv, df1, df2)).collect();
    let r = rank_average(&fstat);
    let log_prob_not_outlier: Vec<f64> = (0..n)
        .map(|i| {
            let log_etp = (n as f64 - r[i] + 0.5).ln() - (n as f64).ln();
            (log_tailp[i] - log_etp).min(0.0)
        })
        .collect();
    let prob_not_outlier: Vec<f64> = log_prob_not_outlier.iter().map(|&v| v.exp()).collect();
    let prob_outlier: Vec<f64> = log_prob_not_outlier.iter().map(|&v| -v.exp_m1()).collect();

    let df2_shrunk = if log_prob_not_outlier.iter().any(|&v| v < 0.0) {
        let min_log_tailp = log_tailp.iter().cloned().fold(f64::INFINITY, f64::min);
        let mut shrunk: Vec<f64> = if min_log_tailp == f64::NEG_INFINITY {
            (0..n).map(|i| prob_not_outlier[i] * df2).collect()
        } else {
            let mut df2_outlier = (0.5_f64).ln() / min_log_tailp * df2;
            let max_fstat = fstat.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
            let new_log_tailp = f_lsf(max_fstat, df1, df2_outlier);
            df2_outlier *= (0.5_f64).ln() / new_log_tailp;
            (0..n)
                .map(|i| prob_not_outlier[i] * df2 + prob_outlier[i] * df2_outlier)
                .collect()
        };
        // Force df2.shrunk to be monotonic in the tail probability.
        let o = order_indices(&log_tailp);
        let mut ordered: Vec<f64> = o.iter().map(|&i| shrunk[i]).collect();
        let mut cum = 0.0;
        let mut imin = 0;
        let mut min_mean = f64::INFINITY;
        let mut means = vec![0.0_f64; n];
        for k in 0..n {
            cum += ordered[k];
            means[k] = cum / (k as f64 + 1.0);
            if means[k] < min_mean {
                min_mean = means[k];
                imin = k;
            }
        }
        let flat = means[imin];
        for item in ordered.iter_mut().take(imin + 1) {
            *item = flat;
        }
        let mut run = f64::NEG_INFINITY;
        for v in ordered.iter_mut() {
            run = run.max(*v);
            *v = run;
        }
        for (k, &i) in o.iter().enumerate() {
            shrunk[i] = ordered[k];
        }
        shrunk
    } else {
        vec![df2; n]
    };

    Ok(RobustFDist {
        scale: s20,
        df2_shrunk,
    })
}

/// Posterior variances and prior hyperparameters. Port of `squeezeVar`.
/// `var` and `df` are per-gene, length `n_genes`. `df_prior` is the per-gene
/// prior d.f.: a constant vector for `robust=false` (the scalar `fitFDist`
/// estimate) and the shrunk vector `df2.shrunk` for `robust=true`. `var_prior`
/// is the prior scale (`fit$scale`): a constant vector without a covariate, and
/// the per-gene spline trend with one (`trend=TRUE`).
pub(crate) struct Squeeze {
    pub(crate) var_post: Array1<f64>,
    pub(crate) var_prior: Array1<f64>,
    pub(crate) df_prior: Array1<f64>,
}

/// Common residual d.f. for the legacy robust path. R's `squeezeVar` picks the
/// legacy hyperparameter method only when the positive d.f. are all equal;
/// otherwise it dispatches to `fitFDistUnequalDF1` (whose core path is ported in
/// [`fit_fdist_unequal_df1`], but whose robust/covariate variants are not).
fn df1_common(df: &Array1<f64>) -> Result<f64> {
    let d0 = df[0];
    if d0 <= 0.0 || d0.is_nan() {
        bail!("robust eBayes requires positive residual degrees of freedom");
    }
    for &d in df.iter() {
        if d != d0 {
            bail!(
                "robust eBayes with unequal residual df (fitFDistUnequalDF1) is not yet supported"
            );
        }
    }
    Ok(d0)
}

/// Posterior variances given hyperparameters. Port of `.squeezeVar`
/// (`squeezeVar.R`): canonical shrinkage where `df.prior` is finite, and the
/// prior scale itself where `df.prior` is infinite. `var_prior` is the
/// (possibly per-gene) prior scale `fit$scale`.
pub(crate) fn squeeze_var_post(
    var: &Array1<f64>,
    df: &Array1<f64>,
    var_prior: &Array1<f64>,
    df_prior: &Array1<f64>,
) -> Array1<f64> {
    let n = var.len();
    let max_dfp = df_prior.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
    let mut var_post = Array1::<f64>::zeros(n);
    if max_dfp.is_finite() {
        for i in 0..n {
            var_post[i] = (df[i] * var[i] + df_prior[i] * var_prior[i]) / (df[i] + df_prior[i]);
        }
        return var_post;
    }
    // Some prior d.f. are infinite: those genes take the prior scale directly.
    var_post.assign(var_prior);
    let min_dfp = df_prior.iter().cloned().fold(f64::INFINITY, f64::min);
    if min_dfp > 1e100 {
        return var_post;
    }
    for i in 0..n {
        if df_prior[i].is_finite() {
            var_post[i] = (df[i] * var[i] + df_prior[i] * var_post[i]) / (df[i] + df_prior[i]);
        }
    }
    var_post
}

pub(crate) fn squeeze_var(
    var: &Array1<f64>,
    df: &Array1<f64>,
    covariate: Option<&Array1<f64>>,
    robust: bool,
) -> Result<Squeeze> {
    squeeze_var_legacy(var, df, covariate, robust, None)
}

/// `squeezeVar` with an explicit override of the legacy/new hyperparameter
/// dispatch. `legacy=None` reproduces R's `legacy=NULL` default (legacy iff the
/// positive residual df are all equal); `Some(false)` forces
/// `fitFDistUnequalDF1` even for equal df (matching e.g. `diffSplice`'s
/// `legacy=FALSE`), `Some(true)` forces the moment estimators.
pub(crate) fn squeeze_var_legacy(
    var: &Array1<f64>,
    df: &Array1<f64>,
    covariate: Option<&Array1<f64>>,
    robust: bool,
    legacy_override: Option<bool>,
) -> Result<Squeeze> {
    let n = var.len();
    if n == 0 {
        bail!("var is empty");
    }
    // Guard zero-df variances.
    let mut v = var.clone();
    if df.len() > 1 {
        for i in 0..n {
            if df[i] == 0.0 {
                v[i] = 0.0;
            }
        }
    }
    // Empirical Bayes has no theoretical advantage with fewer than 3 genes.
    if n < 3 {
        return Ok(Squeeze {
            var_post: v.clone(),
            var_prior: v.clone(),
            df_prior: Array1::zeros(n),
        });
    }

    // R's squeezeVar uses the legacy moment estimators only when the positive
    // residual df are all equal; otherwise it dispatches to fitFDistUnequalDF1.
    let legacy = legacy_override.unwrap_or_else(|| {
        let (mut lo, mut hi) = (f64::INFINITY, f64::NEG_INFINITY);
        for &d in df.iter() {
            if d > 0.0 {
                lo = lo.min(d);
                hi = hi.max(d);
            }
        }
        lo == hi
    });

    let (scale, df_prior): (Array1<f64>, Array1<f64>) = if legacy {
        match (covariate, robust) {
            (Some(cov), true) => {
                let df1 = df1_common(df)?;
                let rf = fit_fdist_robustly_covariate(
                    v.as_slice().unwrap(),
                    df1,
                    Some(cov.as_slice().unwrap()),
                )?;
                (Array1::from(rf.scale), Array1::from(rf.df2_shrunk))
            }
            (Some(cov), false) => {
                let (scale_vec, df2) = fit_fdist_trend(&v, df, cov)?;
                (scale_vec, Array1::from_elem(n, df2))
            }
            (None, true) => {
                let df1 = df1_common(df)?;
                let rf = fit_fdist_robustly(v.as_slice().unwrap(), df1);
                (Array1::from(rf.scale), Array1::from(rf.df2_shrunk))
            }
            (None, false) => {
                let (scale, df2) = fit_fdist(&v, df);
                (Array1::from_elem(n, scale), Array1::from_elem(n, df2))
            }
        }
    } else {
        // Residual df vary across genes -> fitFDistUnequalDF1. The non-robust
        // constant and trended paths are ported; only the robust variant bails.
        match (covariate, robust) {
            (None, false) => {
                let r = fit_fdist_unequal_df1(v.as_slice().unwrap(), df.as_slice().unwrap(), None)?;
                (Array1::from_elem(n, r.scale), Array1::from_elem(n, r.df2))
            }
            (Some(cov), false) => {
                let (scale_vec, df2) = fit_fdist_unequal_df1_trend(
                    v.as_slice().unwrap(),
                    df.as_slice().unwrap(),
                    cov.as_slice().unwrap(),
                    None,
                    None,
                )?;
                (Array1::from(scale_vec), Array1::from_elem(n, df2))
            }
            (cov, true) => {
                let rf = fit_fdist_unequal_df1_robust(
                    v.as_slice().unwrap(),
                    df.as_slice().unwrap(),
                    cov.map(|c| c.as_slice().unwrap()),
                    None,
                )?;
                (Array1::from(rf.scale), Array1::from(rf.df2_shrunk))
            }
        }
    };
    if df_prior.iter().any(|d| d.is_nan()) {
        bail!("could not estimate prior df");
    }

    let var_post = squeeze_var_post(&v, df, &scale, &df_prior);
    Ok(Squeeze {
        var_post,
        var_prior: scale,
        df_prior,
    })
}

/// Estimate the scale `v0` in a mixture of t-distributions, one column of
/// t-statistics. Port of `tmixture_vector`.
pub(crate) fn tmixture_vector(
    tstat: &[f64],
    stdev_unscaled: &[f64],
    df: &[f64],
    proportion: f64,
    v0_lim: (f64, f64),
) -> f64 {
    // Drop missing t-statistics.
    let mut t: Vec<f64> = Vec::new();
    let mut s: Vec<f64> = Vec::new();
    let mut d: Vec<f64> = Vec::new();
    for i in 0..tstat.len() {
        if !tstat[i].is_nan() {
            t.push(tstat[i]);
            s.push(stdev_unscaled[i]);
            d.push(df[i]);
        }
    }
    let ngenes = t.len();
    if ngenes == 0 {
        return f64::NAN;
    }
    let ntarget = (proportion / 2.0 * ngenes as f64).ceil() as usize;
    if ntarget < 1 {
        return f64::NAN;
    }
    let p = (ntarget as f64 / ngenes as f64).max(proportion);

    // Put all t-statistics on the common (maximum) degrees of freedom.
    let max_df = d.iter().cloned().fold(f64::MIN, f64::max);
    let mut tabs: Vec<f64> = t.iter().map(|&v| v.abs()).collect();
    for i in 0..ngenes {
        if d[i] < max_df {
            let tail_p = t_sf(tabs[i], d[i]);
            tabs[i] = t_isf(tail_p, max_df);
        }
    }

    // Top `ntarget` statistics, descending.
    let mut order: Vec<usize> = (0..ngenes).collect();
    order.sort_by(|&a, &b| tabs[b].partial_cmp(&tabs[a]).unwrap());
    let top = &order[..ntarget];

    let mut v0 = vec![0.0_f64; ntarget];
    for (k, &gi) in top.iter().enumerate() {
        let r = (k as f64) + 1.0;
        let tstat_k = tabs[gi];
        let v1 = stdev_unscaled_sq(s[gi]);
        let p0 = 2.0 * t_sf(tstat_k, max_df);
        let ptarget = ((r - 0.5) / ngenes as f64 - (1.0 - p) * p0) / p;
        if ptarget > p0 {
            let qtarget = t_isf(ptarget / 2.0, max_df);
            v0[k] = v1 * ((tstat_k / qtarget).powi(2) - 1.0);
        }
        v0[k] = v0[k].clamp(v0_lim.0, v0_lim.1);
    }
    v0.iter().sum::<f64>() / ntarget as f64
}

#[inline]
fn stdev_unscaled_sq(s: f64) -> f64 {
    s * s
}

/// Overall moderated F-statistic per gene plus its first d.f. Port of the
/// `fstat_only=True` path of `classifyTestsF` (second d.f. is +inf here, so
/// callers use the chi-squared limit for p-values).
fn fstat_only(t: &Array2<f64>, cov_coef: &Array2<f64>) -> (Array1<f64>, f64) {
    let ngenes = t.nrows();
    let ntests = t.ncols();
    if ntests == 1 {
        let f = t.column(0).mapv(|v| v * v);
        return (f, 1.0);
    }

    // cov2cor after lifting any exactly-zero variance to 1.
    let mut cov = cov_coef.clone();
    for i in 0..cov.nrows() {
        if cov[[i, i]] == 0.0 {
            cov[[i, i]] = 1.0;
        }
    }
    let cor = cov2cor(&cov);
    let (evalues, evectors) = eigh(&cor); // ascending
    let e0 = evalues[0];
    let r = evalues.iter().filter(|&&e| e / e0 > 1e-8).count();
    let sqrt_r = (r as f64).sqrt();

    let mut fstat = Array1::<f64>::zeros(ngenes);
    for g in 0..ngenes {
        let mut acc = 0.0;
        for k in 0..r {
            let scale = 1.0 / evalues[k].sqrt() / sqrt_r;
            let mut proj = 0.0;
            for j in 0..ntests {
                proj += t[[g, j]] * evectors[[j, k]];
            }
            proj *= scale;
            acc += proj * proj;
        }
        fstat[g] = acc;
    }
    (fstat, r as f64)
}

/// Empirical Bayes moderation of a fitted model, filling `t`, `p_value`,
/// `lods`, `s2_prior`, `s2_post`, `df_prior`, `df_total`, `var_prior`,
/// `f_stat` and `f_p_value`. Port of `eBayes`. `trend=true` models the prior
/// variance as a spline trend in `fit.amean` (limma's `trend=TRUE`,
/// `covariate=Amean`); `robust=true` selects the robust hyperparameter
/// estimator. The two combine: trended-robust detrends the log-variances with
/// `loessFit` and shrinks per-gene df2 for outliers.
pub fn ebayes(
    fit: &mut MArrayLM,
    proportion: f64,
    stdev_coef_lim: (f64, f64),
    trend: bool,
    robust: bool,
) -> Result<()> {
    let n_genes = fit.n_genes();
    let ncoef = fit.n_coef();

    let df_max = fit.df_residual.iter().cloned().fold(f64::MIN, f64::max);
    if df_max == 0.0 {
        bail!("no residual degrees of freedom in linear model fits");
    }
    if !fit.sigma.iter().any(|v| v.is_finite()) {
        bail!("no finite residual standard deviations");
    }

    // Trend covariate: the average log-expression (`fit$Amean`).
    let covariate = if trend {
        if fit.amean.len() != n_genes {
            bail!("need Amean component in fit to estimate trend");
        }
        Some(&fit.amean)
    } else {
        None
    };

    let var = fit.sigma.mapv(|s| s * s);
    let sq = squeeze_var(&var, &fit.df_residual, covariate, robust)?;
    let s2_prior = sq.var_prior;
    let s2_post = sq.var_post;
    let df_prior = sq.df_prior;
    // `var.prior.lim` and the `var.prior` fallback use the median prior scale;
    // for the untrended path `s2_prior` is constant, so the median is that
    // constant (matching the scalar `out$s2.prior` of `.ebayes`).
    let s2_prior_med = median(s2_prior.as_slice().unwrap());

    // Moderated t-statistics.
    let mut t = Array2::<f64>::zeros((n_genes, ncoef));
    for g in 0..n_genes {
        let denom = s2_post[g].sqrt();
        for j in 0..ncoef {
            t[[g, j]] = fit.coefficients[[g, j]] / fit.stdev_unscaled[[g, j]] / denom;
        }
    }

    // Total degrees of freedom, capped at the pooled residual d.f.
    let df_pooled: f64 = fit.df_residual.iter().filter(|v| v.is_finite()).sum();
    let mut df_total = Array1::<f64>::zeros(n_genes);
    for g in 0..n_genes {
        df_total[g] = (fit.df_residual[g] + df_prior[g]).min(df_pooled);
    }

    // Two-sided p-values on per-gene df_total.
    let mut p_value = Array2::<f64>::zeros((n_genes, ncoef));
    for g in 0..n_genes {
        for j in 0..ncoef {
            p_value[[g, j]] = t_two_sided_pvalue(t[[g, j]], df_total[g]);
        }
    }

    // B-statistic: prior variance of true log-fold-changes per coefficient.
    let v0_lim = (
        stdev_coef_lim.0.powi(2) / s2_prior_med,
        stdev_coef_lim.1.powi(2) / s2_prior_med,
    );
    let mut var_prior = Array1::<f64>::zeros(ncoef);
    for j in 0..ncoef {
        let tcol: Vec<f64> = (0..n_genes).map(|g| t[[g, j]]).collect();
        let scol: Vec<f64> = (0..n_genes).map(|g| fit.stdev_unscaled[[g, j]]).collect();
        let dcol: Vec<f64> = df_total.to_vec();
        let v = tmixture_vector(&tcol, &scol, &dcol, proportion, v0_lim);
        var_prior[j] = if v.is_nan() { 1.0 / s2_prior_med } else { v };
    }

    // Log-odds. The Inf-df kernel is selected per gene (`df.prior > 1e6`),
    // matching the per-row `Infdf` branch in `.ebayes`.
    let log_odds_const = (proportion / (1.0 - proportion)).ln();
    let mut lods = Array2::<f64>::zeros((n_genes, ncoef));
    for g in 0..n_genes {
        let inf_df = df_prior[g] > 1e6;
        for j in 0..ncoef {
            let su2 = fit.stdev_unscaled[[g, j]].powi(2);
            let r = (su2 + var_prior[j]) / su2;
            let t2 = t[[g, j]].powi(2);
            let kernel = if inf_df {
                t2 * (1.0 - 1.0 / r) / 2.0
            } else {
                (1.0 + df_total[g]) / 2.0 * ((t2 + df_total[g]) / (t2 / r + df_total[g])).ln()
            };
            lods[[g, j]] = log_odds_const - r.ln() / 2.0 + kernel;
        }
    }

    // Moderated F-statistics (design is full rank by construction in lmfit).
    // df2 = df.prior + df.residual (per gene), matching `classifyTestsF`; use
    // the F distribution, falling back to its chi-squared limit when df.prior
    // is infinite (the limiting case where df2 -> Inf).
    let (f, df1) = fstat_only(&t, &fit.cov_coefficients);
    let mut f_p_value = Array1::<f64>::zeros(n_genes);
    for g in 0..n_genes {
        let df2 = df_prior[g] + fit.df_residual[g];
        f_p_value[g] = if df2.is_finite() {
            f_sf(f[g], df1, df2)
        } else {
            chi2_sf(df1 * f[g], df1)
        };
    }

    fit.df_prior = Some(df_prior);
    fit.s2_prior = Some(s2_prior);
    fit.var_prior = Some(var_prior);
    fit.proportion = Some(proportion);
    fit.s2_post = Some(s2_post);
    fit.t = Some(t);
    fit.df_total = Some(df_total);
    fit.p_value = Some(p_value);
    fit.lods = Some(lods);
    fit.f_stat = Some(f);
    fit.f_p_value = Some(f_p_value);
    Ok(())
}

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

    /// Deterministic variance vector, identical to the R reference fixture:
    /// `x[i] = ((i*i) mod 97 + 1)/50` (1-based) with two injected outliers.
    fn fixture_x(n: usize) -> Vec<f64> {
        let mut x = vec![0.0_f64; n];
        for i in 1..=n {
            let ii = ((i * i) as i64 % 97) as f64;
            x[i - 1] = (ii + 1.0) / 50.0;
        }
        x[9] *= 30.0;
        x[19] *= 25.0;
        x
    }

    fn rel(a: f64, b: f64) -> f64 {
        ((a - b) / b).abs()
    }

    // Reference: limma 3.68.3 fitFDist(x, df1=13) on the fixture vector.
    #[test]
    fn fit_fdist_matches_r_on_fixture() {
        let x = fixture_x(50);
        let (scale, df2) = fit_fdist(&Array1::from(x), &Array1::from(vec![13.0; 50]));
        assert!(rel(scale, 0.653_342_419_609_763_7) < 1e-11, "scale={scale}");
        assert!(rel(df2, 3.372_566_440_190_4) < 1e-11, "df2={df2}");
    }

    // Reference: limma 3.68.3 fitFDistUnequalDF1(x, df1=...) core path
    // (covariate=NULL, robust=FALSE) via limma:::fitFDistUnequalDF1.
    #[test]
    fn fit_fdist_unequal_df1_matches_r() {
        let nan = f64::NAN;
        // case1: unequal df1, no NA.
        let x1 = [0.8, 1.2, 0.5, 2.1, 0.3, 1.7, 0.9, 3.4, 0.6, 1.1];
        let df1a = [3.0, 5.0, 3.0, 8.0, 4.0, 6.0, 3.0, 10.0, 5.0, 4.0];
        let r1 = fit_fdist_unequal_df1(&x1, &df1a, None).unwrap();
        assert!(
            rel(r1.scale, 1.564_862_016_179_576_5) < 1e-9,
            "scale1={}",
            r1.scale
        );
        assert!(
            rel(r1.df2, 24.131_279_146_203_116) < 1e-6,
            "df21={}",
            r1.df2
        );

        // case2: equal df1 as a unit slice -> df2 lands near the upper bound.
        let x2 = [0.7, 1.5, 0.4, 2.0, 0.9, 1.3, 0.6, 2.8, 1.1, 0.5, 1.9, 0.8];
        let r2 = fit_fdist_unequal_df1(&x2, &[4.0], None).unwrap();
        assert!(
            rel(r2.scale, 1.345_949_488_248_429) < 1e-9,
            "scale2={}",
            r2.scale
        );
        assert!(rel(r2.df2, 8_134.844_780_382_662) < 1e-4, "df22={}", r2.df2);

        // case3: NAs in x (passed as f64::NAN), unequal df1.
        let x3 = [0.8, nan, 0.5, 2.1, 0.3, 1.7, nan, 3.4, 0.6, 1.1];
        let r3 = fit_fdist_unequal_df1(&x3, &df1a, None).unwrap();
        assert!(
            rel(r3.scale, 1.581_632_198_266_799_1) < 1e-9,
            "scale3={}",
            r3.scale
        );
        assert!(rel(r3.df2, 18.285_785_769_935_58) < 1e-6, "df23={}", r3.df2);

        // case4: a tiny df1 (< 0.01) is made un-informative.
        let df1c = [3.0, 5.0, 0.005, 8.0, 4.0, 6.0, 3.0, 10.0, 5.0, 4.0];
        let r4 = fit_fdist_unequal_df1(&x1, &df1c, None).unwrap();
        assert!(
            rel(r4.scale, 1.651_633_468_036_861_2) < 1e-9,
            "scale4={}",
            r4.scale
        );
        assert!(
            rel(r4.df2, 31.759_566_849_299_393) < 1e-6,
            "df24={}",
            r4.df2
        );

        // case5: prior weights, unequal df1.
        let pw5 = [1.0, 0.5, 1.0, 0.2, 1.0, 0.8, 1.0, 0.3, 1.0, 0.9];
        let r5 = fit_fdist_unequal_df1(&x1, &df1a, Some(&pw5)).unwrap();
        assert!(
            rel(r5.scale, 1.178_112_852_205_813_3) < 1e-9,
            "scale5={}",
            r5.scale
        );
        assert!(
            rel(r5.df2, 21.974_339_432_003_557) < 1e-6,
            "df25={}",
            r5.df2
        );

        // case6: fewer than two informative values -> NA.
        let x6 = [0.0, 0.0, 0.5, 0.0, 0.0];
        let r6 = fit_fdist_unequal_df1(&x6, &[3.0], None).unwrap();
        assert!(r6.scale.is_nan() && r6.df2.is_nan());

        // case7: larger unequal-df1 sample (n=30).
        let x7 = [
            0.104_789_932_110_749_54,
            1.392_380_916_692_821_2,
            0.751_353_389_051_935_5,
            1.529_208_238_282_011_4,
            0.514_936_156_173_895_7,
            0.963_397_456_564_684_3,
            0.432_758_433_788_116_6,
            0.543_419_193_926_392_4,
            0.296_464_064_471_563_4,
            0.735_642_598_718_133_5,
            0.331_071_707_567_774_15,
            0.353_013_000_689_392_4,
            0.875_996_261_323_835_3,
            0.315_281_707_651_102_5,
            1.232_284_232_600_559_9,
            0.875_390_966_988_345_4,
            0.538_798_517_365_506_9,
            0.543_189_732_890_084_8,
            2.153_853_344_912_85,
            0.685_024_743_734_213_1,
            1.622_865_451_693_436,
            0.265_130_202_735_927_5,
            0.915_634_720_794_757_8,
            0.312_262_077_723_848_9,
            1.555_920_154_490_626_5,
            1.026_054_830_899_120_2,
            0.602_789_157_706_657_9,
            0.342_758_453_305_209_9,
            0.424_560_016_065_538_2,
            2.434_005_877_511_022_5,
        ];
        let df1f = [
            2.0, 6.0, 2.0, 10.0, 11.0, 5.0, 3.0, 11.0, 2.0, 9.0, 8.0, 5.0, 10.0, 6.0, 5.0, 11.0,
            3.0, 4.0, 10.0, 10.0, 12.0, 5.0, 6.0, 6.0, 5.0, 3.0, 9.0, 4.0, 11.0, 2.0,
        ];
        let r7 = fit_fdist_unequal_df1(&x7, &df1f, None).unwrap();
        assert!(
            rel(r7.scale, 0.806_573_237_883_552_6) < 1e-9,
            "scale7={}",
            r7.scale
        );
        assert!(rel(r7.df2, 23.622_512_955_700_71) < 1e-6, "df27={}", r7.df2);
    }

    /// Compare `fit_fdist_unequal_df1_trend` against
    /// `limma:::fitFDistUnequalDF1(x, df1, covariate, span, robust=FALSE)`.
    fn check_trend_unequal(
        tag: &str,
        x: &[f64],
        df1: &[f64],
        cov: &[f64],
        span: Option<f64>,
        want_scale: &[f64],
        want_df2: f64,
    ) {
        let (scale, df2) = fit_fdist_unequal_df1_trend(x, df1, cov, span, None).unwrap();
        assert_eq!(scale.len(), want_scale.len(), "{tag}: scale length");
        for (i, (&g, &w)) in scale.iter().zip(want_scale.iter()).enumerate() {
            assert!(rel(g, w) < 1e-7, "{tag}: scale[{i}] = {g} vs {w}");
        }
        assert!(
            rel(df2, want_df2) < 1e-6,
            "{tag}: df2 = {df2} vs {want_df2}"
        );
    }

    // Reference: limma 3.68.3 fitFDistUnequalDF1(x, df1, covariate, span,
    // robust=FALSE) — the trended (covariate) path, where the location of log(F)
    // is a weighted loessFit trend so `scale` is per-gene. Covers unequal df1,
    // NA variances, zero-df genes, the two-informative drop-covariate fallback,
    // and an explicit sub-1 span.
    #[test]
    fn fit_fdist_unequal_df1_trend_matches_r() {
        let nan = f64::NAN;

        // U1: unequal df1, finite random covariate, all x positive (span -> 1).
        let x1 = [
            0.801_984_633_686_912_8,
            0.49609033183404244,
            0.905_721_912_086_098_6,
            1.161_295_856_232_137,
            0.930_558_475_187_913_7,
            0.34282888219304597,
            0.801_803_928_440_035_7,
            0.529_970_347_187_704,
            0.303_317_739_429_659_7,
            0.43634578199806073,
            2.872_191_225_200_469,
            1.082735565271354,
            0.554_020_214_658_104_7,
            0.17232615292327957,
            1.4279526086087926,
            0.918_359_347_651_539_4,
            0.18975270382553294,
            0.914_862_033_216_705_5,
            0.500_074_201_902_390_2,
            0.162_317_160_109_331_3,
            0.34258228214197145,
            0.11489383113220689,
            0.19955236342173868,
            0.992_356_731_531_247,
            1.0878808069679275,
            0.950_625_063_398_255_6,
            0.4636445358568187,
            0.820388660004991,
        ];
        let df1a = [
            3.0, 4.0, 5.0, 6.0, 7.0, 3.0, 4.0, 5.0, 6.0, 7.0, 3.0, 4.0, 5.0, 6.0, 7.0, 3.0, 4.0,
            5.0, 6.0, 7.0, 3.0, 4.0, 5.0, 6.0, 7.0, 3.0, 4.0, 5.0,
        ];
        let cov1 = [
            -1.6381258172166453,
            -1.5208792882193276,
            -0.904_509_631_669_526_9,
            -0.852_615_222_673_045,
            -0.814_233_547_594_512_6,
            -0.6596604250247311,
            -0.513_755_332_649_142_2,
            -0.24253337560472144,
            -0.22135494681387732,
            -0.14773505705058476,
            -0.134_460_797_296_282_5,
            -0.11766305700394353,
            -0.11723053018352804,
            0.056830125075783555,
            0.07625359435606488,
            0.078_923_467_779_651_35,
            0.086_897_989_221_599_43,
            0.10964632093950549,
            0.286_428_527_955_366_8,
            0.365_531_126_467_715_8,
            0.457_786_661_261_753_7,
            0.507_460_112_746_367_5,
            1.0348596017885885,
            1.057_763_572_762_019,
            1.5277688307054003,
            1.5788041468947034,
            2.110_755_671_571_493,
            3.4887543381728574,
        ];
        let scale1 = [
            0.964_214_055_847_137,
            0.932_868_069_373_103,
            0.794_184_091_522_826_5,
            0.784_482_936_559_639_7,
            0.777_501_185_020_315_8,
            0.750_993_113_906_070_5,
            0.728_210_279_024_498_1,
            0.691_138_407_955_118_9,
            0.688_512_158_265_063_1,
            0.679_669_472_658_891_4,
            0.678_122_108_223_529,
            0.676_184_534_368_422_9,
            0.676_134_945_723_629_1,
            0.657_387_540_352_499_4,
            0.655_442_938_136_168_3,
            0.655_177_950_324_495_1,
            0.654_389_803_364_642_3,
            0.6521689938653239,
            0.636_329_577_225_599_2,
            0.630_107_746_376_022_2,
            0.6236411784772774,
            0.620_581_944_802_589_5,
            0.623_695_395_508_009_8,
            0.625_319_725_983_554_8,
            0.680_418_298_330_631_2,
            0.687_241_530_950_512_7,
            0.760_685_030_726_688_4,
            0.987_017_616_485_316_8,
        ];
        check_trend_unequal(
            "U1",
            &x1,
            &df1a,
            &cov1,
            None,
            &scale1,
            44.748_552_084_936_58,
        );

        // U2: unequal df1 with NA variances (-> prior weight 0), covariate trend.
        let x2 = [
            0.997_827_271_180_056_4,
            0.365_960_357_657_678_3,
            1.5547356748719023,
            0.47954394487772123,
            0.568_435_096_195_994_5,
            0.31232773956167464,
            nan,
            1.1193342334002578,
            0.10960672354802832,
            2.149_493_457_315_547,
            0.796_222_061_916_288_4,
            0.3767740359181404,
            0.733299842535031,
            0.41302310065282277,
            nan,
            1.2551841332286193,
            1.2272049006422245,
            1.1383466000315061,
            0.42089067100222827,
            0.499_217_334_092_653_7,
            0.255_895_942_270_455,
            1.3662665194530537,
            nan,
            0.9827968020995782,
            0.793_332_318_261_128_9,
            0.519_316_173_913_811,
            0.635_014_486_905_870_8,
            0.37800257930584635,
            0.48141233737263567,
            2.0587086786268594,
        ];
        let df1b = [
            4.0, 5.0, 6.0, 4.0, 5.0, 6.0, 4.0, 5.0, 6.0, 4.0, 5.0, 6.0, 4.0, 5.0, 6.0, 4.0, 5.0,
            6.0, 4.0, 5.0, 6.0, 4.0, 5.0, 6.0, 4.0, 5.0, 6.0, 4.0, 5.0, 6.0,
        ];
        let cov2 = [
            -2.2261639974601546,
            -1.3972219141336129,
            -1.1268572787540512,
            -1.0511252248261194,
            -0.859_770_048_937_119,
            -0.657_198_535_197_769_7,
            -0.527_231_766_897_207,
            -0.4675843559750571,
            -0.45847494387506166,
            -0.45802417478849006,
            -0.2613468347135362,
            -0.25387945051265204,
            -0.088_721_641_174_677_18,
            -0.015045940014992758,
            0.1004052332328539,
            0.13310736259823772,
            0.1718059907420289,
            0.33859292004352803,
            0.38398535622573915,
            0.395_726_956_920_001_9,
            0.518_591_553_346_557_7,
            0.829_617_944_186_455_8,
            1.0913385693953688,
            1.1008530725631553,
            1.1913745681294863,
            1.2611894175389573,
            1.422_484_083_101_07,
            1.6139706522893342,
            1.6304054614305334,
            1.8897472999943068,
        ];
        let scale2 = [
            0.824_406_966_920_449_5,
            0.760_830_977_847_104_4,
            0.745_564_993_774_771_8,
            0.7417418093741357,
            0.733_006_211_876_814_1,
            0.725_458_084_468_013_1,
            0.722_226_857_622_797_8,
            0.721_501_827_787_449_3,
            0.721_447_053_323_366_1,
            0.721444763702765,
            0.725_089_087_237_561_7,
            0.725_423_132_006_360_9,
            0.734_758_603_946_452_5,
            0.739_536_976_804_548_6,
            0.748_843_817_915_800_2,
            0.751_911_723_325_340_6,
            0.755_770_393_701_450_7,
            0.7747532870298679,
            0.780_448_170_351_698_8,
            0.781_950_743_726_128_3,
            0.7982890153801171,
            0.843_472_625_758_987_1,
            0.8848027691870185,
            0.8863512995348547,
            0.901_224_823_067_436_2,
            0.912_855_616_835_919_4,
            0.9401895491824922,
            0.973_449_200_018_210_3,
            0.976_347_380_702_072_6,
            1.0234023796361744,
        ];
        check_trend_unequal("U2", &x2, &df1b, &cov2, None, &scale2, 118.7791641360422);

        // U3: per-gene df1 with some zero (saturated) genes, covariate trend.
        let x3 = [
            0.251_040_123_762_348_4,
            0.946_029_710_379_464_7,
            0.433_370_810_895_133,
            0.19761628916607515,
            0.11599935906047122,
            1.381_026_234_129_436,
            0.15093375322566088,
            0.257_771_191_974_910_2,
            1.1222031731629063,
            0.309_412_625_655_895_7,
            1.029220776290807,
            0.39270648035090616,
            0.14995215009707963,
            0.17777239578489354,
            1.4784162003920487,
            0.483_774_358_826_323_1,
            1.121_198_478_871_973,
            1.0234309766884047,
            0.883_223_797_453_814_2,
            0.13485144365130217,
            0.635_722_498_403_878_1,
            1.939_045_565_103_84,
            1.4499385045760065,
            0.691_323_254_113_878_1,
            1.2859299074525188,
            1.0601974751143668,
        ];
        let df1c = [
            5.0, 5.0, 5.0, 0.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 0.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0,
            5.0, 0.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0,
        ];
        let cov3 = [
            -2.0,
            -1.8,
            -1.6,
            -1.4,
            -1.2,
            -1.0,
            -0.799_999_999_999_999_8,
            -0.599_999_999_999_999_9,
            -0.399_999_999_999_999_9,
            -0.19999999999999996,
            0.0,
            0.20000000000000018,
            0.40000000000000036,
            0.600_000_000_000_000_1,
            0.800_000_000_000_000_3,
            1.0,
            1.2000000000000002,
            1.4000000000000004,
            1.6,
            1.8000000000000003,
            2.0,
            2.2,
            2.4000000000000004,
            2.6000000000000005,
            2.8000000000000007,
            3.0,
        ];
        let scale3 = [
            0.40668992597816683,
            0.41242521568397245,
            0.418_940_566_349_830_1,
            0.42624787518533724,
            0.434_281_670_038_789_4,
            0.44297807062391326,
            0.45237006581909883,
            0.462_370_641_836_691_2,
            0.4729499347484919,
            0.48414136183254186,
            0.495831269947135,
            0.507_996_048_552_190_2,
            0.521_927_576_030_764_6,
            0.548_941_033_557_573_3,
            0.588_138_724_695_406_8,
            0.631_020_760_927_771,
            0.677_645_842_967_236_6,
            0.728_398_476_539_170_2,
            0.7834974061783373,
            0.843_063_743_580_241_2,
            0.907_351_793_618_404_7,
            0.977_045_671_092_314_5,
            1.052954102166904,
            1.135693922221821,
            1.2257551602924397,
            1.3237287641248294,
        ];
        check_trend_unequal("U3", &x3, &df1c, &cov3, None, &scale3, 15.901367270851837);

        // U4: only two informative x -> covariate dropped, scalar broadcast.
        let x4 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.3, 0.7];
        let df1d = [3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 5.0, 6.0];
        let cov4 = [
            -1.9557561982415352,
            -1.292_977_397_553_767,
            0.19421376806649077,
            0.694_854_618_427_691_8,
            0.8076775637491721,
            0.868_494_382_625_113_1,
            0.938_820_403_701_890_8,
            0.9604788337215433,
            1.6274587077895015,
        ];
        let scale4 = [1.884_635_517_965_289e-10; 9];
        check_trend_unequal("U4", &x4, &df1d, &cov4, None, &scale4, 2.0005878237352013);

        // U5: explicit span 0.5 -> a localized weighted-loess window.
        let x5 = [
            0.1289812457101768,
            0.3002112134889286,
            0.550_040_659_074_600_1,
            0.46573099919469607,
            1.0774952902514372,
            1.3182610807088602,
            0.16047799881484517,
            0.518_500_852_710_481_2,
            0.834_275_511_099_624_7,
            0.725_549_432_702_673_6,
            1.2140833796890174,
            1.5551077419721215,
            0.784269556528571,
            1.006_756_111_897_707,
            1.321653444985549,
            0.951_321_286_084_540_5,
            1.5304413970943935,
            0.589_800_355_783_235_8,
            1.036_073_850_772_561,
            1.0277342528581854,
            0.423_439_990_506_606_2,
            1.9646456227922258,
            0.842_085_831_926_768_4,
            1.7878805728286589,
            3.596_988_797_463_115,
            0.30296767259764334,
            1.0045659339077193,
            0.796_699_441_279_788_1,
            1.326_755_222_062_551,
            1.0695105442178703,
            2.144_881_755_828_559,
            0.622_759_395_564_611_2,
            2.8276998203872314,
            2.1903927687989184,
            2.483062167499285,
            0.749_098_603_659_653_3,
            3.907_702_645_448_939,
            0.20615431755893485,
            2.6038277981171647,
            3.3549351981825546,
        ];
        let df1e = [
            4.0, 6.0, 8.0, 4.0, 6.0, 8.0, 4.0, 6.0, 8.0, 4.0, 6.0, 8.0, 4.0, 6.0, 8.0, 4.0, 6.0,
            8.0, 4.0, 6.0, 8.0, 4.0, 6.0, 8.0, 4.0, 6.0, 8.0, 4.0, 6.0, 8.0, 4.0, 6.0, 8.0, 4.0,
            6.0, 8.0, 4.0, 6.0, 8.0, 4.0,
        ];
        let cov5 = [
            -1.5719363672101112,
            -1.4876559802763973,
            -1.205668773420278,
            -0.947_834_452_741_442_7,
            -0.867_496_844_322_951_9,
            -0.7734482255431081,
            -0.633_064_430_586_404_2,
            -0.578_397_572_427_471,
            -0.568_776_328_969_079_9,
            -0.545_470_040_900_026_2,
            -0.45238173772928986,
            -0.40228395729565475,
            -0.37820515326114934,
            -0.3192344865182718,
            -0.281446426167763,
            -0.26463622161444156,
            -0.11758490820349062,
            -0.087_536_387_129_707_54,
            -0.068_746_308_838_045_21,
            -0.0030568754017032697,
            0.012832436423186324,
            0.13596911856262384,
            0.21602846285549585,
            0.323_455_319_680_744_9,
            0.40547638875042125,
            0.432_897_640_553_402_5,
            0.633_472_674_108_249_8,
            0.6775708119851247,
            0.738_835_424_704_705_9,
            0.832_661_427_122_176_1,
            0.925_278_166_330_490_5,
            0.942_147_124_589_913_5,
            1.0805587293931278,
            1.2176667577094684,
            1.225_149_759_397_817,
            1.3975602602488622,
            1.523_537_212_272_87,
            1.672922272062954,
            1.7787585981081968,
            2.1023384239366663,
        ];
        let scale5 = [
            0.31763573095344355,
            0.35495680771975285,
            0.507_762_811_057_252_7,
            0.713_610_650_817_557_4,
            0.802_122_953_936_241_4,
            0.898_524_419_820_902,
            0.998_311_491_740_839_9,
            1.0321407899461048,
            1.0369314022341432,
            1.046_608_644_941_036,
            1.056_490_821_009_27,
            1.0608076920358882,
            1.0631124294236671,
            1.0705923543936118,
            1.0817777624846396,
            1.085_807_667_204_644,
            1.1391615373609474,
            1.1505885393961128,
            1.1555192340512663,
            1.1615575175386534,
            1.1586627709674833,
            1.1503445819382647,
            1.1458850464433143,
            1.164_109_823_638_478,
            1.189_774_858_809_396,
            1.2021255405056286,
            1.355_933_485_866_249,
            1.4042536351786707,
            1.4568194314440612,
            1.5193055045767163,
            1.5450753257905399,
            1.5412433019622132,
            1.5183893273566147,
            1.5524062081822296,
            1.5559793546353367,
            1.656464028255052,
            1.7334149016111868,
            1.8172474003613945,
            1.8699816758773846,
            2.0301683389239225,
        ];
        check_trend_unequal(
            "U5",
            &x5,
            &df1e,
            &cov5,
            Some(0.5),
            &scale5,
            8_134.844_780_382_662,
        );
    }

    /// Compare `fit_fdist_unequal_df1_robust` against
    /// `limma:::fitFDistUnequalDF1(x, df1, covariate, robust=TRUE)`: `scale` is
    /// `fit$scale` (broadcast to length `n` when scalar) and `df_prior` is
    /// `fit$df2.shrunk`, falling back to `fit$df2` broadcast on the early-return
    /// branches.
    fn check_robust_unequal(
        tag: &str,
        x: &[f64],
        df1: &[f64],
        cov: Option<&[f64]>,
        want_scale: &[f64],
        want_df_prior: &[f64],
    ) {
        let rf = fit_fdist_unequal_df1_robust(x, df1, cov, None).unwrap();
        assert_eq!(rf.scale.len(), want_scale.len(), "{tag}: scale length");
        assert_eq!(
            rf.df2_shrunk.len(),
            want_df_prior.len(),
            "{tag}: df_prior length"
        );
        for (i, (&g, &w)) in rf.scale.iter().zip(want_scale).enumerate() {
            assert!(rel(g, w) < 1e-6, "{tag}: scale[{i}] = {g} vs {w}");
        }
        for (i, (&g, &w)) in rf.df2_shrunk.iter().zip(want_df_prior).enumerate() {
            assert!(rel(g, w) < 1e-6, "{tag}: df_prior[{i}] = {g} vs {w}");
        }
    }

    // Reference: limma 3.68.3 fitFDistUnequalDF1(x, df1, covariate, robust=TRUE).
    // Exercises FDR outlier detection, the FDR-weighted Recall refit and the
    // per-gene df2.shrunk machinery: R1 mixes big right and tiny left outliers
    // (the LeftP < 0.001 refinement), R2 is clean data hitting the early df2
    // return, R3 adds a covariate trend, and R4 carries NA variances.
    #[test]
    fn fit_fdist_unequal_df1_robust_matches_r() {
        // R1: unequal df1, no covariate, big right + tiny left outliers.
        let x1 = [
            1.108_151_675_631_431,
            1.0509298052650242,
            0.600_188_499_975_963_4,
            2.965_019_643_949_733,
            0.335_247_401_948_681_1,
            38.575_940_171_816_23,
            0.661_832_021_034_727_8,
            0.729_332_580_660_132_3,
            1.370_924_562_914_53,
            0.565_333_358_945_241_7,
            0.000_891_358_205_870_448_8,
            0.226_015_781_731_191_2,
            0.995_815_680_519_283_3,
            0.21916943903107566,
            1.859_624_487_731_096,
            2.7945079795561356,
            0.657_129_526_693_096_4,
            2.8236981298779793,
            15.694211870925326,
            1.3433221526016594,
            0.47840247772910144,
            1.1840658422439412,
            1.3671832322406985,
            1.2420288840441553,
            0.599_817_535_270_839_1,
            0.803_976_948_632_128_3,
            0.0023237417261068607,
            0.972_652_819_405_507_6,
            0.929_518_912_491_085_1,
            1.0344017922302913,
            1.7457986838269346,
            1.7008009957412349,
            32.231_693_546_203_29,
            2.255368250444707,
            1.0678741460123906,
            0.840_824_740_588_304_1,
            1.1529620359637391,
            2.4040000636039216,
            0.651_835_460_845_672_7,
            0.597_775_628_599_605_4,
        ];
        let df1a = [
            4.0, 8.0, 12.0, 6.0, 10.0, 4.0, 8.0, 12.0, 6.0, 10.0, 4.0, 8.0, 12.0, 6.0, 10.0, 4.0,
            8.0, 12.0, 6.0, 10.0, 4.0, 8.0, 12.0, 6.0, 10.0, 4.0, 8.0, 12.0, 6.0, 10.0, 4.0, 8.0,
            12.0, 6.0, 10.0, 4.0, 8.0, 12.0, 6.0, 10.0,
        ];
        let scale1 = [1.0142735028598235_f64; 40];
        let mut dp1 = [8.115526155129956_f64; 40];
        for k in [5usize, 18, 32] {
            dp1[k] = 2.6473098651492584;
        }
        check_robust_unequal("R1", &x1, &df1a, None, &scale1, &dp1);

        // R2: clean data, FDR finds no outliers -> early df2 return.
        let x2 = [
            0.763_078_304_530_354_4,
            0.565_521_335_131_736_9,
            1.075985312049264,
            0.724_154_327_805_240_5,
            0.549_728_688_483_549,
            0.436_372_438_907_753_8,
            0.721_587_520_982_636,
            1.7674344836815208,
            1.1562115747707173,
            1.3091425739228408,
            1.112_939_888_776_984,
            0.538_081_098_293_272_8,
            0.550_699_286_960_081_4,
            1.8472041376831578,
            0.260_638_187_267_738_6,
            1.6287639169258714,
            0.606_190_554_392_843_5,
            0.631_151_139_455_903_8,
            1.1157241030662073,
            0.977_229_819_457_283_3,
            0.887_046_232_545_975_9,
            0.864_980_324_904_755_7,
            0.909_854_555_540_833_3,
            0.886_196_803_281_563_9,
            1.6346442254821438,
            0.596_092_405_897_667_8,
            1.090_408_615_513_094,
            1.0041175765198984,
            1.109300888123117,
            1.541_110_163_763_897,
        ];
        let df1b = [
            6.0, 8.0, 10.0, 6.0, 8.0, 10.0, 6.0, 8.0, 10.0, 6.0, 8.0, 10.0, 6.0, 8.0, 10.0, 6.0,
            8.0, 10.0, 6.0, 8.0, 10.0, 6.0, 8.0, 10.0, 6.0, 8.0, 10.0, 6.0, 8.0, 10.0,
        ];
        let scale2 = [0.974_552_543_088_280_2_f64; 30];
        let dp2 = [8_134.844_780_382_662_f64; 30];
        check_robust_unequal("R2", &x2, &df1b, None, &scale2, &dp2);

        // R3: covariate trend with right outliers -> per-gene scale.
        let x3 = [
            0.539_044_219_086_069_3,
            0.23661776187622388,
            0.076_201_682_711_961_72,
            0.630_086_743_998_138,
            0.650_987_435_045_026_3,
            0.842_278_820_319_909_2,
            0.491_844_203_672_007_9,
            6.448_427_301_761_955,
            1.2087370452481183,
            0.239_248_212_000_354,
            1.2568382864740935,
            0.37311254857129783,
            0.27671393427180063,
            0.574_850_427_867_031_8,
            0.273_063_551_556_959_1,
            1.080730106245529,
            1.2795702810077614,
            0.611_665_664_175_691_7,
            1.7980634165113292,
            0.663_487_937_134_490_8,
            0.996_807_045_846_497_5,
            18.128946787721024,
            1.7142965169975035,
            1.449_571_042_312_149,
            1.4930161756216276,
            2.187_545_894_935_803,
            1.8192882143145384,
            2.09582235145038,
            1.4058722658383287,
            0.612_227_999_262_316_5,
            14.493_470_715_651_85,
            1.587_104_034_660_905,
            2.1384995210438946,
            0.942_027_524_800_738_4,
            1.7303158602514568,
            0.867_230_571_621_344,
            2.3372774698286336,
            2.0163275306516515,
            1.0986912334131065,
            4.308111430321893,
        ];
        let df1c = [
            5.0, 7.0, 9.0, 5.0, 7.0, 9.0, 5.0, 7.0, 9.0, 5.0, 7.0, 9.0, 5.0, 7.0, 9.0, 5.0, 7.0,
            9.0, 5.0, 7.0, 9.0, 5.0, 7.0, 9.0, 5.0, 7.0, 9.0, 5.0, 7.0, 9.0, 5.0, 7.0, 9.0, 5.0,
            7.0, 9.0, 5.0, 7.0, 9.0, 5.0,
        ];
        let cov3 = [
            -1.6973582475323723,
            -1.653386051863265,
            -1.5197916315627742,
            -1.3994968980185873,
            -1.3770853614073897,
            -0.8685822286424798,
            -0.851_879_364_202_679_1,
            -0.755_372_870_094_767_1,
            -0.693_124_605_092_816_8,
            -0.615_876_899_416_812_8,
            -0.587_287_589_545_757_5,
            -0.525_642_223_659_542_9,
            -0.512_363_475_286_970_1,
            -0.509_108_805_762_841_8,
            -0.47358310536076076,
            -0.17304879792837147,
            -0.163_159_551_650_866_8,
            -0.099_030_693_969_731_35,
            -0.056582774046561445,
            0.008_098_727_936_676_124,
            0.044_483_626_979_542_84,
            0.13289179603471196,
            0.198_091_361_808_575_4,
            0.2614007552112268,
            0.28417729172132883,
            0.35626653938190983,
            0.49054607057569327,
            0.600_740_038_173_888_9,
            0.666_223_070_296_786_5,
            0.748_046_803_709_493_8,
            0.784_220_371_322_524_8,
            0.9592412553828471,
            0.965_502_989_259_626_5,
            0.9808579035910725,
            1.1154412192193877,
            1.2203256484786897,
            1.274_794_669_462_531,
            2.290_936_942_882_79,
            2.5215299924233707,
            2.7975228835504256,
        ];
        let scale3 = [
            0.30716070090370096,
            0.317_044_429_951_319_3,
            0.348_906_926_874_468_6,
            0.380_090_527_948_858_6,
            0.38617781579613586,
            0.551_402_312_268_497_8,
            0.557_818_327_467_101_5,
            0.5962839627315033,
            0.622_396_765_445_769_2,
            0.656_272_944_879_835_6,
            0.6692315210045654,
            0.697_959_466_422_480_9,
            0.704_289_195_891_229_6,
            0.705_848_327_680_508_3,
            0.723_064_434_898_611_9,
            0.883579621760054,
            0.889_319_841_175_482_8,
            0.927_241_684_104_796_7,
            0.952_993_704_816_834,
            0.993_185_060_161_913_9,
            1.0162674207322804,
            1.0736021018381559,
            1.116_799_273_796_704,
            1.1592113970340627,
            1.1745257010455081,
            1.2229745666800398,
            1.312_517_942_704_916,
            1.3745904760509138,
            1.3991833082340708,
            1.426_894_061_899_994,
            1.4389445458368257,
            1.501_236_824_376_662,
            1.5036332566235469,
            1.5095610319005872,
            1.564_432_362_592_114,
            1.610322273457472,
            1.635031881127248,
            2.154_678_058_942_859,
            2.2803239986430137,
            2.433_176_317_272_599,
        ];
        let mut dp3 = [12.208761964849476_f64; 40];
        for k in [7usize, 21, 30] {
            dp3[k] = 2.7408973618244814;
        }
        check_robust_unequal("R3", &x3, &df1c, Some(&cov3), &scale3, &dp3);

        // R4: NA variances + right outliers, no covariate.
        let nan = f64::NAN;
        let x4 = [
            1.023_213_413_751_331,
            0.743480250438045,
            0.444_993_410_271_855_6,
            nan,
            0.050383113595066184,
            2.216_649_724_265_139,
            1.049_203_365_238_416,
            0.695_472_383_043_613_1,
            9.332_089_834_926_714,
            1.200092163119983,
            0.903_454_629_393_505_1,
            2.0236715545973434,
            1.1135951813877558,
            0.977_174_440_448_530_9,
            0.622_307_200_402_812_7,
            0.447_863_402_142_726,
            nan,
            0.49394659653033035,
            2.1605440089890657,
            1.6931692959938356,
            0.564_721_703_951_592,
            0.34583254193314156,
            1.0141939154261466,
            28.808337554746064,
            1.6537857037538832,
            0.267_413_061_360_354_1,
            1.13717505245136,
            0.779_492_650_396_858_7,
            0.619_681_213_679_133_8,
            nan,
            0.900_056_680_281_428_5,
            0.521_829_030_859_030_7,
            1.759024302473003,
            0.21310939936457118,
            0.628_046_189_146_188_1,
            1.1759956568963894,
        ];
        let df1d = [
            4.0, 6.0, 8.0, 10.0, 4.0, 6.0, 8.0, 10.0, 4.0, 6.0, 8.0, 10.0, 4.0, 6.0, 8.0, 10.0,
            4.0, 6.0, 8.0, 10.0, 4.0, 6.0, 8.0, 10.0, 4.0, 6.0, 8.0, 10.0, 4.0, 6.0, 8.0, 10.0,
            4.0, 6.0, 8.0, 10.0,
        ];
        let scale4 = [0.923_414_484_086_539_8_f64; 36];
        let mut dp4 = [17.964906063961386_f64; 36];
        dp4[8] = 7.908450177638727;
        dp4[23] = 2.3144253170237468;
        check_robust_unequal("R4", &x4, &df1d, None, &scale4, &dp4);
    }

    /// Compare `fit_fdist_trend` against `limma:::fitFDist(x, df1, covariate)`.
    fn check_trend(
        tag: &str,
        x: &[f64],
        df1: &[f64],
        cov: &[f64],
        want_scale: &[f64],
        want_df2: f64,
    ) {
        let (scale, df2) = fit_fdist_trend(
            &Array1::from(x.to_vec()),
            &Array1::from(df1.to_vec()),
            &Array1::from(cov.to_vec()),
        )
        .unwrap();
        assert_eq!(scale.len(), want_scale.len(), "{tag}: scale length");
        for (i, (&g, &w)) in scale.iter().zip(want_scale.iter()).enumerate() {
            assert!(rel(g, w) < 1e-7, "{tag}: scale[{i}] = {g} vs {w}");
        }
        if want_df2.is_infinite() {
            assert!(
                df2.is_infinite() && df2 > 0.0,
                "{tag}: df2 = {df2}, want +Inf"
            );
        } else {
            assert!(
                rel(df2, want_df2) < 1e-7,
                "{tag}: df2 = {df2} vs {want_df2}"
            );
        }
    }

    // Reference: limma 3.68.3 fitFDist(x, df1, covariate) on the `notallok`
    // (spline-extrapolation) branch — genes with NA/negative variance or zero
    // df1 are dropped from the spline fit, then the trend is predicted back to
    // them (R's `predict.ns`). Covers interior prediction, extrapolation beyond
    // both boundary knots, and ±Inf covariate clamping.
    #[test]
    fn fit_fdist_trend_notallok_matches_r() {
        let na = f64::NAN;

        // T1: NA variances dropped, scalar df1, nok=27 -> evar<=0 (df2=Inf).
        let x1 = [
            0.668_050_109_540_809_2,
            1.300_486_703_906_259_8,
            0.591_123_993_662_068_2,
            1.412_182_169_606_417_6,
            na,
            1.709_392_079_592_422,
            0.561_697_540_152_080_8,
            0.331_753_569_556_590_74,
            0.702_404_178_640_605_4,
            0.177_388_217_223_870_2,
            1.551_393_054_288_455,
            0.831_676_142_545_982_7,
            0.773_268_079_682_792_1,
            na,
            0.911_149_742_257_311_5,
            0.495_051_838_979_765_1,
            1.129_617_165_840_778_5,
            0.098_872_546_868_354_41,
            0.672_670_029_898_832,
            0.735_263_172_100_904_2,
            0.486_887_217_737_688_56,
            na,
            1.234_402_923_334_334_3,
            4.127_830_275_907_668,
            1.388_065_239_430_386,
            0.394_749_041_568_978_05,
            0.216_216_083_567_096_82,
            0.801_846_236_298_862_3,
            0.554_442_030_140_825_8,
            1.569_493_488_930_550_5,
        ];
        let cov1 = [
            -1.818_934_500_797_971,
            -1.728_927_283_634_95,
            -1.411_173_043_656_684_3,
            -1.277_946_167_246_393_2,
            -1.051_186_696_922_955,
            -1.037_444_582_736_065_9,
            -0.687_798_326_082_019_1,
            -0.641_357_554_081_644_3,
            -0.543_881_103_817_260,
            -0.511_375_321_765_840_1,
            -0.427_863_231_850_420_5,
            -0.259_802_107_839_475_27,
            -0.077_603_595_448_649_43,
            -0.074_823_364_922_382_83,
            -0.057_649_748_482_495_21,
            -0.050_984_123_761_055_82,
            0.112_457_509_152_447_9,
            0.138_339_047_584_972_37,
            0.148_902_488_663_004_62,
            0.165_380_870_625_706_3,
            0.302_492_245_643_733_7,
            0.386_835_291_215_675,
            0.422_604_331_138_563_4,
            1.111_675_269_537_822_5,
            1.129_809_120_495_224_5,
            1.153_158_166_984_797,
            1.173_722_464_235_718_2,
            1.509_897_438_118_779,
            1.619_937_007_939_922_5,
            1.852_147_574_930_227_6,
        ];
        let s1 = [
            1.348_567_234_047_103_8,
            1.296_347_719_097_414,
            1.130_914_668_754_533_7,
            1.070_993_415_947_367_5,
            0.982_372_808_960_233_8,
            0.977_547_292_047_697,
            0.875_717_577_057_818_7,
            0.865_220_263_201_161_5,
            0.845_500_198_322_507_9,
            0.839_623_038_183_936_6,
            0.826_136_748_026_259_4,
            0.806_136_936_846_043_6,
            0.795_621_441_356_589,
            0.795_554_680_444_406_7,
            0.795_206_588_456_107_5,
            0.795_101_427_402_649_9,
            0.797_644_530_437_386_9,
            0.798_921_842_451_763,
            0.799_509_937_488_619_4,
            0.800_504_103_827_121_2,
            0.812_340_081_012_985_6,
            0.822_714_674_823_542_3,
            0.827_813_167_232_808_7,
            1.005_772_610_665_043_5,
            1.012_504_564_481_419,
            1.021_328_158_076_887_8,
            1.029_244_429_861_418_7,
            1.177_795_277_758_329_8,
            1.234_096_945_183_988_7,
            1.364_490_970_760_690_7,
        ];
        check_trend("T1", &x1, &[4.0], &cov1, &s1, f64::INFINITY);

        // T2: per-gene df1 with four zeros dropped (nok=21), finite df2.
        let x2 = [
            0.393_968_498_746_662,
            0.511_726_219_974_134_5,
            0.212_270_353_770_316_2,
            0.534_864_711_251_273_9,
            0.232_128_306_222_922_24,
            0.047_712_491_375_574_21,
            1.248_271_133_535_308_7,
            0.263_067_344_344_326_2,
            0.790_753_407_450_641_7,
            0.469_837_994_117_615_34,
            0.315_809_258_921_149_5,
            0.348_058_507_399_016_5,
            1.206_002_476_299_393_3,
            0.781_211_309_940_489_5,
            1.665_074_463_600_860_5,
            0.347_507_409_684_856_27,
            0.621_200_042_189_179_8,
            0.729_814_667_722_752_8,
            0.608_097_395_142_041_8,
            1.309_015_590_730_340_2,
            0.804_008_546_611_541_3,
            1.194_150_451_327_983,
            1.040_493_996_890_242_7,
            0.394_479_284_163_874_8,
            0.880_330_595_306_281,
        ];
        let df2v = [
            5.0, 5.0, 0.0, 5.0, 5.0, 5.0, 5.0, 0.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 0.0,
            5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 0.0, 5.0,
        ];
        let cov2 = [
            -2.0,
            -1.791_666_666_666_666_7,
            -1.583_333_333_333_333_3,
            -1.375,
            -1.166_666_666_666_666_5,
            -0.958_333_333_333_333_3,
            -0.75,
            -0.541_666_666_666_666_5,
            -0.333_333_333_333_333_26,
            -0.125,
            0.083_333_333_333_333_48,
            0.291_666_666_666_666_96,
            0.5,
            0.708_333_333_333_333_5,
            0.916_666_666_666_667,
            1.125,
            1.333_333_333_333_333_5,
            1.541_666_666_666_667,
            1.75,
            1.958_333_333_333_333_5,
            2.166_666_666_666_667,
            2.375,
            2.583_333_333_333_334,
            2.791_666_666_666_667,
            3.0,
        ];
        let s2 = [
            0.375_866_847_655_144_7,
            0.397_045_606_572_135_6,
            0.419_412_836_004_638_63,
            0.443_029_800_802_343_94,
            0.467_960_301_921_244_6,
            0.494_270_719_668_742_7,
            0.522_030_051_141_576_6,
            0.551_309_941_054_917_8,
            0.582_184_705_109_587,
            0.614_731_344_992_770_4,
            0.649_029_554_056_260_4,
            0.685_161_712_665_588_2,
            0.723_212_872_163_913_8,
            0.763_273_684_911_177_8,
            0.805_450_545_345_930_1,
            0.849_859_176_853_651_3,
            0.896_622_442_547_367_9,
            0.945_870_847_585_434_6,
            0.997_743_080_190_036_7,
            1.052_386_594_585_068,
            1.109_958_239_362_042_4,
            1.170_624_935_100_712_7,
            1.234_564_405_420_003,
            1.301_965_966_017_903_3,
            1.373_031_376_679_661,
        ];
        check_trend("T2", &x2, &df2v, &cov2, &s2, 60.374_728_318_686_31);

        // T3: the four dropped genes sit at the covariate extremes, so the trend
        // is *extrapolated* linearly beyond both boundary knots.
        let x3 = [
            na,
            na,
            0.691_711_843_436_850_6,
            0.622_630_017_386_997_3,
            0.494_243_503_394_316,
            0.938_759_308_781_728_3,
            0.460_653_170_893_032_2,
            0.880_526_577_745_147_4,
            0.755_940_839_358_170_9,
            0.626_842_012_588_787_4,
            0.952_557_139_139_847_7,
            1.279_121_040_856_989_7,
            0.794_193_190_170_537_2,
            1.255_700_888_024_271_3,
            0.655_995_846_018_213_5,
            1.755_978_006_353_944_5,
            1.652_006_575_519_735,
            1.079_323_278_760_472_3,
            1.693_050_241_693_432,
            1.402_149_156_974_183_5,
            0.958_589_262_459_760,
            1.829_534_294_870_645,
            1.550_436_258_458_263_2,
            1.575_031_982_009_485_8,
            2.535_692_986_895_965_5,
            1.416_895_207_861_961_7,
            na,
            na,
        ];
        let cov3 = [
            -1.891_182_416_943_259_6,
            -1.666_109_282_557_862_5,
            -1.476_381_394_207_286_8,
            -1.476_272_416_371_422_7,
            -1.091_494_505_528_445_7,
            -1.090_290_133_251_563_8,
            -1.017_177_492_279_665_3,
            -0.867_315_092_720_626_8,
            -0.745_408_652_670_449_8,
            -0.741_015_812_857_175_1,
            -0.677_748_882_644_702_1,
            -0.347_012_052_112_275_37,
            -0.303_889_624_932_322_97,
            0.102_873_317_075_009_8,
            0.221_423_360_187_956_86,
            0.418_257_493_593_159_67,
            0.537_119_740_704_388_8,
            0.576_701_769_185_424_5,
            0.833_412_786_547_817_4,
            0.845_660_571_442_754_3,
            0.911_991_506_245_756_2,
            1.161_349_619_383_321,
            1.168_829_982_095_781,
            1.273_770_716_555_102_2,
            1.415_588_600_613_468_6,
            1.417_529_081_968_673_2,
            1.532_444_749_036_084_2,
            1.970_839_814_188_167_3,
        ];
        let s3 = [
            0.637_131_884_149_240_9,
            0.690_785_590_177_825_7,
            0.739_508_117_708_049_5,
            0.739_537_068_543_450_7,
            0.849_310_923_288_389_5,
            0.849_679_873_790_225,
            0.872_398_347_921_150_6,
            0.921_021_203_919_749_7,
            0.962_744_435_484_967_7,
            0.964_286_148_350_636_5,
            0.986_795_924_296_255_7,
            1.114_512_738_865_558,
            1.132_516_050_809_078_4,
            1.320_232_229_914_269_3,
            1.381_772_413_551_737_7,
            1.491_667_167_764_201,
            1.562_951_502_145_707_2,
            1.587_548_008_022_439_8,
            1.757_993_189_042_720_7,
            1.766_614_969_623_39,
            1.814_112_145_534_830_1,
            2.005_233_486_689_923_5,
            2.011_283_694_592_009,
            2.098_176_111_347_731_5,
            2.221_722_436_895_536,
            2.223_462_779_155_565_4,
            2.328_994_415_482_243,
            2.779_685_778_433_139,
        ];
        check_trend("T3", &x3, &[5.0], &cov3, &s3, f64::INFINITY);

        // T4: all genes usable, but covariate has -Inf and +Inf -> clamp to the
        // finite range ±1 before fitting.
        let x4 = [
            1.460_684_658_159_952,
            0.929_320_089_464_292_1,
            1.004_561_388_965_886_2,
            1.359_002_066_955_692_6,
            0.937_218_928_892_790_1,
            1.937_957_094_827_315_7,
            1.100_890_963_710_844_6,
            1.285_515_962_577_706_4,
            1.400_028_989_125_639_9,
            1.948_970_758_855_427_7,
            0.597_534_506_383_883_4,
            0.588_098_119_136_306_6,
            1.247_633_370_475_978_6,
            1.013_165_829_348_660_7,
            0.873_404_557_452_135,
            0.466_607_779_576_223_46,
            1.331_846_007_727_750_3,
            1.226_078_586_611_032_6,
            1.456_974_504_888_057_8,
            0.935_437_091_071_497_2,
            1.878_429_661_609_096,
            0.930_717_449_286_479_7,
            1.035_465_770_520_301_6,
            0.957_973_578_912_317_4,
            1.854_466_168_248_283_3,
            0.763_294_910_176_882_9,
        ];
        let cov4 = [
            0.773_639_322_679_174_7,
            0.059_688_352_991_138_07,
            -0.317_072_866_123_578_55,
            -0.509_246_573_590_155_6,
            -1.153_036_622_104_934,
            1.052_241_278_586_400_4,
            f64::NEG_INFINITY,
            0.575_769_841_586_934,
            -0.242_343_833_265_866_95,
            1.636_407_646_700_849_5,
            -0.208_163_180_594_472_72,
            -0.879_915_899_267_383_5,
            -0.572_403_354_610_630_7,
            -0.012_328_118_511_838_04,
            -1.433_772_725_748_814,
            -0.984_699_712_984_733_5,
            0.163_666_275_174_588_95,
            1.260_036_251_887_696_4,
            f64::INFINITY,
            -0.722_422_472_389_862_6,
            1.608_378_354_296_734_7,
            -0.347_530_581_042_343_86,
            -1.074_226_364_095_721_7,
            0.697_988_268_711_317_5,
            0.897_482_808_214_027_4,
            -1.141_727_632_668_820_6,
        ];
        let s4 = [
            1.659_304_790_599_606,
            1.410_300_900_749_76,
            1.314_788_110_712_838_3,
            1.274_489_288_337_632_4,
            1.170_487_066_419_548_6,
            1.783_268_130_752_830_6,
            1.041_775_265_812_246_8,
            1.580_764_750_466_000_3,
            1.331_907_893_243_926_3,
            2.098_221_380_490_413,
            1.340_030_458_519_045_2,
            1.209_624_934_861_425_4,
            1.262_320_411_922_209_6,
            1.390_235_928_160_515_4,
            1.136_099_857_058_694,
            1.193_853_528_831_997_1,
            1.440_830_421_860_32,
            1.886_460_481_436_271,
            2.831_039_106_202_752_5,
            1.235_330_911_132_471_3,
            2.081_305_738_731_465,
            1.308_053_044_955_874_7,
            1.181_141_042_233_640_8,
            1.628_370_489_755_436,
            1.712_428_617_774_108_4,
            1.171_986_411_871_814_2,
        ];
        check_trend("T4", &x4, &[4.0], &cov4, &s4, f64::INFINITY);

        // T5: a negative variance (idx3) and a zero df1 (idx10) dropped — both
        // interior to the covariate range, so they are interpolated.
        let x5 = [
            0.429_178_288_219_323_63,
            0.361_629_848_232_283_06,
            0.491_374_690_982_892_3,
            -0.001,
            0.501_012_501_661_693_6,
            0.879_948_740_279_662_7,
            1.497_207_667_953_139_6,
            1.250_541_607_815_663_5,
            0.952_132_854_010_289,
            0.871_035_936_902_676_4,
            0.783_574_923_522_921_8,
            1.688_002_308_919_733_1,
            0.677_018_707_358_007_7,
            1.628_738_342_835_098,
            1.302_988_874_606_991_5,
            1.848_202_373_630_821,
            0.356_945_916_106_020_1,
            1.171_857_737_395_632_5,
            0.699_580_841_179_741_8,
            1.081_768_203_975_455_3,
        ];
        let df5v = [
            8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 0.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0,
            8.0, 8.0, 8.0,
        ];
        let cov5 = [
            0.0,
            0.263_157_894_736_842_1,
            0.526_315_789_473_684_2,
            0.789_473_684_210_526_3,
            1.052_631_578_947_368_4,
            1.315_789_473_684_210_4,
            1.578_947_368_421_052_7,
            1.842_105_263_157_894_7,
            2.105_263_157_894_736_7,
            2.368_421_052_631_578_8,
            2.631_578_947_368_421,
            2.894_736_842_105_263,
            3.157_894_736_842_105_3,
            3.421_052_631_578_947_3,
            3.684_210_526_315_789_4,
            3.947_368_421_052_631_4,
            4.210_526_315_789_473_5,
            4.473_684_210_526_316,
            4.736_842_105_263_157_5,
            5.0,
        ];
        let s5 = [
            0.421_945_688_265_181_64,
            0.496_991_060_626_983_4,
            0.583_814_156_179_060_1,
            0.682_132_444_698_899_4,
            0.790_614_642_193_487_8,
            0.906_561_103_796_838_2,
            1.025_650_538_328_756_9,
            1.141_841_611_559_546,
            1.247_528_460_987_988,
            1.334_034_995_926_419_3,
            1.392_483_850_295_140_7,
            1.416_327_109_876_026_8,
            1.406_606_780_107_036_5,
            1.368_085_106_724_575_8,
            1.307_014_262_587_435_2,
            1.230_183_433_811_377_7,
            1.144_135_050_993_384,
            1.054_625_112_265_449_5,
            0.966_335_335_515_517_2,
            0.882_799_502_530_530_9,
        ];
        check_trend("T5", &x5, &df5v, &cov5, &s5, f64::INFINITY);
    }

    // Reference: limma 3.68.3 squeezeVar(var, df) with unequal df, which routes
    // through fitFDistUnequalDF1 (non-legacy, non-robust) rather than fitFDist.
    #[test]
    fn squeeze_var_unequal_df_matches_r() {
        let var = Array1::from(vec![
            0.8, 1.2, 0.5, 2.1, 0.3, 1.7, 0.9, 3.4, 0.6, 1.1, 1.45, 0.72,
        ]);
        let df = Array1::from(vec![
            3.0, 5.0, 3.0, 8.0, 4.0, 6.0, 3.0, 10.0, 5.0, 4.0, 7.0, 6.0,
        ]);
        let sq = squeeze_var(&var, &df, None, false).unwrap();
        assert!(
            rel(sq.var_prior[0], 1.484_165_903_978_432) < 1e-9,
            "scale={}",
            sq.var_prior[0]
        );
        assert!(
            rel(sq.df_prior[0], 27.213_028_142_992_464) < 1e-6,
            "df={}",
            sq.df_prior[0]
        );
        let want_post = [
            1.416_231_710_086_273_4,
            1.440_058_609_451_970_7,
            1.386_443_236_195_458_2,
            1.624_076_415_172_368_3,
            1.332_413_129_649_257,
            1.523_156_765_352_288_4,
            1.426_161_201_383_212,
            1.998_994_766_778_821,
            1.346_928_588_061_770_7,
            1.434_934_422_531_838_8,
            1.477_175_545_602_393_2,
            1.346_117_804_174_621_6,
        ];
        for (i, &w) in want_post.iter().enumerate() {
            assert!(
                rel(sq.var_post[i], w) < 1e-6,
                "post[{i}]={}",
                sq.var_post[i]
            );
        }
    }

    // Reference: limma 3.68.3 fitFDistRobustly(x, df1=13) on the fixture vector.
    #[test]
    fn fit_fdist_robustly_matches_r_on_fixture() {
        let x = fixture_x(50);
        let fit = fit_fdist_robustly(&x, 13.0);
        assert!(
            rel(fit.scale[0], 0.691_643_229_362_652_4) < 1e-10,
            "scale={}",
            fit.scale[0]
        );

        let sum: f64 = fit.df2_shrunk.iter().sum();
        let min = fit.df2_shrunk.iter().cloned().fold(f64::INFINITY, f64::min);
        // Non-outlier genes sit at the global robust df2 (limma f$df2 = 3.3726).
        let max = fit
            .df2_shrunk
            .iter()
            .cloned()
            .fold(f64::NEG_INFINITY, f64::max);
        assert!(rel(sum, 152.336_194_363_376) < 1e-10, "sum={sum}");
        assert!(rel(min, 2.966_970_204_216_085_7) < 1e-10, "min={min}");
        assert!(rel(max, 3.372_566_440_190_4) < 1e-10, "max={max}");

        // Injected outliers shrink to the minimum df2; gene 1 stays at df2.
        assert!(rel(fit.df2_shrunk[9], 2.966_970_204_216_085_7) < 1e-10);
        assert!(rel(fit.df2_shrunk[19], 2.966_970_204_216_085_7) < 1e-10);
        assert!(rel(fit.df2_shrunk[0], 3.372_566_440_190_4) < 1e-10);
        assert!(rel(fit.df2_shrunk[49], 2.966_970_204_216_085_7) < 1e-10);
    }

    // Exercises the uniroot branch: a smooth base with strong outliers makes the
    // non-robust df2 (2.276) too small, so robust estimation pushes it back up.
    // Reference: limma 3.68.3 fitFDistRobustly(xf, df1=13).
    #[test]
    fn fit_fdist_robustly_uniroot_branch() {
        let n = 60;
        let mut x = vec![0.0_f64; n];
        for i in 1..=n {
            x[i - 1] = (((i % 11) as f64 - 5.0) / 9.0).exp();
        }
        for &k in &[6usize, 16, 26, 36, 46, 56] {
            x[k - 1] *= 60.0;
        }
        let fit = fit_fdist_robustly(&x, 13.0);
        assert!(
            rel(fit.scale[0], 1.144_677_907_473_284_5) < 1e-9,
            "scale={}",
            fit.scale[0]
        );
        // Non-outlier genes sit at the global robust df2 (limma f$df2 = 26.5086).
        let max = fit
            .df2_shrunk
            .iter()
            .cloned()
            .fold(f64::NEG_INFINITY, f64::max);
        assert!(rel(max, 26.508_557_634_070_502) < 1e-9, "max={max}");
        let sum: f64 = fit.df2_shrunk.iter().sum();
        assert!(rel(sum, 1_433.095_052_835_895) < 1e-9, "sum={sum}");
        // Non-outliers stay at df2; the six outliers shrink to ~0.272.
        assert!(rel(fit.df2_shrunk[0], 26.508_557_634_070_502) < 1e-9);
        assert!(rel(fit.df2_shrunk[29], 26.508_557_634_070_502) < 1e-9);
        assert!(rel(fit.df2_shrunk[5], 0.272_156_765_986_153) < 1e-7);
        assert!(rel(fit.df2_shrunk[15], 0.272_156_765_986_621) < 1e-7);
        assert!(rel(fit.df2_shrunk[55], 0.272_156_766_106_002) < 1e-7);
    }

    // Exercises the df2 = Inf branch with outlier shrinkage: a near-constant base
    // (under-dispersed -> df2 = Inf) with three huge outliers whose posterior df2
    // shrinks to extremely small finite values via the deep chi-square tail.
    // Reference: limma 3.68.3 fitFDistRobustly(xh, df1=13).
    #[test]
    fn fit_fdist_robustly_inf_branch_with_outliers() {
        let n = 60;
        let mut x = vec![0.0_f64; n];
        for i in 1..=n {
            x[i - 1] = 0.45 + 0.10 * ((i % 7) as f64) / 7.0;
        }
        for &k in &[5usize, 25, 45] {
            x[k - 1] *= 100.0;
        }
        let fit = fit_fdist_robustly(&x, 13.0);
        assert!(
            rel(fit.scale[0], 0.535_669_567_990_593_8) < 1e-9,
            "scale={}",
            fit.scale[0]
        );
        // Global robust df2 = Inf, so non-outlier genes keep an infinite df2.
        let n_finite = fit.df2_shrunk.iter().filter(|v| v.is_finite()).count();
        assert_eq!(n_finite, 3, "n_finite={n_finite}");
        assert!(fit.df2_shrunk[0].is_infinite());
        assert!(rel(fit.df2_shrunk[4], 1.367_489_954_146_42e-257) < 1e-7);
        assert!(rel(fit.df2_shrunk[24], 1.321_157_077_384_01e-250) < 1e-7);
        assert!(rel(fit.df2_shrunk[44], 2.287_533_079_578_86e-243) < 1e-7);
    }
}