bixverse-rs 0.4.5

//! Statistical helpers

use faer::{Mat, linalg::solvers::DenseSolveCore};
use rand::rngs::StdRng;
use rand::{Rng, SeedableRng};
use rayon::prelude::*;
use statrs::distribution::FisherSnedecor;
use statrs::distribution::{Continuous, ContinuousCDF, Normal};
use statrs::function::gamma::ln_gamma;

use crate::core::math::vector_helpers::*;
use crate::prelude::*;

//////////////////
// Effect sizes //
//////////////////

/// A type alias representing effect size results
///
/// ### Fields
///
/// * `0` - The calculated effect sizes
/// * `1` - The corresponding standard errors
pub type EffectSizeRes<T> = (Vec<T>, Vec<T>);

/// Calculate the Hedge's g effect size and its standard error
///
/// ### Params
///
/// * `mean_a` - The mean values of group a.
/// * `mean_b` - The mean values of group b.
/// * `std_a` - The standard deviations of group a.
/// * `std_b` - The standard deviations of group b.
/// * `n_a` - Number of samples in a.
/// * `n_b` - Number of samples in b.
/// * `small_sample_correction` - Apply a small sample correction? Recommended
///   when `n_a` + `n_b` ≤ 35.
///
/// ### Returns
///
/// A tuple with the effect sizes being the first element, and the standard
/// errors the second element.
pub fn hedge_g_effect<T>(
    mean_a: &[T],
    mean_b: &[T],
    std_a: &[T],
    std_b: &[T],
    n_a: usize,
    n_b: usize,
    small_sample_correction: bool,
) -> EffectSizeRes<T>
where
    T: BixverseFloat,
{
    assert_same_len!(mean_a, mean_b, std_a, std_b);

    let n_a_t = T::from_usize(n_a).unwrap();
    let n_b_t = T::from_usize(n_b).unwrap();
    let total_n = T::from_usize(n_a + n_b).unwrap();
    let two = T::from_usize(2).unwrap();
    let three = T::from_usize(3).unwrap();

    let (effect_sizes, standard_errors): (Vec<T>, Vec<T>) = mean_a
        .par_iter()
        .zip(mean_b.par_iter())
        .zip(std_a.par_iter())
        .zip(std_b.par_iter())
        .map(|(((mean_a, mean_b), std_a), std_b)| {
            let pooled_sd = (((n_a_t - T::one()) * std_a.powi(2)
                + (n_b_t - T::one()) * std_b.powi(2))
                / (total_n - two))
                .sqrt();

            let mut effect_size = (*mean_a - *mean_b) / pooled_sd;

            if small_sample_correction {
                let correction_factor = ((total_n - three)
                    / (total_n - T::from_f64(2.25).unwrap()))
                    * ((total_n - two) / total_n).sqrt();
                effect_size = correction_factor * effect_size;
            }

            let standard_error =
                ((total_n / (n_a_t * n_b_t)) + (effect_size.powi(2) / (two * total_n))).sqrt();

            (effect_size, standard_error)
        })
        .unzip();

    (effect_sizes, standard_errors)
}

///////////////////////
// Statistical tests //
///////////////////////

/// Test alternatives for different statistical tests
#[derive(Clone, Debug, Default)]
pub enum TestAlternative {
    /// Two sided test for the Z-score
    #[default]
    TwoSided,
    /// One-sided test for greater than
    Greater,
    /// One-sided test for lesser than
    Less,
}

/// Helper function to get the test alternative
///
/// ### Params
///
/// * `s` - String, type of test to run.
///
/// ### Returns
///
/// Option of the `TestAlternative`
pub fn get_test_alternative(s: &str) -> Option<TestAlternative> {
    match s.to_lowercase().as_str() {
        "twosided" => Some(TestAlternative::TwoSided),
        "greater" => Some(TestAlternative::Greater),
        "less" => Some(TestAlternative::Less),
        _ => None,
    }
}

