gam-sae 0.3.149

//! Genuine FDR control for the sparse-dict / curved chart lanes via
//! universal-inference split-LR e-values + full-family e-BH (issue #2246).
//!
//! # Why the descriptive BIC margin is not an e-value
//!
//! The lanes ([`super::block_chart`], [`super::block_stream`], [`super::cofit`],
//! [`crate::manifold::inframe_curved`], [`crate::tiered`]) score a curved chart
//! against a rank-1 linear shell with a held-out **BIC margin**
//! `M = deviance_gain − ½·d_eff·ln n_eff`. A BIC/Laplace margin is NOT a valid
//! log-e-value: `E = max(exp M, 1) ≥ 1` always and `E > 1` with positive
//! probability under any nondegenerate null, so `E_{H0}[E] > 1`. No calibration
//! of a BIC margin fixes this. Feeding such margins to e-BH (as an earlier build
//! did) voids the FDR guarantee — the 2026-07-09 external math audit's finding 4.
//! Those claims were scrubbed; the margins are now documented as descriptive
//! model-selection verdicts only. This module is the designed replacement: a
//! genuinely FDR-controlled discovery gate.
//!
//! # Universal inference (Wasserman–Ramdas–Balakrishnan)
//!
//! The null (rank-1 linear shell) and the alternative (radial / curved chart)
//! are **non-nested**, so a classical likelihood-ratio test has no calibrated
//! reference distribution. Universal inference sidesteps this with a data split:
//!
//! 1. Split the candidate's rows into a fit fold `D0` and an evaluation fold
//!    `D1` (disjoint).
//! 2. Fit the alternative density `p_alt(·; θ̂_{D0})` on `D0` ONLY — so its
//!    parameter `θ̂_{D0}` is independent of `D1`.
//! 3. Form the split e-value
//!    `E = p_alt(D1; θ̂_{D0}) / sup_{θ0 ∈ H0} p_null(D1; θ0)` — the null is
//!    re-maximized (MLE) on the held-out `D1`, INCLUDING every nuisance
//!    parameter (mean, subspace, variances).
//!
//! ## Provable `E_{H0}[E] ≤ 1`
//!
//! Suppose the data are generated by some `θ0* ∈ H0`. Bounding the re-maximized
//! null denominator below by the density at the true parameter,
//! `sup_{θ0} p_null(D1; θ0) ≥ p_null(D1; θ0*)`, gives
//!
//! ```text
//!   E ≤ p_alt(D1; θ̂_{D0}) / p_null(D1; θ0*).
//! ```
//!
//! Taking the `H0` expectation and integrating over `D1` first (`θ̂_{D0}` is a
//! constant given `D0`, and both `p_alt`, `p_null` are PROPER densities that
//! integrate to 1):
//!
//! ```text
//!   E_{H0}[E] ≤ E_{D0} E_{D1|D0}[ p_alt(D1; θ̂_{D0}) / p_null(D1; θ0*) ]
//!            = E_{D0} ∫ p_alt(D1; θ̂_{D0}) dD1
//!            = 1.
//! ```
//!
//! This holds regardless of non-nestedness, of whether the alternative is
//! correctly specified, and of the quality of `θ̂_{D0}` — a bad fit only costs
//! power, never validity. The one non-negotiable requirement is that both
//! `p_alt` and `p_null` are genuine normalized densities on the same space; the
//! lane's `radial_predict` / `pca_reconstruct` reconstruction losses are
//! self-referential projections (the "prediction" of a point uses the point
//! itself) and are NOT densities, which is exactly why this module fits proper
//! parametric densities instead of reusing them: a 1-factor probabilistic-PCA
//! Gaussian for the linear-shell null and a Bessel-`I0` radial-ring marginal for
//! the curved alternative (see [`fit_ppca1`], [`fit_ring`]).
//!
//! ## Cross-fit for power
//!
//! A single split wastes half the rows on each side. We swap the roles of every
//! fold and AVERAGE the resulting e-values: `Ē = (1/K)·Σ_f E_f`. Each `E_f` is a
//! valid e-value (its alternative is fit on the complement of fold `f`, its null
//! re-maximized on fold `f`), and the mean of valid e-values is a valid e-value
//! by linearity — `E_{H0}[Ē] = (1/K)·Σ_f E_{H0}[E_f] ≤ 1` — so cross-fitting
//! buys power at zero cost to validity.
//!
//! # Full-family e-BH
//!
//! The screened candidate family (energy top-k blocks, score-thresholded pairs)
//! must FEED the e-BH family, not gate it after the evidence is seen: gating on
//! the same statistic post hoc redefines the family and voids the guarantee.
//! Every screened candidate contributes one `log_e`, and the whole vector is
//! handed to [`gam_terms::inference::structure_evidence::e_benjamini_hochberg`]
//! (Wang–Ramdas e-BH, FDR ≤ α under ARBITRARY dependence — the co-firing charts
//! share rows, so the PRDS assumption p-value BH needs is violated and only e-BH
//! is legal). There is no permutation-style budget cap `B ≥ m/α`: the UI
//! construction yields unbounded e-values, so a single strong candidate can clear
//! the `ln(m/(α·k))` bar on its own.

