limma-rust 0.1.0

//! Multiple-testing decisions. Port of limma's `decideTests` (`decidetests.R`):
//! the `separate`, `global`, `hierarchical` and `nestedF` methods, plus the
//! `p.adjust` family limma uses (`none`, `bonferroni`, `holm`, `BH`/`fdr`,
//! `BY`) and the row-wise `.classifyTestsP` classifier ([`classify_tests_p`]).
//! For `separate`/`global`, NA p-values (missing-value data) are treated as
//! NotSig rather than R's `NA` outcome; `nestedF` errors on any NA
//! `F.p.value`, as in limma. Complete-data inputs match R exactly.

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

use crate::classifytestsf::classify_tests_f;
use crate::fit::MArrayLM;
use crate::special::{ln_norm_cdf, t_two_sided_pvalue};

/// p-value adjustment method, matching the subset of R's `p.adjust.methods`
/// accepted by `decideTests`.
#[derive(Clone, Copy, PartialEq, Eq, Debug)]
pub enum Adjust {
    None,
    Bonferroni,
    Holm,
    BH,
    BY,
}

impl Adjust {
    pub fn parse(s: &str) -> Result<Self> {
        Ok(match s.to_ascii_lowercase().as_str() {
            "none" => Adjust::None,
            "bonferroni" => Adjust::Bonferroni,
            "holm" => Adjust::Holm,
            "bh" | "fdr" => Adjust::BH,
            "by" => Adjust::BY,
            other => bail!(
                "invalid adjust.method '{}' (none|bonferroni|holm|BH|fdr|BY)",
                other
            ),
        })
    }
}

/// `decideTests` classification method (the full `c("separate","global",
/// "hierarchical","nestedF")` set). Which methods are supported depends on the
/// entry point: [`decide_tests`] (MArrayLM) handles `separate`/`global`/
/// `nestedF` and defers `hierarchical` to the p-value path; whereas
/// [`decide_tests_pvalues`] (p-value matrix) handles `separate`/`global`/
/// `hierarchical` and rejects `nestedF` (which requires a fit), mirroring limma.
#[derive(Clone, Copy, PartialEq, Eq, Debug)]
pub enum DecideMethod {
    Separate,
    Global,
    Hierarchical,
    NestedF,
}

impl DecideMethod {
    pub fn parse(s: &str) -> Result<Self> {
        Ok(match s.to_ascii_lowercase().as_str() {
            "separate" => DecideMethod::Separate,
            "global" => DecideMethod::Global,
            "hierarchical" => DecideMethod::Hierarchical,
            "nestedf" => DecideMethod::NestedF,
            other => bail!(
                "decideTests method '{}' not supported (separate|global|hierarchical|nestedF)",
                other
            ),
        })
    }
}

/// Outcome matrix from `decideTests`: -1 (down), 0 (not significant), +1 (up).
/// Port of limma's `TestResults` for the `c(-1, 0, 1)` / `c("Down","NotSig","Up")`
/// level set.
#[derive(Clone, Debug)]
pub struct TestResults {
    pub data: Array2<i8>, // n_genes x n_coef
    pub gene_names: Vec<String>,
    pub coef_names: Vec<String>,
}

impl TestResults {
    /// Per-column (down, not-sig, up) counts, matching `summary.TestResults`.
    pub fn counts(&self) -> Vec<(usize, usize, usize)> {
        (0..self.data.ncols())
            .map(|j| {
                let mut down = 0;
                let mut up = 0;
                let mut notsig = 0;
                for g in 0..self.data.nrows() {
                    match self.data[[g, j]] {
                        v if v < 0 => down += 1,
                        v if v > 0 => up += 1,
                        _ => notsig += 1,
                    }
                }
                (down, notsig, up)
            })
            .collect()
    }

    /// `summary.TestResults`: the labelled count table with one row per level
    /// (`Down` = -1, `NotSig` = 0, `Up` = +1) and one column per coefficient,
    /// giving how many genes fall at each level. Mirrors limma's table with the
    /// `c("Down","NotSig","Up")` labels attached by `decideTests`.
    pub fn summary(&self) -> TestResultsSummary {
        let mut down = Vec::with_capacity(self.coef_names.len());
        let mut notsig = Vec::with_capacity(self.coef_names.len());
        let mut up = Vec::with_capacity(self.coef_names.len());
        for (d, n, u) in self.counts() {
            down.push(d);
            notsig.push(n);
            up.push(u);
        }
        TestResultsSummary {
            coef_names: self.coef_names.clone(),
            down,
            notsig,
            up,
        }
    }
}

/// Labelled count table returned by [`TestResults::summary`], the analog of
/// `summary.TestResults`. The three count vectors are aligned with
/// [`Self::LABELS`] / [`Self::LEVELS`] and indexed by coefficient (column).
#[derive(Clone, Debug)]
pub struct TestResultsSummary {
    pub coef_names: Vec<String>,
    pub down: Vec<usize>,
    pub notsig: Vec<usize>,
    pub up: Vec<usize>,
}