/// Transform Z-scores into p-values (assuming normality).
///
/// ### Params
///
/// * `z_scores` - The Z scores to transform to p-values
///
/// ### Returns
///
/// The p-value vector based on the Z scores (two sided)
pub fn z_scores_to_pval<T>(z_scores: &[T], test_alternative: &str) -> Vec<T>
where
    T: BixverseFloat,
{
    let test_alternative = get_test_alternative(test_alternative).unwrap_or_default();

    let normal = Normal::new(0.0, 1.0).unwrap();

    let one = T::one();
    let two = T::from_usize(2).unwrap();
    let six = T::from_usize(6).unwrap();

    z_scores
        .iter()
        .map(|&z| match test_alternative {
            TestAlternative::TwoSided => {
                let abs_z = z.abs();
                if abs_z > six {
                    let abs_z_f64 = abs_z.to_f64().unwrap();
                    let pdf = T::from_f64(normal.pdf(abs_z_f64)).unwrap();
                    let p = pdf / abs_z * (one - one / (abs_z * abs_z));
                    two * p
                } else {
                    let abs_z_f64 = abs_z.to_f64().unwrap();
                    let cdf = T::from_f64(normal.cdf(abs_z_f64)).unwrap();
                    two * (one - cdf)
                }
            }
            TestAlternative::Greater => {
                if z > six {
                    let z_f64 = z.to_f64().unwrap();
                    let pdf = T::from_f64(normal.pdf(z_f64)).unwrap();
                    pdf / z * (one - one / (z * z))
                } else {
                    let z_f64 = z.to_f64().unwrap();
                    let cdf = T::from_f64(normal.cdf(z_f64)).unwrap();
                    one - cdf
                }
            }
            TestAlternative::Less => {
                let neg_six = -six;
                if z < neg_six {
                    let abs_z = z.abs();
                    let abs_z_f64 = abs_z.to_f64().unwrap();
                    let pdf = T::from_f64(normal.pdf(abs_z_f64)).unwrap();
                    pdf / abs_z * (one - one / (abs_z * abs_z))
                } else {
                    let z_f64 = z.to_f64().unwrap();
                    T::from_f64(normal.cdf(z_f64)).unwrap()
                }
            }
        })
        .collect()
}

/// Calculate the p-value of a hypergeometric test.
///
/// ### Params
///
/// * `q` - Number of white balls drawn
/// * `m` - Number of white balls in the urn
/// * `n` - Number of black balls in the urn
/// * `k` - Number of balls drawn from the urn
///
/// ### Return
///
/// The p-value of the hypergeometric test
pub fn hypergeom_pval<T>(q: usize, m: usize, n: usize, k: usize) -> T
where
    T: BixverseFloat,
{
    // `q == 0` is P(X >= 1), a legitimate query, and the loop below computes it.
    // Returning 1.0 here meant every gene set with exactly one hit scored p = 1.
    // An empty summation range is still handled, at the `log_probs.is_empty()`
    // check further down.
    let population = m + n;
    let (n_f, m_f, k_f) = (
        T::from_usize(n).unwrap(),
        T::from_usize(m).unwrap(),
        T::from_usize(k).unwrap(),
    );
    let population_f = T::from_usize(population).unwrap();

    let upper = k.min(m);
    let mut log_probs = Vec::new();

    for i in (q + 1)..=upper {
        let i_f = T::from_usize(i).unwrap();

        // ln_gamma likely only supports f64, so convert
        let log_pmf_f64 = ln_gamma(m_f.to_f64().unwrap() + 1.0)
            - ln_gamma(i_f.to_f64().unwrap() + 1.0)
            - ln_gamma((m_f - i_f).to_f64().unwrap() + 1.0)
            + ln_gamma(n_f.to_f64().unwrap() + 1.0)
            - ln_gamma((k_f - i_f).to_f64().unwrap() + 1.0)
            - ln_gamma((n_f - (k_f - i_f)).to_f64().unwrap() + 1.0)
            - (ln_gamma(population_f.to_f64().unwrap() + 1.0)
                - ln_gamma(k_f.to_f64().unwrap() + 1.0)
                - ln_gamma((population_f - k_f).to_f64().unwrap() + 1.0));

        log_probs.push(T::from_f64(log_pmf_f64).unwrap());
    }

    if log_probs.is_empty() {
        return T::zero();
    }

    let max_log_prob = log_probs
        .iter()
        .cloned()
        .fold(T::neg_infinity(), |a, b| a.max(b));

    let mut sum = T::zero();
    for log_p in log_probs {
        sum += (log_p - max_log_prob).exp();
    }

    sum * max_log_prob.exp()
}