use gam_terms::inference::structure_evidence::e_benjamini_hochberg;
use ndarray::Array2;

use super::block_chart::jacobi_eigh;

/// The FDR-controlled discovery certificate for one screened candidate family.
#[derive(Clone, Debug)]
pub struct FdrCertificate {
    /// The declared target FDR level the certificate controls at.
    pub alpha: f64,
    /// Per-candidate universal-inference split-LR log-e-value, one entry per
    /// screened candidate in the family's order. `−∞` marks a candidate that
    /// could not support a curved claim (too few rows, or fewer than two
    /// coordinate channels) — it banks no evidence and is never rejected.
    pub log_e: Vec<f64>,
    /// Indices (into `log_e`) of the candidates e-BH confirms as genuine curved
    /// structure at FDR ≤ `alpha`.
    pub rejected: Vec<usize>,
}

/// Run full-family e-BH over a screened candidate family's log-e-values.
///
/// The caller is responsible for having computed EVERY screened candidate's
/// `log_e` (via [`crossfit_ui_log_evalue`] / [`shell_vs_ring_log_evalue`]) and
/// passing them ALL — the screening feeds this family, it must not gate it.
pub fn family_fdr_certificate(log_e: Vec<f64>, alpha: f64) -> FdrCertificate {
    let rejected = e_benjamini_hochberg(&log_e, alpha);
    FdrCertificate {
        alpha,
        log_e,
        rejected,
    }
}