impl TestResultsSummary {
    /// Row labels, matching `attr(results, "labels")` set by `decideTests`.
    pub const LABELS: [&'static str; 3] = ["Down", "NotSig", "Up"];
    /// Row levels, matching `attr(results, "levels")` set by `decideTests`.
    pub const LEVELS: [i8; 3] = [-1, 0, 1];
}

/// Adjust a vector of (assumed finite) p-values. Port of `stats::p.adjust` for
/// the methods limma exposes through `decideTests`. The caller filters NAs, so
/// `n == p.len()`, matching R's `o <- !is.na(p)` pre-filtering in `decideTests`.
pub fn p_adjust(p: &[f64], method: Adjust) -> Vec<f64> {
    let n = p.len();
    if n == 0 {
        return Vec::new();
    }
    let nf = n as f64;
    match method {
        Adjust::None => p.to_vec(),
        Adjust::Bonferroni => p.iter().map(|&x| (nf * x).min(1.0)).collect(),
        Adjust::Holm => {
            // Ascending order; multiply the i-th smallest by (n + 1 - i), take a
            // running maximum, clamp to 1. Matches `cummax((n+1-i)*p[o])`.
            let mut o: Vec<usize> = (0..n).collect();
            o.sort_by(|&a, &b| p[a].partial_cmp(&p[b]).unwrap());
            let mut out = vec![0.0_f64; n];
            let mut cummax = f64::NEG_INFINITY;
            for (rank0, &idx) in o.iter().enumerate() {
                let i = (rank0 + 1) as f64;
                let val = (nf + 1.0 - i) * p[idx];
                cummax = cummax.max(val);
                out[idx] = cummax.min(1.0);
            }
            out
        }
        Adjust::BH | Adjust::BY => {
            // Descending order; for the rank-i value use factor n/i, accumulate
            // a running minimum, clamp to 1. BY scales by q = sum(1/k, k=1..n).
            let mut o: Vec<usize> = (0..n).collect();
            o.sort_by(|&a, &b| p[b].partial_cmp(&p[a]).unwrap());
            let q = if matches!(method, Adjust::BY) {
                (1..=n).map(|k| 1.0 / k as f64).sum::<f64>()
            } else {
                1.0
            };
            let mut out = vec![0.0_f64; n];
            let mut cummin = f64::INFINITY;
            for (rank0, &idx) in o.iter().enumerate() {
                let i = (n - rank0) as f64; // i runs n, n-1, ..., 1
                let val = q * (nf / i) * p[idx];
                cummin = cummin.min(val);
                out[idx] = cummin.min(1.0);
            }
            out
        }
    }
}

/// Sign of a coefficient as `-1 / 0 / 1` (R's `sign`, with `0.0 -> 0`).
fn sign_i8(x: f64) -> i8 {
    if x > 0.0 {
        1
    } else if x < 0.0 {
        -1
    } else {
        0
    }
}

/// Correlation matrix of the coefficients for `classifyTestsF`, matching limma's
/// guard: any coefficient variance exactly zero (e.g. an all-zero contrast) is
/// reset to 1 before `cov2cor`, avoiding a divide-by-zero in the whitening step.
fn cor_from_cov(cov: &Array2<f64>) -> Array2<f64> {
    let mut c = cov.clone();
    for i in 0..c.nrows() {
        if c[[i, i]] == 0.0 {
            c[[i, i]] = 1.0;
        }
    }
    crate::linalg::cov2cor(&c)
}

/// Classify each (gene, contrast) as up/down/not-significant. Port of
/// `decideTests.MArrayLM` for `method %in% c("separate","global","nestedF")`.
/// Requires a fit that has been through `eBayes` (`separate`/`global` need
/// `p.value`; `nestedF` needs `t` and `F.p.value`).
pub fn decide_tests(
    fit: &MArrayLM,
    method: DecideMethod,
    adjust: Adjust,
    p_value: f64,
    lfc: f64,
) -> Result<TestResults> {
    let coef = &fit.coefficients;
    let (ng, nc) = (coef.nrows(), coef.ncols());
    let mut data = Array2::<i8>::zeros((ng, nc));

    match method {
        DecideMethod::Separate => {
            let p = fit
                .p_value
                .as_ref()
                .ok_or_else(|| anyhow::anyhow!("need to run eBayes first (p.value missing)"))?;
            for j in 0..nc {
                let mut idx = Vec::new();
                let mut vals = Vec::new();
                for g in 0..ng {
                    let v = p[[g, j]];
                    if v.is_finite() {
                        idx.push(g);
                        vals.push(v);
                    }
                }
                let adj = p_adjust(&vals, adjust);
                for (k, &g) in idx.iter().enumerate() {
                    if adj[k] < p_value {
                        data[[g, j]] = sign_i8(coef[[g, j]]);
                    }
                }
            }
        }
        DecideMethod::Global => {
            let p = fit
                .p_value
                .as_ref()
                .ok_or_else(|| anyhow::anyhow!("need to run eBayes first (p.value missing)"))?;
            let mut idx = Vec::new();
            let mut vals = Vec::new();
            for j in 0..nc {
                for g in 0..ng {
                    let v = p[[g, j]];
                    if v.is_finite() {
                        idx.push((g, j));
                        vals.push(v);
                    }
                }
            }
            let adj = p_adjust(&vals, adjust);
            for (k, &(g, j)) in idx.iter().enumerate() {
                if adj[k] < p_value {
                    data[[g, j]] = sign_i8(coef[[g, j]]);
                }
            }
        }
        DecideMethod::Hierarchical => {
            bail!(
                "decideTests(method=\"hierarchical\") on a fit needs the stepwise \
                 .classifyTestsP classifier, which is not yet ported; use \
                 decide_tests_pvalues for the p-value-matrix hierarchical method"
            );
        }
        DecideMethod::NestedF => {
            let t = fit
                .t
                .as_ref()
                .ok_or_else(|| anyhow::anyhow!("need to run eBayes first (t missing)"))?;
            let fp = fit
                .f_p_value
                .as_ref()
                .ok_or_else(|| anyhow::anyhow!("need to run eBayes first (F.p.value missing)"))?;
            if fp.iter().any(|v| v.is_nan()) {
                bail!("nestedF method can't handle NA p-values");
            }
            // Gene-level selection on the (adjusted) overall F p-value.
            let adj = p_adjust(&fp.to_vec(), adjust);
            let sel: Vec<usize> = (0..ng).filter(|&g| adj[g] < p_value).collect();
            let i = sel.len();
            if i > 0 {
                // Adjustment factor `a` shrinks the per-gene threshold so the
                // family-wise/false-discovery control carries into classifyTestsF.
                let n = ng as f64;
                let a = match adjust {
                    Adjust::None => 1.0,
                    Adjust::Bonferroni => 1.0 / n,
                    Adjust::Holm => 1.0 / (n - i as f64 + 1.0),
                    Adjust::BH => i as f64 / n,
                    Adjust::BY => {
                        let harmonic: f64 = (1..=ng).map(|k| 1.0 / k as f64).sum();
                        i as f64 / n / harmonic
                    }
                };
                let cor = cor_from_cov(&fit.cov_coefficients);
                let mut tsel = Array2::<f64>::zeros((i, nc));
                let mut dfsel = Vec::with_capacity(i);
                for (k, &g) in sel.iter().enumerate() {
                    for j in 0..nc {
                        tsel[[k, j]] = t[[g, j]];
                    }
                    // df = df.prior + df.residual (per gene); Inf when no prior.
                    let dfg = match &fit.df_prior {
                        Some(dp) if dp.len() == 1 => dp[0] + fit.df_residual[g],
                        Some(dp) => dp[g] + fit.df_residual[g],
                        None => f64::INFINITY,
                    };
                    dfsel.push(dfg);
                }
                let rsel = classify_tests_f(&tsel, Some(&cor), &dfsel, p_value * a);
                for (k, &g) in sel.iter().enumerate() {
                    for j in 0..nc {
                        data[[g, j]] = rsel[[k, j]];
                    }
                }
            }
        }
    }

    if lfc > 0.0 {
        for j in 0..nc {
            for g in 0..ng {
                if coef[[g, j]].abs() <= lfc {
                    data[[g, j]] = 0;
                }
            }
        }
    }

    Ok(TestResults {
        data,
        gene_names: fit.gene_names.clone(),
        coef_names: fit.coef_names.clone(),
    })
}

/// Classify a matrix of p-values (genes × contrasts) into up/down/not-sig. Port
/// of `decideTests.default`. Supports `separate`, `global` and the Simes
/// `hierarchical` method; `nestedF` is rejected (it requires a fit). When
/// `coefficients` is supplied its sign sets the Up/Down direction and `lfc`
/// zeroes entries with `|coef| < lfc`; otherwise significant entries are
/// reported as `1` (significant, no direction). Uses the `p <= cutoff` rule of
/// `decideTests.default` (the MArrayLM path uses strict `<`).
pub fn decide_tests_pvalues(
    p: &Array2<f64>,
    method: DecideMethod,
    adjust: Adjust,
    p_value: f64,
    lfc: f64,
    coefficients: Option<&Array2<f64>>,
) -> Result<TestResults> {
    if p.iter().any(|v| v.is_nan()) {
        bail!("NA p-values not supported");
    }
    if p.iter().any(|&v| !(0.0..=1.0).contains(&v)) {
        bail!("object doesn't appear to be a matrix of p-values");
    }
    let (ng, nc) = (p.nrows(), p.ncols());
    if let Some(co) = coefficients {
        if co.nrows() != ng || co.ncols() != nc {
            bail!("dim(object) disagrees with dim(coefficients)");
        }
    }

    let mut padj = p.clone();
    let mut cutoff = p_value;

    match method {
        DecideMethod::Separate => {
            for j in 0..nc {
                let col: Vec<f64> = (0..ng).map(|g| p[[g, j]]).collect();
                let adj = p_adjust(&col, adjust);
                for g in 0..ng {
                    padj[[g, j]] = adj[g];
                }
            }
        }
        DecideMethod::Global => {
            let all: Vec<f64> = p.iter().copied().collect();
            let adj = p_adjust(&all, adjust);
            for (v, &a) in padj.iter_mut().zip(adj.iter()) {
                *v = a;
            }
        }
        DecideMethod::Hierarchical => {
            // Simes' method per row to get gene-wise p-values.
            let mult: Vec<f64> = (1..=nc).map(|k| nc as f64 / k as f64).collect();
            let mut genewise = vec![0.0_f64; ng];
            for g in 0..ng {
                let mut row: Vec<f64> = (0..nc).map(|j| p[[g, j]]).collect();
                row.sort_by(|a, b| a.partial_cmp(b).unwrap());
                genewise[g] = (0..nc)
                    .map(|k| row[k] * mult[k])
                    .fold(f64::INFINITY, f64::min);
            }
            let gadj = p_adjust(&genewise, adjust);
            let de: Vec<bool> = gadj.iter().map(|&v| v <= p_value).collect();
            let nde = de.iter().filter(|&&b| b).count();
            for g in 0..ng {
                if de[g] {
                    let row: Vec<f64> = (0..nc).map(|j| p[[g, j]]).collect();
                    let adj = p_adjust(&row, adjust);
                    for j in 0..nc {
                        padj[[g, j]] = adj[j];
                    }
                } else {
                    for j in 0..nc {
                        padj[[g, j]] = 1.0;
                    }
                }
            }
            let a = match adjust {
                Adjust::None => 1.0,
                Adjust::Bonferroni => 1.0 / ng as f64,
                Adjust::Holm => 1.0 / (ng - nde + 1) as f64,
                Adjust::BH => nde as f64 / ng as f64,
                Adjust::BY => {
                    let harmonic: f64 = (1..=ng).map(|k| 1.0 / k as f64).sum();
                    nde as f64 / ng as f64 / harmonic
                }
            };
            cutoff = a * p_value;
        }
        DecideMethod::NestedF => {
            bail!("nestedF adjust method requires an MArrayLM object");
        }
    }

    let mut data = Array2::<i8>::zeros((ng, nc));
    for g in 0..ng {
        for j in 0..nc {
            if padj[[g, j]] <= cutoff {
                let mut v: i8 = 1;
                if let Some(co) = coefficients {
                    if co[[g, j]] < 0.0 {
                        v = -1;
                    }
                    if lfc > 0.0 && co[[g, j]].abs() < lfc {
                        v = 0;
                    }
                }
                data[[g, j]] = v;
            }
        }
    }

    Ok(TestResults {
        data,
        gene_names: vec![String::new(); ng],
        coef_names: vec![String::new(); nc],
    })
}

/// `.classifyTestsP`: classify each gene (row) of a t-statistic matrix as
/// up (`+1`), down (`-1`) or not significant (`0`) by adjusting the two-sided
/// p-values *across the contrasts within each gene* and thresholding at
/// `p_value` (strict `<`).
///
/// `df` is the reference degrees of freedom: pass a single value to share it
/// across genes (limma's default is `Inf`, i.e. a normal reference) or one
/// value per gene (e.g. `df.residual + df.prior` from a fit). The `method`
/// adjustment is applied row-wise via [`p_adjust`]. This is the row-wise
/// companion to [`classify_tests_f`].
pub fn classify_tests_p(
    tstat: &Array2<f64>,
    df: &[f64],
    p_value: f64,
    method: Adjust,
) -> Array2<i8> {
    let (ng, nc) = tstat.dim();
    let mut result = Array2::<i8>::zeros((ng, nc));
    for i in 0..ng {
        let df_i = if df.len() == 1 { df[0] } else { df[i] };
        let prow: Vec<f64> = (0..nc).map(|j| two_sided_p(tstat[[i, j]], df_i)).collect();
        let adj = p_adjust(&prow, method);
        for j in 0..nc {
            if adj[j] < p_value {
                let t = tstat[[i, j]];
                result[[i, j]] = if t > 0.0 {
                    1
                } else if t < 0.0 {
                    -1
                } else {
                    0
                };
            }
        }
    }
    result
}

/// Two-sided Student-t p-value `2*pt(-|t|, df)`, using the normal limit when
/// `df` is infinite.
fn two_sided_p(t: f64, df: f64) -> f64 {
    if t.is_nan() {
        return f64::NAN;
    }
    if df.is_infinite() {
        2.0 * ln_norm_cdf(-t.abs()).exp()
    } else {
        t_two_sided_pvalue(t, df)
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::ebayes::ebayes;
    use crate::fit::lmfit;
    use crate::toptable::p_adjust_bh;
    use ndarray::array;

    /// lmFit + eBayes on `y` (genes × samples) against `design`, defaults.
    fn fit_of(y: &Array2<f64>, design: &Array2<f64>) -> MArrayLM {
        let ng = y.nrows();
        let mut fit = lmfit(
            y,
            design,
            vec![String::new(); ng],
            vec![String::new(); design.ncols()],
        )
        .unwrap();
        ebayes(&mut fit, 0.01, (0.1, 4.0), false, false).unwrap();
        fit
    }

    fn nestedf(fit: &MArrayLM, adjust: Adjust) -> Array2<i8> {
        decide_tests(fit, DecideMethod::NestedF, adjust, 0.1, 0.0)
            .unwrap()
            .data
    }

    #[test]
    fn summary_testresults_matches_r() {
        // Columns C1 = [-1,-1,0,0,1,1], C2 = [1,1,1,0,-1,0]; see
        // scratch/summary_testresults_ref.R.
        let data = array![[-1i8, 1], [-1, 1], [0, 1], [0, 0], [1, -1], [1, 0]];
        let res = TestResults {
            data,
            gene_names: vec![String::new(); 6],
            coef_names: vec!["C1".into(), "C2".into()],
        };
        let s = res.summary();
        assert_eq!(s.down, vec![2, 1]);
        assert_eq!(s.notsig, vec![2, 2]);
        assert_eq!(s.up, vec![2, 3]);
        assert_eq!(s.coef_names, vec!["C1".to_string(), "C2".to_string()]);
        assert_eq!(TestResultsSummary::LABELS, ["Down", "NotSig", "Up"]);
        assert_eq!(TestResultsSummary::LEVELS, [-1, 0, 1]);
    }

    // Case 1: orthogonal group-means design -> cov2cor is the identity and
    // df.prior is Inf, so classifyTestsF takes the chi-squared limit.
    #[test]
    fn nestedf_orthogonal_df_inf_matches_r() {
        let y = array![
            [
                -0.326036490515386,
                0.526448098873642,
                -0.163755665941423,
                0.894937200540396,
                0.482458809858766,
                -1.15025528126137,
                -0.259802107839475,
                1.50989743811878,
                1.85214757493023
            ],
            [
                0.552461855419139,
                -0.794844435415055,
                0.708522104156376,
                3.27915199641374,
                3.758213782399,
                2.72552883818354,
                -3.41117304365668,
                -0.380062992060078,
                -0.888324730462178
            ],
            [
                -0.67494384395583,
                1.42775554468304,
                -0.267980546181617,
                1.00786575031424,
                -2.31932736896749,
                0.577901002861171,
                -0.641357554081644,
                1.1531581669848,
                -0.51137532176584
            ],
            [
                0.214359459043425,
                -1.46681969416347,
                -1.46392175987102,
                -2.07310649057379,
                -0.459504780196151,
                -1.39690264657703,
                0.112457509152448,
                -0.0776035954486494,
                -0.54388110381726
            ],
            [
                2.3107692173136,
                1.7633166213971,
                2.74443582287532,
                1.18985338378285,
                -1.10538366697927,
                0.74905771616478,
                -1.57739566886144,
                -3.81893450079797,
                -3.72892728363495
            ],
            [
                1.1739662875627,
                -0.1933379649975,
                -1.4103901810052,
                -0.724374218369974,
                0.402928273406297,
                -1.05118669692296,
                0.386835291215675,
                -1.03744458273607,
                0.470749538951291
            ],
            [
                0.618789855625968,
                -0.849754740333862,
                0.467067606015246,
                0.167983772348879,
                0.568934917844047,
                0.165380870625706,
                -0.687798326082019,
                0.302492245643734,
                0.00538712217177903
            ],
            [
                -3.11273431475421,
                -2.9415345021505,
                -3.11932010769425,
                3.92033515858084,
                2.29391671757039,
                4.12980912049522,
                0.148902488663005,
                -1.27794616724639,
                1.34804578555275
            ],
            [
                0.917028289512712,
                -0.817670355875958,
                0.467238961765515,
                -1.67160481202292,
                -0.290090625815567,
                1.17372246423572,
                -0.0576497484824952,
                0.138339047584972,
                0.724096713299373
            ],
            [
                -0.223259364627263,
                -2.05030781563963,
                0.498135564435914,
                0.448469069614306,
                -1.48387807160556,
                -0.427863231850421,
                -0.0748233649223828,
                -0.0509841237610558,
                1.5525491646626
            ],
        ];
        let d = array![
            [1.0, 0.0, 0.0],
            [1.0, 0.0, 0.0],
            [1.0, 0.0, 0.0],
            [0.0, 1.0, 0.0],
            [0.0, 1.0, 0.0],
            [0.0, 1.0, 0.0],
            [0.0, 0.0, 1.0],
            [0.0, 0.0, 1.0],
            [0.0, 0.0, 1.0],
        ];
        let fit = fit_of(&y, &d);
        let none = array![
            [0i8, 0, 0],
            [0, 1, -1],
            [0, 0, 0],
            [0, -1, 0],
            [1, 0, -1],
            [0, 0, 0],
            [0, 0, 0],
            [-1, 1, 0],
            [0, 0, 0],
            [0, 0, 0],
        ];
        let bh = array![
            [0i8, 0, 0],
            [0, 1, -1],
            [0, 0, 0],
            [0, 0, 0],
            [1, 0, -1],
            [0, 0, 0],
            [0, 0, 0],
            [-1, 1, 0],
            [0, 0, 0],
            [0, 0, 0],
        ];
        assert_eq!(nestedf(&fit, Adjust::None), none);
        assert_eq!(nestedf(&fit, Adjust::BH), bh);
        // bonferroni/holm/BY coincide with BH for this dataset.
        assert_eq!(nestedf(&fit, Adjust::Bonferroni), bh);
        assert_eq!(nestedf(&fit, Adjust::Holm), bh);
        assert_eq!(nestedf(&fit, Adjust::BY), bh);
    }

    // Case 2: treatment-contrast design -> non-diagonal cov.coefficients (so
    // cov2cor + eigen-whitening run for real) with heteroscedastic genes giving
    // a finite df.prior; the five adjust methods split into three patterns.
    #[test]
    fn nestedf_treatment_finite_df_matches_r() {
        let y = array![
            [
                -0.56589213263584,
                -0.129111406171985,
                0.883844216372846,
                0.165824405266157,
                0.510226918101866,
                -0.296387565274696,
                -0.492528835489312,
                -0.648942826369362,
                0.112440116314965
            ],
            [
                -0.440493509393273,
                1.97452284169111,
                0.943716053263893,
                1.85296247687763,
                0.918879051549637,
                -1.08806877466116,
                0.6005901507507,
                2.44390249684897,
                0.758510828034637
            ],
            [
                -0.672800308611205,
                -1.12260239699623,
                -1.571493652712,
                0.758358865990787,
                3.16762861462062,
                -1.71682047480086,
                0.335444413313772,
                -3.25946041926619,
                -4.68105427064873
            ],
            [
                -0.423524520837027,
                0.478683534928917,
                -0.705633425785031,
                0.125863310320817,
                -0.90874639309216,
                -1.00280943895942,
                -0.501117549685322,
                0.366958043550944,
                -0.664545331122461
            ],
            [
                -0.14602743813943,
                1.44508592518769,
                -0.358918727489447,
                -0.66242135593651,
                0.0751840739467855,
                -0.0948756636231239,
                0.210626203050984,
                -0.501720689149703,
                -0.550236280526646
            ],
            [
                -2.86109866273184,
                0.678893878924788,
                0.991825247632336,
                -1.79353490652056,
                -0.455885146494197,
                -0.173438703743574,
                -1.08355120304816,
                -1.76993918388511,
                1.66298603162855
            ],
            [
                -0.388300259051768,
                0.440476081563293,
                -0.429356748501314,
                3.60592978345959,
                3.11410184476029,
                3.47987710551828,
                -2.49264794186162,
                -4.16894137418337,
                -3.0883195139039
            ],
            [
                -1.93331260617645,
                0.903537424251507,
                0.51244916074171,
                0.0196320667605606,
                -0.977619008137272,
                -0.288571232816508,
                0.452976738959059,
                -0.486617267289267,
                0.617269628542488
            ],
            [
                0.626091910478746,
                1.49985602027188,
                -3.49951251341427,
                4.4610353076871,
                -0.860897362008685,
                -0.326803782997329,
                -0.412492792768815,
                2.55297326743676,
                2.42184097256313
            ],
            [
                0.234326050551141,
                0.17952148545252,
                -0.36977969915983,
                0.254021081017535,
                0.453273832028051,
                -1.04934374531001,
                0.609930131193006,
                0.760402655291936,
                0.446361110207894
            ],
            [
                -0.138440069780317,
                0.252561692089675,
                1.83598463357393,
                3.70833902421073,
                3.23737777957835,
                4.13127266288387,
                -3.98118598742833,
                -3.18961650350616,
                -4.23194940211263
            ],
            [
                2.35976493177311,
                1.15449268217921,
                1.19387673975458,
                0.180160065150074,
                -1.36989243177605,
                2.64622305514602,
                -0.0689664241349938,
                -0.929549543153579,
                -1.06316624908704
            ],
            [
                -0.387941487034651,
                0.209175055219415,
                -0.0544107817969502,
                0.0691119664745673,
                0.0142895174023705,
                0.0671937469525734,
                -0.462459018396492,
                0.417995706160492,
                -0.31269377190439
            ],
            [
                -0.107005800520593,
                -0.661074953263737,
                0.516719009994363,
                3.46523024323356,
                2.90939820969627,
                3.58148939547231,
                -3.04144579252648,
                -4.38912509195675,
                -3.23472505001669
            ],
            [
                3.33934094640944,
                5.49246633658689,
                2.17312161159718,
                -0.214606718731824,
                -0.198951093031657,
                -1.06259138423522,
                -0.670897848735379,
                -3.35869824696667,
                0.697941907044954
            ],
        ];
        let d = array![
            [1.0, 0.0, 0.0],
            [1.0, 0.0, 0.0],
            [1.0, 0.0, 0.0],
            [1.0, 1.0, 0.0],
            [1.0, 1.0, 0.0],
            [1.0, 1.0, 0.0],
            [1.0, 0.0, 1.0],
            [1.0, 0.0, 1.0],
            [1.0, 0.0, 1.0],
        ];
        let fit = fit_of(&y, &d);
        let none_bh = array![
            [0i8, 0, 0],
            [0, 0, 0],
            [0, 0, 0],
            [0, 0, 0],
            [0, 0, 0],
            [0, 0, 0],
            [0, 1, -1],
            [0, 0, 0],
            [0, 0, 0],
            [0, 0, 0],
            [1, 1, -1],
            [0, 0, 0],
            [0, 0, 0],
            [0, 1, -1],
            [1, -1, -1],
        ];
        let bonf = array![
            [0i8, 0, 0],
            [0, 0, 0],
            [0, 0, 0],
            [0, 0, 0],
            [0, 0, 0],
            [0, 0, 0],
            [0, 1, -1],
            [0, 0, 0],
            [0, 0, 0],
            [0, 0, 0],
            [0, 1, -1],
            [0, 0, 0],
            [0, 0, 0],
            [0, 1, -1],
            [0, 0, 0],
        ];
        let holm_by = array![
            [0i8, 0, 0],
            [0, 0, 0],
            [0, 0, 0],
            [0, 0, 0],
            [0, 0, 0],
            [0, 0, 0],
            [0, 1, -1],
            [0, 0, 0],
            [0, 0, 0],
            [0, 0, 0],
            [0, 1, -1],
            [0, 0, 0],
            [0, 0, 0],
            [0, 1, -1],
            [1, 0, 0],
        ];
        assert_eq!(nestedf(&fit, Adjust::None), none_bh);
        assert_eq!(nestedf(&fit, Adjust::BH), none_bh);
        assert_eq!(nestedf(&fit, Adjust::Bonferroni), bonf);
        assert_eq!(nestedf(&fit, Adjust::Holm), holm_by);
        assert_eq!(nestedf(&fit, Adjust::BY), holm_by);
    }

    fn pdata() -> (Array2<f64>, Array2<f64>) {
        let p = array![
            [0.0008, 0.50, 0.30],
            [0.0001, 0.002, 0.80],
            [0.20, 0.40, 0.95],
            [0.030, 0.045, 0.60],
            [0.70, 0.65, 0.55],
            [0.004, 0.009, 0.012],
            [0.50, 0.0006, 0.40],
            [0.10, 0.11, 0.0009],
        ];
        let coef = array![
            [1.5, -0.3, 0.8],
            [-2.1, 1.7, -0.2],
            [0.4, -0.6, 0.1],
            [1.2, 1.1, -0.5],
            [-0.2, 0.3, -0.4],
            [2.0, -1.8, 1.6],
            [-0.5, -2.4, 0.6],
            [0.9, 1.0, -2.2],
        ];
        (p, coef)
    }

    #[test]
    fn decide_tests_pvalues_matches_r() {
        let (p, coef) = pdata();
        let dt = |m, a| {
            decide_tests_pvalues(&p, m, a, 0.05, 0.0, Some(&coef))
                .unwrap()
                .data
        };
        let sep_none = array![
            [1i8, 0, 0],
            [-1, 1, 0],
            [0, 0, 0],
            [1, 1, 0],
            [0, 0, 0],
            [1, -1, 1],
            [0, -1, 0],
            [0, 0, -1],
        ];
        let sep_bh = array![
            [1i8, 0, 0],
            [-1, 1, 0],
            [0, 0, 0],
            [0, 0, 0],
            [0, 0, 0],
            [1, -1, 1],
            [0, -1, 0],
            [0, 0, -1],
        ];
        let sep_bonf = array![
            [1i8, 0, 0],
            [-1, 1, 0],
            [0, 0, 0],
            [0, 0, 0],
            [0, 0, 0],
            [1, 0, 0],
            [0, -1, 0],
            [0, 0, -1],
        ];
        let glob_bonf = array![
            [1i8, 0, 0],
            [-1, 1, 0],
            [0, 0, 0],
            [0, 0, 0],
            [0, 0, 0],
            [0, 0, 0],
            [0, -1, 0],
            [0, 0, -1],
        ];
        let hier_none = sep_bh.clone(); // row4 not gene-wise significant
        let hier_bonf = glob_bonf.clone();
        let hier_holm = sep_bonf.clone();

        // separate
        assert_eq!(dt(DecideMethod::Separate, Adjust::None), sep_none);
        assert_eq!(dt(DecideMethod::Separate, Adjust::BH), sep_bh);
        assert_eq!(dt(DecideMethod::Separate, Adjust::Bonferroni), sep_bonf);
        assert_eq!(dt(DecideMethod::Separate, Adjust::Holm), sep_bonf);
        assert_eq!(dt(DecideMethod::Separate, Adjust::BY), sep_bonf);
        // global (none/BH coincide with separate; multiplicity is matrix-wide)
        assert_eq!(dt(DecideMethod::Global, Adjust::None), sep_none);
        assert_eq!(dt(DecideMethod::Global, Adjust::BH), sep_bh);
        assert_eq!(dt(DecideMethod::Global, Adjust::Bonferroni), glob_bonf);
        assert_eq!(dt(DecideMethod::Global, Adjust::Holm), glob_bonf);
        assert_eq!(dt(DecideMethod::Global, Adjust::BY), glob_bonf);
        // hierarchical (Simes)
        assert_eq!(dt(DecideMethod::Hierarchical, Adjust::None), hier_none);
        assert_eq!(dt(DecideMethod::Hierarchical, Adjust::BH), hier_none);
        assert_eq!(
            dt(DecideMethod::Hierarchical, Adjust::Bonferroni),
            hier_bonf
        );
        assert_eq!(dt(DecideMethod::Hierarchical, Adjust::Holm), hier_holm);
        assert_eq!(dt(DecideMethod::Hierarchical, Adjust::BY), hier_bonf);
    }

    #[test]
    fn decide_tests_pvalues_lfc_and_nocoef_match_r() {
        let (p, coef) = pdata();
        // lfc = 1.3 zeroes entries with |coef| < 1.3 (row4's 1.2 and 1.1).
        let lfc_res = array![
            [1i8, 0, 0],
            [-1, 1, 0],
            [0, 0, 0],
            [0, 0, 0],
            [0, 0, 0],
            [1, -1, 1],
            [0, -1, 0],
            [0, 0, -1],
        ];
        let got = decide_tests_pvalues(
            &p,
            DecideMethod::Separate,
            Adjust::None,
            0.05,
            1.3,
            Some(&coef),
        )
        .unwrap()
        .data;
        assert_eq!(got, lfc_res);

        // No coefficients: significant entries are 1 (no direction).
        let nocoef = array![
            [1i8, 0, 0],
            [1, 1, 0],
            [0, 0, 0],
            [0, 0, 0],
            [0, 0, 0],
            [1, 1, 1],
            [0, 1, 0],
            [0, 0, 1],
        ];
        let got = decide_tests_pvalues(&p, DecideMethod::Global, Adjust::BH, 0.05, 0.0, None)
            .unwrap()
            .data;
        assert_eq!(got, nocoef);
    }

    #[test]
    fn nestedf_on_pvalue_matrix_errors() {
        let (p, _) = pdata();
        assert!(
            decide_tests_pvalues(&p, DecideMethod::NestedF, Adjust::BH, 0.05, 0.0, None).is_err()
        );
    }

    #[test]
    fn bh_matches_toptable_impl() {
        let p = [0.01, 0.04, 0.03, 0.20, 0.005, 0.9, 0.5];
        let a = p_adjust(&p, Adjust::BH);
        let b = p_adjust_bh(&p);
        for (x, y) in a.iter().zip(b.iter()) {
            assert!((x - y).abs() < 1e-15, "{x} vs {y}");
        }
    }

    #[test]
    fn bonferroni_clamps_at_one() {
        let p = [0.1, 0.5, 0.9];
        let a = p_adjust(&p, Adjust::Bonferroni);
        assert!((a[0] - 0.3).abs() < 1e-15);
        assert!((a[1] - 1.0).abs() < 1e-15);
        assert!((a[2] - 1.0).abs() < 1e-15);
    }

    #[test]
    fn holm_matches_r() {
        // R: p.adjust(c(0.01,0.04,0.03,0.005), "holm") =
        //    0.03 0.06 0.06 0.02
        let p = [0.01, 0.04, 0.03, 0.005];
        let a = p_adjust(&p, Adjust::Holm);
        let want = [0.03, 0.06, 0.06, 0.02];
        for (x, y) in a.iter().zip(want.iter()) {
            assert!((x - y).abs() < 1e-12, "{x} vs {y}");
        }
    }

    #[test]
    fn by_is_bh_times_harmonic() {
        // BY(p) relates to BH via the harmonic factor q = sum(1/k).
        let p = [0.01, 0.04, 0.03, 0.20, 0.005];
        let by = p_adjust(&p, Adjust::BY);
        // R: p.adjust(c(0.01,0.04,0.03,0.20,0.005), "BY") =
        //    0.05708333 0.11416667 0.11416667 0.45666667 0.05708333
        let want = [
            0.0570833333,
            0.1141666667,
            0.1141666667,
            0.4566666667,
            0.0570833333,
        ];
        for (x, y) in by.iter().zip(want.iter()) {
            assert!((x - y).abs() < 1e-9, "{x} vs {y}");
        }
    }

    #[test]
    fn classify_tests_p_matches_r() {
        let tstat = array![
            [3.2, -1.1, 0.4],
            [2.6, 2.9, -2.7],
            [-0.3, 0.8, 1.0],
            [4.1, -3.8, 2.2],
        ];
        // limma:::.classifyTestsP(tstat, df, p.value, method); negative zeros in
        // R collapse to 0 here.
        let holm_inf = array![[1, 0, 0], [1, 1, -1], [0, 0, 0], [1, -1, 1]];
        assert_eq!(
            classify_tests_p(&tstat, &[f64::INFINITY], 0.05, Adjust::Holm),
            holm_inf
        );
        assert_eq!(
            classify_tests_p(&tstat, &[f64::INFINITY], 0.10, Adjust::BH),
            holm_inf
        );
        let df10 = array![[1, 0, 0], [1, 1, -1], [0, 0, 0], [1, -1, 0]];
        assert_eq!(classify_tests_p(&tstat, &[10.0], 0.05, Adjust::Holm), df10);
        assert_eq!(classify_tests_p(&tstat, &[5.0], 0.05, Adjust::None), df10);
    }
}