/// Calculate the FDR
///
/// ### Params
///
/// * `pvals` - P-values for which to calculate the FDR
///
/// ### Returns
///
/// The calculated FDRs
pub fn calc_fdr<T>(pvals: &[T]) -> Vec<T>
where
    T: BixverseFloat,
{
    let n = pvals.len();
    // The monotonicity pass below indexes `[n - 1]`, so an empty input panicked.
    // Reachable from `hypergeom_helper` with no gene sets and from
    // `finalise_go_res` when every level is empty.
    if n == 0 {
        return Vec::new();
    }
    let n_t = T::from_usize(n).unwrap();
    let one = T::one();

    let mut indexed_pval: Vec<(usize, T)> =
        pvals.par_iter().enumerate().map(|(i, &x)| (i, x)).collect();

    indexed_pval
        .sort_unstable_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal));

    let adj_pvals_tmp: Vec<T> = indexed_pval
        .par_iter()
        .enumerate()
        .map(|(i, (_, p))| {
            let i_t = T::from_usize(i + 1).unwrap();
            (n_t / i_t) * *p
        })
        .collect();

    let mut current_min = adj_pvals_tmp[n - 1].min(one);
    let mut monotonic_adj = vec![current_min; n];

    for i in (0..n - 1).rev() {
        current_min = current_min.min(adj_pvals_tmp[i]).min(one);
        monotonic_adj[i] = current_min;
    }

    let mut adj_pvals = vec![T::zero(); n];

    for (i, &(original_idx, _)) in indexed_pval.iter().enumerate() {
        adj_pvals[original_idx] = monotonic_adj[i];
    }

    adj_pvals
}

////////////
// MANOVA //
////////////

/// ManovaResults
#[derive(Debug, Clone)]
#[allow(dead_code)]
pub struct ManovaResult<T>
where
    T: BixverseFloat,
{
    /// Between-groups SSCP matrix
    pub sscp_between: Mat<T>,
    /// Within-groups SSCP matrix
    pub sscp_within: Mat<T>,
    /// Total SSCP matrix
    pub sscp_total: Mat<T>,
    /// Degrees of freedom between groups
    pub df_between: usize,
    /// Degrees of freedom within groups
    pub df_within: usize,
    /// Total degrees of freedom
    pub df_total: usize,
    /// Number of variables
    pub n_vars: usize,
    /// Means for each group
    pub group_means: Vec<Vec<T>>,
    /// Overall means
    pub overall_mean: Vec<T>,
}