/// Cross-fit universal-inference split-LR log-e-value from fitter closures.
///
/// Generic over the alternative's fitted state `A`. For each of `folds` folds,
/// rows with `i % folds == f` are the evaluation fold `D1` and the rest are the
/// fit fold `D0`. The alternative is fit on `D0` (`fit_alt`) and scored on `D1`
/// (`alt_loglik`); the null is re-maximized and scored in-sample on `D1`
/// (`null_sup_loglik`). The per-fold e-values are averaged (see module doc: the
/// mean of valid e-values is a valid e-value). Returns `ln Ē`.
///
/// - `fit_alt(train)`: fit the alternative on the train rows, returning `None`
///   if the train fold is too degenerate to support a curved claim (that fold
///   then banks the trivially-valid e-value `0`, conservatively).
/// - `alt_loglik(alt, eval)`: the held-out log-density `Σ_{i∈eval} ln p_alt(y_i;
///   θ̂_{D0})` under the `D0`-fit alternative.
/// - `null_sup_loglik(eval)`: the constrained supremum `sup_{θ0∈H0} Σ_{i∈eval}
///   ln p_null(y_i; θ0)` — the null density's own MLE on the eval fold. An
///   under-maximized null inflates the e-value and voids validity.
pub fn crossfit_ui_log_evalue<A>(
    n: usize,
    folds: usize,
    mut fit_alt: impl FnMut(&[usize]) -> Result<Option<A>, String>,
    mut alt_loglik: impl FnMut(&A, &[usize]) -> Result<f64, String>,
    mut null_sup_loglik: impl FnMut(&[usize]) -> Result<f64, String>,
) -> Result<f64, String> {
    if n < 2 {
        return Ok(f64::NEG_INFINITY);
    }
    let folds = folds.max(2).min(n);
    let mut fold_log_e: Vec<f64> = Vec::with_capacity(folds);
    for f in 0..folds {
        let mut train = Vec::new();
        let mut eval = Vec::new();
        for i in 0..n {
            if i % folds == f {
                eval.push(i);
            } else {
                train.push(i);
            }
        }
        if train.is_empty() || eval.is_empty() {
            continue;
        }
        // A fold whose train side cannot support the curved alternative banks
        // the trivially-valid e-value 0 (`ln 0 = −∞`): `E_{H0}[0] = 0 ≤ 1`, so
        // validity is preserved and the fold only withholds power.
        let log_e_f = match fit_alt(&train)? {
            Some(alt) => {
                let alt_ll = alt_loglik(&alt, &eval)?;
                let null_ll = null_sup_loglik(&eval)?;
                alt_ll - null_ll
            }
            None => f64::NEG_INFINITY,
        };
        fold_log_e.push(log_e_f);
    }
    if fold_log_e.is_empty() {
        return Ok(f64::NEG_INFINITY);
    }
    // ln Ē = ln( (1/K)·Σ_f exp(log_e_f) ) = logsumexp(log_e_f) − ln K.
    Ok(logsumexp(&fold_log_e) - (fold_log_e.len() as f64).ln())
}

/// Concrete universal-inference split-LR log-e-value for one screened candidate:
/// linear-shell null vs radial-ring alternative on the candidate's coordinate
/// rows.
///
/// `coords` is `[n_rows, q]`, the block-coordinate observations for the rows
/// where the candidate's block(s) fire. The null density is a 1-factor
/// probabilistic-PCA Gaussian ([`fit_ppca1`], the rank-1 linear shell as a
/// proper density); the alternative is the Bessel-`I0` radial-ring marginal in
/// the train data's top-2 eigen-plane ([`fit_ring`]). A candidate with fewer
/// than two coordinate channels cannot support a ring and returns `−∞`.
pub fn shell_vs_ring_log_evalue(
    coords: &Array2<f64>,
    folds: usize,
    ridge: f64,
) -> Result<f64, String> {
    let n = coords.nrows();
    let q = coords.ncols();
    if q < 2 {
        return Ok(f64::NEG_INFINITY);
    }
    crossfit_ui_log_evalue(
        n,
        folds,
        |train| {
            if train.len() < 3 {
                return Ok(None);
            }
            Ok(Some(fit_ring(coords, train, ridge)))
        },
        |ring, eval| Ok(ring_loglik(coords, eval, ring)),
        |eval| {
            let null = fit_ppca1(coords, eval, ridge)?;
            Ok(ppca1_loglik(coords, eval, &null))
        },
    )
}

/// A fitted 1-factor probabilistic-PCA density `N(μ, σ0²·I + λ0·v vᵀ)` — the
/// rank-1 linear shell null as a PROPER Gaussian density.
struct Ppca1 {
    mean: Vec<f64>,
    /// Unit top eigenvector `v`.
    axis: Vec<f64>,
    /// Isotropic residual variance `σ0²` (mean of the trailing eigenvalues).
    sigma2: f64,
    /// Factor loading variance `λ0 = max(ℓ1 − σ0², 0)`.
    lambda: f64,
    logdet: f64,
}

/// Maximum-likelihood fit of the 1-factor PPCA density on `rows`. This is the
/// honest constrained supremum over the linear-shell null: mean, principal axis,
/// factor variance, and isotropic residual variance are all re-estimated.
fn fit_ppca1(coords: &Array2<f64>, rows: &[usize], ridge: f64) -> Result<Ppca1, String> {
    let n = rows.len();
    let q = coords.ncols();
    if n == 0 || q == 0 {
        return Err("ppca1 fit: empty rows or columns".to_string());
    }
    let mut mean = vec![0.0f64; q];
    for &i in rows {
        for j in 0..q {
            mean[j] += coords[[i, j]];
        }
    }
    for m in &mut mean {
        *m /= n as f64;
    }
    let mut cov = vec![0.0f64; q * q];
    for &i in rows {
        for a in 0..q {
            let va = coords[[i, a]] - mean[a];
            for b in 0..q {
                cov[a * q + b] += va * (coords[[i, b]] - mean[b]);
            }
        }
    }
    // Biased (MLE) covariance normaliser 1/n, so the density is the genuine
    // constrained MLE. A relative ridge keeps the trailing variance strictly
    // positive on a perfectly-degenerate fold.
    for v in &mut cov {
        *v /= n as f64;
    }
    let (vals, vecs) = jacobi_eigh(cov, q)?;
    let mut order: Vec<usize> = (0..q).collect();
    order.sort_by(|&a, &b| {
        vals[b]
            .partial_cmp(&vals[a])
            .unwrap_or(std::cmp::Ordering::Equal)
    });
    let trace: f64 = vals.iter().sum();
    let floor = (ridge * trace / q as f64).max(1.0e-300);
    let top = order[0];
    let ell1 = vals[top].max(floor);
    let axis: Vec<f64> = (0..q).map(|j| vecs[j * q + top]).collect();
    // σ0² = mean of the q−1 trailing eigenvalues (the isotropic PPCA residual).
    let sigma2 = if q > 1 {
        let rest: f64 = order.iter().skip(1).map(|&k| vals[k].max(0.0)).sum();
        (rest / (q - 1) as f64).max(floor)
    } else {
        ell1
    };
    let lambda = (ell1 - sigma2).max(0.0);
    let logdet = (q as f64 - 1.0) * sigma2.ln() + (sigma2 + lambda).ln();
    Ok(Ppca1 {
        mean,
        axis,
        sigma2,
        lambda,
        logdet,
    })
}

/// In-sample total log-density `Σ_{i∈rows} ln p_null(y_i)` under a fitted PPCA1.
/// Sherman–Morrison gives the quadratic form without forming `C0⁻¹`.
fn ppca1_loglik(coords: &Array2<f64>, rows: &[usize], m: &Ppca1) -> f64 {
    let q = coords.ncols();
    let log_2pi = (2.0 * std::f64::consts::PI).ln();
    let sm = m.lambda / (m.sigma2 + m.lambda);
    let mut total = 0.0;
    for &i in rows {
        let mut norm2 = 0.0;
        let mut proj = 0.0;
        for j in 0..q {
            let r = coords[[i, j]] - m.mean[j];
            norm2 += r * r;
            proj += r * m.axis[j];
        }
        let quad = (norm2 - sm * proj * proj) / m.sigma2;
        total += -0.5 * (q as f64 * log_2pi + m.logdet + quad);
    }
    total
}

/// A fitted radial-ring density: a Bessel-`I0` ring in the top-2 eigen-plane of
/// the train fold times an isotropic Gaussian off that plane. A PROPER density.
struct Ring {
    mean: Vec<f64>,
    /// Orthonormal top-2 eigen-plane basis `(e1, e2)`.
    e1: Vec<f64>,
    e2: Vec<f64>,
    /// Ring radius `R = mean_i ‖z_i‖` (in-plane).
    radius: f64,
    /// In-plane isotropic noise variance `σ²` (radial spread about the ring).
    sigma2: f64,
    /// Off-plane isotropic variance `σ⊥²`; `None` when `q == 2`.
    sigma_perp2: Option<f64>,
}