#[allow(dead_code)]
impl<T> ManovaResult<T>
where
    T: BixverseFloat + std::iter::Sum,
{
    /// Non-zero eigenvalue of E⁻¹H for two-group MANOVA.
    ///
    /// For rank-1 H, the unique non-zero generalised eigenvalue equals
    /// trace(`E^-1 x H`), avoiding `det(E) / det(E+H)` which overflows f64 for
    /// moderately large `p (entries of E scale ~ n^3 / 12` for ranked data).
    ///
    /// ### Returns
    ///
    /// The eigenvalue as `T`, clamped to zero to handle numerical noise.
    fn rank1_eigenvalue(&self) -> T {
        debug_assert_eq!(
            self.df_between, 1,
            "rank1_eigenvalue assumes two-group MANOVA (df_between == 1)"
        );
        let e_inv = self.sscp_within.partial_piv_lu().inverse();
        let prod = &e_inv * &self.sscp_between;
        let trace: T = prod.diagonal().column_vector().iter().copied().sum::<T>();
        // Numerical noise can make this slightly negative for ill-conditioned E.
        trace.max(T::zero())
    }

    /// Wilks' Lambda for two-group MANOVA.
    ///
    /// Λ = 1 / (1 + λ), where λ is the rank-1 eigenvalue of E⁻¹H.
    ///
    /// ### Returns
    ///
    /// Wilks' Λ as `T`.
    pub fn wilks_lambda(&self) -> T {
        let lambda = self.rank1_eigenvalue();
        T::one() / (T::one() + lambda)
    }

    /// Pillai's trace for two-group MANOVA.
    ///
    /// V = λ / (1 + λ), where λ is the rank-1 eigenvalue of `E^-1 x H`.
    ///
    /// ### Returns
    ///
    /// Pillai's V as `T`.
    pub fn pillai_trace(&self) -> T {
        let lambda = self.rank1_eigenvalue();
        lambda / (T::one() + lambda)
    }

    /// Exact F-test for two-group MANOVA (Hotelling's T^2 form).
    ///
    /// F = ((n - p - 1) / p) * λ,   df = (p, n - p - 1)
    ///
    /// Returns (NaN, 1.0) when the inputs are degenerate
    /// (n - p - 1 <= 0, non-finite eigenvalue, etc.).
    fn two_group_f_test(&self) -> (T, T) {
        let lambda = self.rank1_eigenvalue();
        let p = T::from_usize(self.n_vars).unwrap();
        let n = T::from_usize(self.df_within + self.df_between + 1).unwrap();
        let df2 = n - p - T::one();

        if !lambda.is_finite() || df2 <= T::zero() {
            return (T::nan(), T::one());
        }

        let f_stat = (df2 / p) * lambda;

        let df1_f64 = p.to_f64().unwrap();
        let df2_f64 = df2.to_f64().unwrap();
        let f_f64 = f_stat.to_f64().unwrap();

        if !f_f64.is_finite() || f_f64 < 0.0 {
            return (f_stat, T::one());
        }

        let f_dist = match FisherSnedecor::new(df1_f64, df2_f64) {
            Ok(d) => d,
            Err(_) => return (f_stat, T::one()),
        };
        let p_value = T::from_f64(1.0 - f_dist.cdf(f_f64)).unwrap();
        (f_stat, p_value)
    }

    /// F-statistic and p-value derived from Wilks' Lambda.
    ///
    /// For two-group MANOVA this coincides with the Pillai F-test.
    ///
    /// ### Returns
    ///
    /// A tuple of `(f_stat, p_value)` as `(T, T)`.
    pub fn wilks_f_test(&self) -> (T, T) {
        self.two_group_f_test()
    }

    /// F-statistic and p-value derived from Pillai's trace.
    ///
    /// For two-group MANOVA this coincides with the Wilks F-test.
    ///
    /// ### Returns
    ///
    /// A tuple of `(f_stat, p_value)` as `(T, T)`.
    pub fn pillai_f_test(&self) -> (T, T) {
        self.two_group_f_test()
    }
}

/// ManovaSummary
#[derive(Debug)]
#[allow(dead_code)]
pub struct ManovaSummary<T>
where
    T: BixverseFloat,
{
    /// Wilks' lambda value
    pub wilks_lambda: T,
    /// Pillai's trace value
    pub pillai_trace: T,
    /// Degrees of freedom between groups
    pub df_between: usize,
    /// Degrees of freedom within groups
    pub df_within: usize,
    /// F statistic according to Wilk
    pub f_stat_wilk: T,
    /// P-value according to Wilk
    pub p_val_wilk: T,
    /// F statistic according to Pillai
    pub f_stat_pillai: T,
    /// P-value according to Pillai
    pub p_val_pillai: T,
}

impl<T> ManovaSummary<T>
where
    T: BixverseFloat + std::iter::Sum,
{
    /// Get the summary results from a `ManovaRes`
    ///
    /// ### Params
    ///
    /// * `res` - The calculated ManovaResults
    ///
    /// ### Returns
    ///
    /// The `ManovaSummary`.
    pub fn from_manova_res(res: &ManovaResult<T>) -> Self {
        let (f_stat_wilk, p_val_wilk) = res.wilks_f_test();
        let (f_stat_pillai, p_val_pillai) = res.pillai_f_test();

        ManovaSummary {
            wilks_lambda: res.wilks_lambda(),
            pillai_trace: res.pillai_trace(),
            df_between: res.df_between,
            df_within: res.df_within,
            f_stat_wilk,
            p_val_wilk,
            f_stat_pillai,
            p_val_pillai,
        }
    }
}