/// Fit the radial-ring alternative on the train rows. All parameters (mean,
/// plane, radius, both variances) are functions of the train fold ONLY, so the
/// resulting density is `D0`-measurable — the predictability the UI e-value
/// needs.
fn fit_ring(coords: &Array2<f64>, rows: &[usize], ridge: f64) -> Ring {
    let n = rows.len();
    let q = coords.ncols();
    let mut mean = vec![0.0f64; q];
    for &i in rows {
        for j in 0..q {
            mean[j] += coords[[i, j]];
        }
    }
    for m in &mut mean {
        *m /= n as f64;
    }
    let mut cov = vec![0.0f64; q * q];
    for &i in rows {
        for a in 0..q {
            let va = coords[[i, a]] - mean[a];
            for b in 0..q {
                cov[a * q + b] += va * (coords[[i, b]] - mean[b]);
            }
        }
    }
    for v in &mut cov {
        *v /= n as f64;
    }
    let (vals, vecs) = match jacobi_eigh(cov, q) {
        Ok(pair) => pair,
        Err(_) => (vec![0.0; q], identity_flat(q)),
    };
    let mut order: Vec<usize> = (0..q).collect();
    order.sort_by(|&a, &b| {
        vals[b]
            .partial_cmp(&vals[a])
            .unwrap_or(std::cmp::Ordering::Equal)
    });
    let k1 = order[0];
    let k2 = order[1];
    let e1: Vec<f64> = (0..q).map(|j| vecs[j * q + k1]).collect();
    let e2: Vec<f64> = (0..q).map(|j| vecs[j * q + k2]).collect();
    let trace: f64 = vals.iter().sum();
    let floor = (ridge * trace / q as f64).max(1.0e-12);
    let mut radius = 0.0;
    let mut radial_var = 0.0;
    let mut perp_ss = 0.0;
    let mut rhos = Vec::with_capacity(n);
    for &i in rows {
        let mut z1 = 0.0;
        let mut z2 = 0.0;
        let mut norm2 = 0.0;
        for j in 0..q {
            let r = coords[[i, j]] - mean[j];
            z1 += r * e1[j];
            z2 += r * e2[j];
            norm2 += r * r;
        }
        let rho = (z1 * z1 + z2 * z2).sqrt();
        rhos.push(rho);
        radius += rho;
        perp_ss += (norm2 - z1 * z1 - z2 * z2).max(0.0);
    }
    radius /= n as f64;
    for &rho in &rhos {
        radial_var += (rho - radius) * (rho - radius);
    }
    let sigma2 = (radial_var / n as f64).max(floor);
    let sigma_perp2 = if q > 2 {
        Some((perp_ss / (n as f64 * (q - 2) as f64)).max(floor))
    } else {
        None
    };
    Ring {
        mean,
        e1,
        e2,
        radius,
        sigma2,
        sigma_perp2,
    }
}

/// In-sample total log-density `Σ_{i∈rows} ln p_alt(y_i)` under a fitted Ring.
/// In-plane: the ring marginal `(2πσ²)⁻¹ exp(−(‖z‖²+R²)/2σ²) I0(R‖z‖/σ²)`;
/// off-plane: `N(0, σ⊥²·I_{q−2})`.
fn ring_loglik(coords: &Array2<f64>, rows: &[usize], m: &Ring) -> f64 {
    let q = coords.ncols();
    let log_2pi = (2.0 * std::f64::consts::PI).ln();
    let mut total = 0.0;
    for &i in rows {
        let mut z1 = 0.0;
        let mut z2 = 0.0;
        let mut norm2 = 0.0;
        for j in 0..q {
            let r = coords[[i, j]] - m.mean[j];
            z1 += r * m.e1[j];
            z2 += r * m.e2[j];
            norm2 += r * r;
        }
        let rho = (z1 * z1 + z2 * z2).sqrt();
        let plane = -(log_2pi + m.sigma2.ln())
            - (rho * rho + m.radius * m.radius) / (2.0 * m.sigma2)
            + ln_i0(m.radius * rho / m.sigma2);
        total += plane;
        if let Some(sp) = m.sigma_perp2 {
            let perp_ss = (norm2 - z1 * z1 - z2 * z2).max(0.0);
            let k = (q - 2) as f64;
            total += -0.5 * k * (log_2pi + sp.ln()) - perp_ss / (2.0 * sp);
        }
    }
    total
}