///////////
// ANOVA //
///////////

/// AnovaSummary (based on MANOVA models)
#[derive(Debug, Clone)]
#[allow(dead_code)]
pub struct AnovaSummary<T>
where
    T: BixverseFloat,
{
    /// Variable index
    pub variable_index: usize,
    /// Sum of squares between groups
    pub ss_between: T,
    /// Sum of squares within groups
    pub ss_within: T,
    /// Mean square between groups
    pub ms_between: T,
    /// Mean square within groups
    pub ms_within: T,
    /// F statistic
    pub f_stat: T,
    /// P-value
    pub p_val: T,
}

/// Generates from MANOVE results the AnovaSummary
///
/// ### Params
///
/// * `res` - The MANOVA result to analyse
///
/// ### Returns
///
/// A vector of AnovaSummaries
pub fn summary_aov<T>(res: &ManovaResult<T>) -> Vec<AnovaSummary<T>>
where
    T: BixverseFloat,
{
    let mut aov_res = Vec::with_capacity(res.n_vars);
    let df_between_t = T::from_usize(res.df_between).unwrap();
    let df_within_t = T::from_usize(res.df_within).unwrap();

    let f_dist = FisherSnedecor::new(
        df_between_t.to_f64().unwrap(),
        df_within_t.to_f64().unwrap(),
    )
    .ok();

    for var_idx in 0..res.n_vars {
        let ss_between = res.sscp_between[(var_idx, var_idx)];
        let ss_within = res.sscp_within[(var_idx, var_idx)];
        let ms_between = ss_between / df_between_t;
        let ms_within = ss_within / df_within_t;

        let (f_stat, p_val) = if ms_within <= T::zero() {
            (T::nan(), T::one())
        } else {
            let f = ms_between / ms_within;
            let f_f64 = f.to_f64().unwrap();
            let pv = match (&f_dist, f_f64.is_finite() && f_f64 >= 0.0) {
                (Some(d), true) => T::from_f64(1.0 - d.cdf(f_f64)).unwrap(),
                _ => T::one(),
            };
            (f, pv)
        };

        aov_res.push(AnovaSummary {
            variable_index: var_idx,
            ss_between,
            ss_within,
            ms_between,
            ms_within,
            f_stat,
            p_val,
        });
    }

    aov_res
}

///////////////////
// Probabilities //
///////////////////

/// Implementation of the trigamma function (second derivative of ln(gamma(x)))
///
/// ### Params
///
/// * `x` - The value for which to calculate the trigamma function.
///
/// ### Returns
///
/// The trigamma value for the given input.
pub fn trigamma<T: BixverseFloat>(x: T) -> T {
    let mut x = x;
    let mut result = T::zero();

    if x <= T::from_f64(5.0).unwrap() {
        while x < T::from_f64(5.0).unwrap() {
            result += T::one() / (x * x);
            x += T::one();
        }
    }

    let xx = x * x;
    result += T::one() / x
        + T::one() / (T::from_f64(2.0).unwrap() * xx)
        + T::one() / (T::from_f64(6.0).unwrap() * xx * x);

    let xxx = xx * x;
    result += -T::one() / (T::from_f64(30.0).unwrap() * xxx * x)
        + T::one() / (T::from_f64(42.0).unwrap() * xxx * xx * x)
        - T::one() / (T::from_f64(30.0).unwrap() * xxx * xxx * x);

    result
}

///////////
// Other //
///////////

/// Logit function
///
/// ### Params
///
/// * `p` - Probability value (must be in (0, 1))
///
/// ### Returns
///
/// Log-odds: ln(p / (1-p))
pub fn logit<T>(p: T) -> T
where
    T: BixverseFloat,
{
    (p / (T::one() - p)).ln()
}

/// Inverse logit (sigmoid) function
///
/// ### Params
///
/// * `q` - Log-odds value
///
/// ### Returns
///
/// Probability: exp(q) / (1 + exp(q))
pub fn inv_logit<T>(q: T) -> T
where
    T: BixverseFloat,
{
    q.exp() / (T::one() + q.exp())
}

/////////////////////
// Critical values //
/////////////////////

/// Calculate the critical value using bootstrap resampling
///
/// ### Params
///
/// * `values` - Slice of values to resample from.
/// * `sample_size` - Number of samples to draw in the bootstrap sample.
/// * `alpha` - The significance level for the critical value.
/// * `seed` - Random seed for reproducibility.
///
/// ### Returns
///
/// The critical value at the specified alpha level.
pub fn calculate_critval<T: BixverseFloat>(
    values: &[T],
    sample_size: usize,
    alpha: &T,
    seed: usize,
) -> T {
    let mut rng = StdRng::seed_from_u64(seed as u64);
    let mut random_sample: Vec<T> = (0..sample_size)
        .map(|_| {
            let index = rng.random_range(0..values.len());
            values[index]
        })
        .collect();
    random_sample.sort_by(|a, b| b.partial_cmp(a).unwrap());
    let index = (*alpha * T::from_usize(random_sample.len()).unwrap())
        .ceil()
        .to_usize()
        .unwrap();
    random_sample[index + 1]
}

//////////////
// Outliers //
//////////////

/// Type of outlier detection for MAD thresholding
#[derive(Clone, Debug, Default)]
pub enum OutlierDirection {
    /// Check if outlier is below OR above the thresholds
    #[default]
    Both,
    /// Check if outlier is below the threshold
    Below,
    /// Check if outlier is above the threshold
    Above,
}

/// Helper function to get the outlier detection method
///
/// ### Params
///
/// * `s` - String, type of test to run.
///
/// ### Returns
///
/// Option of the `OutlierDirection`
pub fn parse_outlier_type(s: &str) -> Option<OutlierDirection> {
    match s.to_lowercase().as_str() {
        "below" => Some(OutlierDirection::Below),
        "above" => Some(OutlierDirection::Above),
        "twosided" => Some(OutlierDirection::Both),
        _ => None,
    }
}

/// MAD outlier detection
///
/// ### Params
///
/// * `x` - Slice of values to check for outliers.
/// * `threshold` - Number of MADs to accept as not being an outlier.
/// * `direction` - Direction to check for outliers (below, above, or both).
///
/// ### Returns
///
/// A tuple with a vector of booleans indicating whether each value in the input
/// is an outlier and the value of applied margin.
pub fn mad_outlier<T>(x: &[T], threshold: T, direction: OutlierDirection) -> (Vec<bool>, T)
where
    T: BixverseFloat,
{
    let median_val = match median(x) {
        Some(m) => m,
        None => return (vec![], T::zero()),
    };
    let mad_val = match mad(x, None) {
        Some(m) => m,
        None => return (vec![], T::zero()),
    };
    let margin = threshold * mad_val;
    let res = x
        .iter()
        .map(|&v| match direction {
            OutlierDirection::Below => v < median_val - margin,
            OutlierDirection::Above => v > median_val + margin,
            OutlierDirection::Both => (v - median_val).abs() > margin,
        })
        .collect::<Vec<bool>>();
    (res, margin)
}