fn identity_flat(q: usize) -> Vec<f64> {
    let mut v = vec![0.0f64; q * q];
    for i in 0..q {
        v[i * q + i] = 1.0;
    }
    v
}

fn logsumexp(xs: &[f64]) -> f64 {
    let m = xs.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
    if !m.is_finite() {
        return m;
    }
    let s: f64 = xs.iter().map(|&x| (x - m).exp()).sum();
    m + s.ln()
}

/// Overflow-safe `ln I0(x)` for `x ≥ 0` (Abramowitz & Stegun 9.8.1 / 9.8.2). The
/// large branch factors the `e^x/√x` envelope into the log so it never
/// materialises `e^x`.
fn ln_i0(x: f64) -> f64 {
    let ax = x.abs();
    if ax < 3.75 {
        let t = ax / 3.75;
        let t2 = t * t;
        (1.0 + t2
            * (3.5156229
                + t2 * (3.0899424
                    + t2 * (1.2067492 + t2 * (0.2659732 + t2 * (0.0360768 + t2 * 0.0045813))))))
            .ln()
    } else {
        let y = 3.75 / ax;
        let poly = 0.39894228
            + y * (0.01328592
                + y * (0.00225319
                    + y * (-0.00157565
                        + y * (0.00916281
                            + y * (-0.02057706
                                + y * (0.02635537 + y * (-0.01647633 + y * 0.00392377)))))));
        ax - 0.5 * ax.ln() + poly.ln()
    }
}

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

    // Deterministic, allocation-free xorshift PRNG so the null-battery FDR
    // simulation is reproducible without a rand dependency.
    struct Rng(u64);
    impl Rng {
        fn u01(&mut self) -> f64 {
            let mut x = self.0;
            x ^= x << 13;
            x ^= x >> 7;
            x ^= x << 17;
            self.0 = x;
            // 53-bit mantissa in [0,1).
            ((x >> 11) as f64) / ((1u64 << 53) as f64)
        }
        fn normal(&mut self) -> f64 {
            // Box–Muller.
            let u1 = self.u01().max(1.0e-300);
            let u2 = self.u01();
            (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos()
        }
    }

    /// Draw `n` rows from the rank-1 LINEAR SHELL null: `y = t·v + noise`, a
    /// 1-D affine subspace with isotropic Gaussian scatter. This is the null the
    /// acceptance criterion requires the battery to simulate.
    fn draw_line(rng: &mut Rng, n: usize, q: usize, spread: f64, noise: f64) -> Array2<f64> {
        let axis: Vec<f64> = {
            let mut a: Vec<f64> = (0..q).map(|_| rng.normal()).collect();
            let nrm = (a.iter().map(|x| x * x).sum::<f64>()).sqrt().max(1.0e-12);
            for x in &mut a {
                *x /= nrm;
            }
            a
        };
        let mut out = Array2::<f64>::zeros((n, q));
        for i in 0..n {
            let t = spread * rng.normal();
            for j in 0..q {
                out[[i, j]] = t * axis[j] + noise * rng.normal();
            }
        }
        out
    }

    /// Draw `n` rows from a genuine CURVED RING in the first two coordinates:
    /// `(R cosθ, R sinθ) + noise`, θ uniform — the alternative is true here.
    fn draw_ring(rng: &mut Rng, n: usize, q: usize, radius: f64, noise: f64) -> Array2<f64> {
        let mut out = Array2::<f64>::zeros((n, q));
        for i in 0..n {
            let theta = 2.0 * std::f64::consts::PI * rng.u01();
            out[[i, 0]] = radius * theta.cos() + noise * rng.normal();
            out[[i, 1]] = radius * theta.sin() + noise * rng.normal();
            for j in 2..q {
                out[[i, j]] = noise * rng.normal();
            }
        }
        out
    }

    #[test]
    fn evalue_null_expectation_leq_one() {
        // Monte-Carlo estimate of E_{H0}[E] for the concrete construction under
        // the linear-shell null. The UI guarantee is E_{H0}[E] ≤ 1; we allow a
        // Monte-Carlo margin.
        let mut rng = Rng(0x1234_5678_9abc_def0);
        let trials = 400;
        let mut sum_e = 0.0;
        for _ in 0..trials {
            let coords = draw_line(&mut rng, 80, 3, 2.0, 0.3);
            let log_e = shell_vs_ring_log_evalue(&coords, 2, 1.0e-6).unwrap();
            sum_e += log_e.exp();
        }
        let mean_e = sum_e / trials as f64;
        assert!(
            mean_e <= 1.25,
            "UI e-value must satisfy E_H0[E] <= 1 (MC margin); got mean E = {mean_e}"
        );
    }

    #[test]
    fn null_battery_rank1_shell_fdr_controlled() {
        // Acceptance criterion: simulate the rank-1 shell null, run the full
        // pipeline (family of candidates -> e-values -> e-BH), and verify the
        // empirical FDR <= alpha. Under the GLOBAL null every candidate is
        // false, so any rejection is a false discovery and FDR = P(reject >= 1).
        let alpha = 0.1;
        let mut rng = Rng(0xdead_beef_0bad_f00d);
        let sims = 300;
        let family_size = 12;
        let mut false_discovery_sims = 0;
        for _ in 0..sims {
            let log_e: Vec<f64> = (0..family_size)
                .map(|_| {
                    let coords = draw_line(&mut rng, 60, 3, 2.0, 0.3);
                    shell_vs_ring_log_evalue(&coords, 2, 1.0e-6).unwrap()
                })
                .collect();
            let cert = family_fdr_certificate(log_e, alpha);
            if !cert.rejected.is_empty() {
                false_discovery_sims += 1;
            }
        }
        let empirical_fdr = false_discovery_sims as f64 / sims as f64;
        assert!(
            empirical_fdr <= alpha + 0.03,
            "empirical FDR {empirical_fdr} must be <= alpha {alpha} (MC margin)"
        );
    }

    #[test]
    fn power_under_curved_signal() {
        // A family where one candidate is a genuine ring among linear-shell
        // decoys: e-BH should discover the ring (power > 0) while the guarantee
        // still holds. This checks the construction is not trivially inert.
        let alpha = 0.1;
        let mut rng = Rng(0x0123_4567_89ab_cdef);
        let sims = 60;
        let mut discovered = 0;
        for _ in 0..sims {
            let mut log_e = Vec::new();
            // Index 0 is the true ring; the rest are linear-shell nulls.
            let ring = draw_ring(&mut rng, 240, 3, 3.0, 0.25);
            log_e.push(shell_vs_ring_log_evalue(&ring, 2, 1.0e-6).unwrap());
            for _ in 0..6 {
                let line = draw_line(&mut rng, 240, 3, 2.0, 0.3);
                log_e.push(shell_vs_ring_log_evalue(&line, 2, 1.0e-6).unwrap());
            }
            let cert = family_fdr_certificate(log_e, alpha);
            if cert.rejected.contains(&0) {
                discovered += 1;
            }
        }
        let power = discovered as f64 / sims as f64;
        assert!(
            power > 0.5,
            "the genuine ring should be discovered with real power; got {power}"
        );
    }

    #[test]
    fn degenerate_candidate_never_rejected() {
        // A single-channel candidate cannot support a ring and must bank -inf,
        // never contributing a false discovery.
        let coords =
            Array2::<f64>::from_shape_vec((10, 1), (0..10).map(|x| x as f64).collect()).unwrap();
        let log_e = shell_vs_ring_log_evalue(&coords, 2, 1.0e-6).unwrap();
        assert!(
            log_e.is_infinite() && log_e < 0.0,
            "q<2 must give -inf log_e"
        );
        let cert = family_fdr_certificate(vec![log_e], 0.1);
        assert!(cert.rejected.is_empty());
    }
}