///////////
// Tests //
///////////

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

    /// `logit` and `inv_logit` are inverses, with p = 0.5 sitting at zero log-odds.
    #[test]
    fn test_logit_and_inv_logit() {
        let p: f64 = 0.5;
        let log_odds: f64 = logit(p);
        assert!((log_odds - 0.0).abs() < 1e-6); // ln(0.5 / 0.5) = ln(1) = 0.0

        let recovered_p = inv_logit(log_odds);
        assert!((recovered_p - 0.5).abs() < 1e-6);
    }

    /// Two-sided p-values at the two z-scores whose answers are known by heart.
    #[test]
    fn test_z_scores_to_pval() {
        let z_scores: Vec<f64> = vec![0.0, 1.95996398454]; // ~ 95% confidence interval

        let pvals_two_sided = z_scores_to_pval(&z_scores, "twosided");
        // Z = 0 -> p = 1.0
        assert!((pvals_two_sided[0] - 1.0).abs() < 1e-6);
        // Z = 1.96 -> p ≈ 0.05
        assert!((pvals_two_sided[1] - 0.05).abs() < 1e-5);
    }

    /// Benjamini-Hochberg adjustment stays monotonic and comes back in the input order.
    #[test]
    fn test_calc_fdr() {
        // Classic Benjamini-Hochberg test case
        let pvals: Vec<f64> = vec![0.01, 0.04, 0.03];
        // Sorted: 0.01 (idx 0), 0.03 (idx 2), 0.04 (idx 1)
        // Adjusted tmp: 0.01*(3/1)=0.03, 0.03*(3/2)=0.045, 0.04*(3/3)=0.04
        // Monotonic min backwards: 0.04, min(0.04, 0.045)=0.04, min(0.04, 0.03)=0.03
        // Result: [0.03, 0.04, 0.04]
        let fdr = calc_fdr(&pvals);

        assert!((fdr[0] - 0.03).abs() < 1e-6);
        assert!((fdr[1] - 0.04).abs() < 1e-6);
        assert!((fdr[2] - 0.04).abs() < 1e-6);
    }

    /// Only points beyond the median plus `threshold` MADs count as outliers.
    #[test]
    fn test_mad_outlier() {
        let vec: Vec<f64> = vec![1.0, 2.0, 3.0, 4.0, 100.0];
        // Median is 3.0
        // Deviations: [2.0, 1.0, 0.0, 1.0, 97.0]
        // MAD (median of deviations) is 1.0
        // Threshold = 3.0, so margin = 3.0
        // Acceptable range: 3.0 ± 3.0 = [0.0, 6.0]
        let (outliers, margin) = mad_outlier(&vec, 3.0, OutlierDirection::Both);

        assert!((margin - 3.0).abs() < 1e-6);
        assert_eq!(outliers, vec![false, false, false, false, true]); // 100.0 is an outlier
    }

    /// The alternative and outlier-direction parsers ignore case.
    #[test]
    fn test_parse_helpers() {
        assert!(matches!(
            get_test_alternative("twosided"),
            Some(TestAlternative::TwoSided)
        ));
        assert!(matches!(
            get_test_alternative("GREATER"),
            Some(TestAlternative::Greater)
        ));

        assert!(matches!(
            parse_outlier_type("below"),
            Some(OutlierDirection::Below)
        ));
        assert!(matches!(
            parse_outlier_type("Twosided"),
            Some(OutlierDirection::Both)
        ));
    }

    /// Two-group MANOVA against hand-computed Wilks, Pillai and F statistics.
    #[test]
    fn test_manova_two_group_analytic() {
        // Group 0: rows 0, 1; Group 1: rows 2, 3
        // x = [[1, 4], [2, 1], [5, 7], [6, 3]]
        //
        // Hand-computed:
        //   sscp_within = [[1.0, -3.5], [-3.5, 12.5]]
        //   d = mean_1 - mean_0 = [4.0, 2.5]
        //   c = n1*n2/n = 1
        //   lambda = c * d' * E^-1 * d = 1105
        //   Wilks = 1/1106, Pillai = 1105/1106, Wilks + Pillai = 1
        //   F = ((n - p - 1)/p) * lambda = (1/2)*1105 = 552.5, df = (2, 1)
        let mat = Mat::from_fn(4, 2, |i, j| match (i, j) {
            (0, 0) => 1.0f64,
            (0, 1) => 4.0,
            (1, 0) => 2.0,
            (1, 1) => 1.0,
            (2, 0) => 5.0,
            (2, 1) => 7.0,
            (3, 0) => 6.0,
            (3, 1) => 3.0,
            _ => unreachable!(),
        });
        let res = crate::enrichment::mitch::manova_mitch(mat.as_ref(), &[2, 3]);

        let wilks = res.wilks_lambda();
        let pillai = res.pillai_trace();
        let recovered_lambda = (1.0 - wilks) / wilks;

        assert!((recovered_lambda - 1105.0).abs() < 1e-6);
        assert!((wilks + pillai - 1.0).abs() < 1e-12);
        assert!((wilks - 1.0 / 1106.0).abs() < 1e-9);
        assert!((pillai - 1105.0 / 1106.0).abs() < 1e-9);

        let (f_w, _) = res.wilks_f_test();
        let (f_p, p_p) = res.pillai_f_test();
        assert!((f_w - 552.5).abs() < 1e-6);
        assert!((f_p - 552.5).abs() < 1e-6);
        assert!((0.0..=1.0).contains(&p_p));
    }

    /// Regression: at Mitch-scale magnitudes `det()` overflowed to NaN and panicked downstream.
    #[test]
    fn test_manova_no_overflow_large_p() {
        // Regression test: previously panicked via det() overflow → NaN →
        // FisherSnedecor::cdf XOutOfRange. Mimics Mitch ranked-data magnitudes
        // (n ~ 3000, sscp diagonals ~ n^3/12 ~ 2e9).
        let p = 80;
        let n_total = 2949usize;
        let sigma = (n_total as f64).powi(3) / 12.0;

        let sscp_within = Mat::from_fn(p, p, |i, j| if i == j { sigma } else { sigma * 0.3 });
        // Small between-group signal; H entries ~ sigma * 1e-3
        let sscp_between = Mat::from_fn(
            p,
            p,
            |i, j| {
                if i == j { sigma * 1.5e-3 } else { sigma * 5e-4 }
            },
        );
        let sscp_total = &sscp_within + &sscp_between;

        let res = ManovaResult::<f64> {
            sscp_between,
            sscp_within,
            sscp_total,
            df_between: 1,
            df_within: n_total - 2,
            df_total: n_total - 1,
            n_vars: p,
            group_means: vec![vec![0.0; p], vec![0.0; p]],
            overall_mean: vec![0.0; p],
        };

        let wilks = res.wilks_lambda();
        let pillai = res.pillai_trace();
        let (f_w, p_w) = res.wilks_f_test();
        let (f_p, p_p) = res.pillai_f_test();

        assert!(wilks.is_finite() && wilks > 0.0 && wilks <= 1.0);
        assert!((0.0..1.0).contains(&pillai) && pillai.is_finite());
        assert!((wilks + pillai - 1.0).abs() < 1e-9);
        assert!(f_w.is_finite() && f_w >= 0.0);
        assert!(f_p.is_finite() && f_p >= 0.0);
        assert!((0.0..=1.0).contains(&p_w));
        assert!((0.0..=1.0).contains(&p_p));
        // Wilks and Pillai F-tests are identical for two groups.
        assert!((f_w - f_p).abs() < 1e-9);
        assert!((p_w - p_p).abs() < 1e-9);
    }

    /// Regression: a constant column has no within-group variance and used to panic.
    #[test]
    fn test_summary_aov_constant_column() {
        // ss_within == 0 used to feed NaN into FisherSnedecor and panic.
        let p = 2;
        let sscp_within = Mat::from_fn(p, p, |i, j| {
            if i == 0 && j == 0 {
                0.0
            } else if i == j {
                10.0
            } else {
                0.0
            }
        });
        let sscp_between = Mat::from_fn(p, p, |i, j| if i == j { 5.0 } else { 0.0 });
        let sscp_total = &sscp_within + &sscp_between;

        let res = ManovaResult::<f64> {
            sscp_between,
            sscp_within,
            sscp_total,
            df_between: 1,
            df_within: 10,
            df_total: 11,
            n_vars: p,
            group_means: vec![vec![0.0; p], vec![0.0; p]],
            overall_mean: vec![0.0; p],
        };

        let aov = summary_aov(&res);
        assert_eq!(aov.len(), 2);
        assert!(aov[0].f_stat.is_nan());
        assert!((aov[0].p_val - 1.0).abs() < 1e-12);
        assert!(aov[1].f_stat.is_finite() && aov[1].f_stat > 0.0);
        assert!((0.0..=1.0).contains(&aov[1].p_val));
    }
}