gam-sae 0.3.149

//! Standing null battery and spike-in calibration for topology/nerve claims.
//!
//! A Betti, nerve, or conditionality statistic is only interpretable together
//! with negative controls that preserve mundane structure while destroying the
//! claimed one. This module supplies the reusable harness: generate standing
//! nulls, run an arbitrary scalar audit on observed and null data, and attach a
//! spike-in power curve for a manufactured circle inside real residual noise.

use ndarray::{Array1, Array2, ArrayView1, ArrayView2};
use rand::RngExt;
use rand::SeedableRng;
use rand::rngs::StdRng;
use rand::seq::SliceRandom;
use std::f64::consts::PI;

/// Direction of the claim statistic under the null.
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
pub enum Tail {
    Larger,
    Smaller,
}

impl Tail {
    pub fn as_str(self) -> &'static str {
        match self {
            Self::Larger => "larger",
            Self::Smaller => "smaller",
        }
    }
}

/// Plus-one-corrected Monte Carlo tail probability and its sampling uncertainty.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct EmpiricalPValue {
    pub p_value: f64,
    pub monte_carlo_standard_error: f64,
    pub extreme_draws: usize,
    pub draws: usize,
}

/// Evaluate a scalar statistic against an explicit Monte Carlo null distribution.
///
/// This is the single plus-one correction used by every native null-calibrated
/// report. Ties count as extreme in the requested tail, so a null statistic equal
/// to the observation is evidence against rejection rather than a free win.
pub fn empirical_p_value(
    observed: f64,
    null_samples: &[f64],
    tail: Tail,
) -> Result<EmpiricalPValue, String> {
    require_finite(observed, "observed statistic")?;
    if null_samples.is_empty() {
        return Err("empirical p-value requires at least one null draw".to_string());
    }
    for (draw, &sample) in null_samples.iter().enumerate() {
        require_finite(sample, &format!("null statistic draw {draw}"))?;
    }
    let extreme_draws = null_samples
        .iter()
        .filter(|&&sample| match tail {
            Tail::Larger => sample >= observed,
            Tail::Smaller => sample <= observed,
        })
        .count();
    let draws = null_samples.len();
    let p_value = (extreme_draws as f64 + 1.0) / (draws as f64 + 1.0);
    let monte_carlo_standard_error = (p_value * (1.0 - p_value) / (draws as f64 + 1.0)).sqrt();
    Ok(EmpiricalPValue {
        p_value,
        monte_carlo_standard_error,
        extreme_draws,
        draws,
    })
}

/// One member of the standing null battery.
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
pub enum NullKind {
    /// Randomize each column's Fourier phases along token order, preserving its
    /// one-dimensional power spectrum while destroying coherent ordered shape.
    PhaseRandomized,
    /// Apply a Haar-like random orthogonal basis change. Coordinate-dependent
    /// claims should fail; basis-invariant claims should be unchanged.
    RandomRotation,
    /// Permute token rows, preserving the activation cloud and column marginals
    /// while destroying row-order structure.
    TokenShuffle,
    /// Permute every column independently. This preserves each one-dimensional
    /// empirical marginal exactly while destroying cross-column geometry such
    /// as an unordered circle or ring of clusters.
    PerDimensionShuffle,
    /// Draw independent Gaussian principal-component scores with the same mean
    /// and variance in every supplied PC column, preserving the eigenspectrum
    /// while destroying cyclic/non-Gaussian structure.
    MatchedSpectrumGaussian,
    /// Draw a multivariate Gaussian with the observed mean and full covariance.
    /// Unlike `MatchedSpectrumGaussian`, this is expressed in the supplied
    /// coordinates and preserves their cross-covariance as well as their scales.
    CovarianceMatchedGaussian,
    /// Use separately supplied random-weight activations from the same
    /// architecture, moment-matched to the observed matrix shape.
    ArchitectureMatchedRandomWeight,
}

impl NullKind {
    pub fn as_str(self) -> &'static str {
        match self {
            Self::PhaseRandomized => "phase_randomized",
            Self::RandomRotation => "random_rotation",
            Self::TokenShuffle => "token_shuffle",
            Self::PerDimensionShuffle => "per_dimension_shuffle",
            Self::MatchedSpectrumGaussian => "matched_spectrum_gaussian",
            Self::CovarianceMatchedGaussian => "covariance_matched_gaussian",
            Self::ArchitectureMatchedRandomWeight => "architecture_matched_random_weight",
        }
    }
}

/// Configuration for [`run_null_battery`].
#[derive(Clone, Debug)]
pub struct NullBatteryConfig {
    pub replicates: usize,
    pub seed: u64,
    pub kinds: Vec<NullKind>,
    pub tail: Tail,
}

impl NullBatteryConfig {
    pub fn standing(replicates: usize, seed: u64) -> Self {
        Self {
            replicates,
            seed,
            kinds: vec![
                NullKind::PhaseRandomized,
                NullKind::RandomRotation,
                NullKind::TokenShuffle,
                NullKind::PerDimensionShuffle,
                NullKind::MatchedSpectrumGaussian,
                NullKind::CovarianceMatchedGaussian,
                NullKind::ArchitectureMatchedRandomWeight,
            ],
            tail: Tail::Larger,
        }
    }
}

/// Null distribution summary for a single null kind.
#[derive(Clone, Debug, PartialEq)]
pub struct NullSummary {
    pub kind: NullKind,
    pub tail: Tail,
    pub observed: f64,
    pub n: usize,
    pub mean: f64,
    pub sd: f64,
    pub min: f64,
    pub q25: f64,
    pub median: f64,
    pub q75: f64,
    pub max: f64,
    pub z: f64,
    pub p_value: f64,
    pub monte_carlo_standard_error: f64,
    pub extreme_draws: usize,
    /// Statistic values in draw order. Keeping seed order (rather than sorting
    /// this ledger for quantiles) makes a persisted null artifact replayable.
    pub samples: Vec<f64>,
}

/// Complete null-battery result for one scalar audit statistic.
#[derive(Clone, Debug)]
pub struct NullBatteryReport {
    pub observed: f64,
    pub summaries: Vec<NullSummary>,
}

/// Synthetic topology planted into residual activations for spike-in calibration.
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
pub enum SpikeInShape {
    Circle,
    Torus,
}

impl SpikeInShape {
    fn signal_rank(self) -> usize {
        match self {
            SpikeInShape::Circle => 2,
            SpikeInShape::Torus => 4,
        }
    }

    fn expected_betti(self) -> (usize, usize, usize) {
        match self {
            SpikeInShape::Circle => (1, 1, 0),
            SpikeInShape::Torus => (1, 2, 1),
        }
    }
}

/// How trial noise is drawn before a synthetic topology is injected.
#[derive(Clone, Copy, Debug, PartialEq)]
pub enum SpikeInNoiseMode {
    /// Resample centered residual rows with replacement. This keeps the actual
    /// post-sink-peel residual covariance and heavy-tailed row distribution in
    /// the calibration loop.
    EmpiricalResidualBootstrap,
    /// Preserve each residual channel's one-dimensional spectrum while
    /// destroying coherent topology.
    PhaseRandomizedResidual,
    /// Treat the supplied matrix as a moment donor and draw a covariance-
    /// matched heavy-tailed surrogate with the donor's average marginal excess
    /// kurtosis. Use this when real residual activations are unavailable.
    CovarianceHeavyTailSurrogate,
}

/// Configuration for [`spike_in_roc_curve`].
#[derive(Clone, Debug)]
pub struct SpikeInRocConfig {
    pub shape: SpikeInShape,
    pub snrs: Vec<f64>,
    pub trials: usize,
    pub seed: u64,
    pub fpr_levels: Vec<f64>,
    pub noise_mode: SpikeInNoiseMode,
}

impl SpikeInRocConfig {
    pub fn circle(snrs: Vec<f64>, trials: usize, seed: u64) -> Self {
        Self {
            shape: SpikeInShape::Circle,
            snrs,
            trials,
            seed,
            // reporting-only: FPR operating points at which spike-in ROC TPR is reported.
            fpr_levels: vec![0.01, 0.05, 0.10],
            noise_mode: SpikeInNoiseMode::PhaseRandomizedResidual,
        }
    }

    pub fn torus(snrs: Vec<f64>, trials: usize, seed: u64) -> Self {
        Self {
            shape: SpikeInShape::Torus,
            snrs,
            trials,
            seed,
            // reporting-only: FPR operating points at which spike-in ROC TPR is reported.
            fpr_levels: vec![0.01, 0.05, 0.10],
            noise_mode: SpikeInNoiseMode::EmpiricalResidualBootstrap,
        }
    }
}

/// Explicit covariance/heavy-tail moments for generating a surrogate residual
/// activation matrix when real post-sink-peel activations are unavailable.
#[derive(Clone, Debug)]
pub struct ResidualMomentSpec {
    pub rows: usize,
    pub mean: Array1<f64>,
    pub covariance: Array2<f64>,
    pub excess_kurtosis: f64,
}

/// Block-chart promotion plus topology-audit verdict for one detection run.
#[derive(Clone, Debug)]
pub struct DetectionPipelineReport {
    pub shape: SpikeInShape,
    pub statistic: f64,
    pub promoted: bool,
    pub topology: TopologyAuditReport,
    pub detected: bool,
}

/// Lightweight topology audit payload carried by spike-in ROC points. The
/// measured Betti numbers are the detector's auditable claim for the supplied
/// residual matrix, not a latched topology label.
#[derive(Clone, Debug)]
pub struct TopologyAuditReport {
    pub expected_betti0: usize,
    pub expected_betti1: usize,
    pub expected_betti2: usize,
    pub measured_betti0: usize,
    pub measured_betti1: usize,
    pub measured_betti2: usize,
    pub rank_energy: f64,
    pub spectral_balance: f64,
    pub residual_tail_energy: f64,
    pub accepted: bool,
}

/// One ROC operating point at a fixed false-positive-rate threshold.
#[derive(Clone, Debug)]
pub struct SpikeInRocThreshold {
    pub false_positive_rate: f64,
    pub threshold: f64,
    pub true_positive_rate: f64,
}

/// ROC payload for one injected SNR.
#[derive(Clone, Debug)]
pub struct SpikeInRocPoint {
    pub snr: f64,
    pub trials: usize,
    pub mean_stat: f64,
    pub promoted_fraction: f64,
    pub topology_accept_fraction: f64,
    pub roc: Vec<SpikeInRocThreshold>,
}

/// Claim payload carrying the full spike-in ROC calibration. This is the shape
/// topology, Betti, and conditionality claims should ship when the detector has
/// an explicit operating FPR.
#[derive(Clone, Debug)]
pub struct CalibratedRocClaimReport {
    pub claim: String,
    pub claimed_snr: f64,
    pub claimed_false_positive_rate: f64,
    pub nulls: NullBatteryReport,
    pub spike_in_roc: Vec<SpikeInRocPoint>,
    pub claimed_snr_power: f64,
}

/// Compact null-calibrated audit payload attached directly to emitted claims.
#[derive(Clone, Debug)]
pub struct ClaimNullCalibration {
    pub claim: String,
    pub observed_statistic: f64,
    pub null_pvalue: f64,
    pub null_z: f64,
    pub claimed_snr: f64,
    pub claimed_false_positive_rate: f64,
    pub spikein_power: f64,
    pub null_distribution: Vec<NullSummary>,
    pub spike_in_roc: Vec<SpikeInRocPoint>,
}

impl ClaimNullCalibration {
    pub fn from_calibrated_roc(report: CalibratedRocClaimReport) -> Self {
        let null_pvalue = primary_null_pvalue(&report.nulls);
        let null_z = primary_null_z(&report.nulls);
        Self {
            claim: report.claim,
            observed_statistic: report.nulls.observed,
            null_pvalue,
            null_z,
            claimed_snr: report.claimed_snr,
            claimed_false_positive_rate: report.claimed_false_positive_rate,
            spikein_power: report.claimed_snr_power,
            null_distribution: report.nulls.summaries,
            spike_in_roc: report.spike_in_roc,
        }
    }
}

pub fn primary_null_pvalue(report: &NullBatteryReport) -> f64 {
    let mut selected = Vec::new();
    for summary in &report.summaries {
        if matches!(
            summary.kind,
            NullKind::PhaseRandomized
                | NullKind::MatchedSpectrumGaussian
                | NullKind::CovarianceMatchedGaussian
                | NullKind::ArchitectureMatchedRandomWeight
        ) {
            selected.push(summary.p_value);
        }
    }
    if selected.is_empty() {
        selected.extend(report.summaries.iter().map(|summary| summary.p_value));
    }
    selected.into_iter().fold(0.0_f64, f64::max)
}

pub fn primary_null_z(report: &NullBatteryReport) -> f64 {
    let mut selected = Vec::new();
    for summary in &report.summaries {
        if matches!(
            summary.kind,
            NullKind::PhaseRandomized
                | NullKind::MatchedSpectrumGaussian
                | NullKind::CovarianceMatchedGaussian
                | NullKind::ArchitectureMatchedRandomWeight
        ) {
            selected.push(summary.z);
        }
    }
    if selected.is_empty() {
        selected.extend(report.summaries.iter().map(|summary| summary.z));
    }
    selected.into_iter().fold(f64::INFINITY, f64::min)
}

/// Run an arbitrary scalar audit on the observed matrix and each requested null.
///
/// If [`NullKind::ArchitectureMatchedRandomWeight`] is requested,
/// `random_weight` must be a real activation matrix harvested from a random-
/// weight model with the same architecture family. The harness only resamples
/// and moment-matches it; it does not fabricate an architecture null.
pub fn run_null_battery<F>(
    data: ArrayView2<'_, f64>,
    random_weight: Option<ArrayView2<'_, f64>>,
    config: &NullBatteryConfig,
    mut audit: F,
) -> Result<NullBatteryReport, String>
where
    F: FnMut(ArrayView2<'_, f64>) -> Result<f64, String>,
{
    validate_matrix(data, "observed")?;
    if config.replicates == 0 {
        return Err("null battery requires at least one replicate".to_string());
    }
    let observed = audit(data)?;
    require_finite(observed, "observed statistic")?;
    let mut summaries = Vec::with_capacity(config.kinds.len());

    for &kind in &config.kinds {
        let mut samples = Vec::with_capacity(config.replicates);
        for rep in 0..config.replicates {
            let rep_seed = mix_seed(config.seed, kind_seed(kind), rep as u64);
            let null = match kind {
                NullKind::PhaseRandomized => phase_randomized_surrogate(data, rep_seed)?,
                NullKind::RandomRotation => random_rotation_null(data, rep_seed)?,
                NullKind::TokenShuffle => token_shuffle_null(data, rep_seed)?,
                NullKind::PerDimensionShuffle => per_dimension_shuffle_null(data, rep_seed)?,
                NullKind::MatchedSpectrumGaussian => {
                    matched_spectrum_gaussian_null(data, rep_seed)?
                }
                NullKind::CovarianceMatchedGaussian => {
                    covariance_matched_gaussian_null(data, rep_seed)?
                }
                NullKind::ArchitectureMatchedRandomWeight => {
                    let Some(rw) = random_weight else {
                        return Err(
                            "architecture-matched random-weight null requested without random_weight activations"
                                .to_string(),
                        );
                    };
                    architecture_matched_random_weight_null(data, rw, rep_seed)?
                }
            };
            let stat = audit(null.view())?;
            require_finite(stat, "null statistic")?;
            samples.push(stat);
        }
        summaries.push(summarize_null_distribution(
            kind,
            observed,
            samples,
            config.tail,
        )?);
    }

    Ok(NullBatteryReport {
        observed,
        summaries,
    })
}

/// Fourier phase-randomized surrogate, independently per activation channel.
pub fn phase_randomized_surrogate(
    data: ArrayView2<'_, f64>,
    seed: u64,
) -> Result<Array2<f64>, String> {
    validate_matrix(data, "phase-randomized input")?;
    let n = data.nrows();
    let p = data.ncols();
    let mut out = Array2::<f64>::zeros((n, p));
    let mut rng = StdRng::seed_from_u64(seed);
    for col in 0..p {
        let mut re = vec![0.0_f64; n];
        let mut im = vec![0.0_f64; n];
        for k in 0..n {
            for t in 0..n {
                let angle = -2.0 * PI * (k as f64) * (t as f64) / (n as f64);
                let x = data[[t, col]];
                re[k] += x * angle.cos();
                im[k] += x * angle.sin();
            }
        }
        let last_paired = (n.saturating_sub(1)) / 2;
        for k in 1..=last_paired {
            let amp = (re[k] * re[k] + im[k] * im[k]).sqrt();
            let phase = rng.random_range(0.0..(2.0 * PI));
            re[k] = amp * phase.cos();
            im[k] = amp * phase.sin();
            let mirror = n - k;
            re[mirror] = re[k];
            im[mirror] = -im[k];
        }
        if n % 2 == 0 {
            let nyquist = n / 2;
            im[nyquist] = 0.0;
        }
        for t in 0..n {
            let mut acc = 0.0_f64;
            for k in 0..n {
                let angle = 2.0 * PI * (k as f64) * (t as f64) / (n as f64);
                acc += re[k] * angle.cos() - im[k] * angle.sin();
            }
            out[[t, col]] = acc / (n as f64);
        }
    }
    Ok(out)
}

/// Multiply the activation matrix by a seeded random orthogonal matrix.
pub fn random_rotation_null(data: ArrayView2<'_, f64>, seed: u64) -> Result<Array2<f64>, String> {
    validate_matrix(data, "random-rotation input")?;
    let p = data.ncols();
    let q = random_orthogonal(p, seed)?;
    Ok(data.dot(&q))
}

/// Deterministically seeded token-row shuffle.
pub fn token_shuffle_null(data: ArrayView2<'_, f64>, seed: u64) -> Result<Array2<f64>, String> {
    validate_matrix(data, "token-shuffle input")?;
    let mut order: Vec<usize> = (0..data.nrows()).collect();
    let mut rng = StdRng::seed_from_u64(seed);
    order.shuffle(&mut rng);
    let mut out = Array2::<f64>::zeros(data.raw_dim());
    for (dst, &src) in order.iter().enumerate() {
        out.row_mut(dst).assign(&data.row(src));
    }
    Ok(out)
}

/// Independently permute each dimension of an activation/coordinate matrix.
///
/// Every output column is a bit-for-bit permutation of the corresponding input
/// column, so all empirical marginal moments and tails are preserved exactly.
/// The independently seeded permutations break cross-dimensional geometry; in
/// particular, this is the matched shuffle control required for an unordered
/// circle census, whereas [`token_shuffle_null`] leaves that point cloud intact.
pub fn per_dimension_shuffle_null(
    data: ArrayView2<'_, f64>,
    seed: u64,
) -> Result<Array2<f64>, String> {
    validate_matrix(data, "per-dimension-shuffle input")?;
    let mut rng = StdRng::seed_from_u64(seed);
    let mut out = Array2::<f64>::zeros(data.raw_dim());
    for col in 0..data.ncols() {
        let mut order: Vec<usize> = (0..data.nrows()).collect();
        order.shuffle(&mut rng);
        for (dst, &src) in order.iter().enumerate() {
            out[[dst, col]] = data[[src, col]];
        }
    }
    Ok(out)
}

/// Gaussian matched-spectrum null for a matrix already expressed in principal-
/// component coordinates.
///
/// Each output PC column is an independent normal draw with the observed
/// column's mean and sample standard deviation. Consequently the per-PC
/// eigenspectrum is preserved while cyclic ordering, higher moments, and
/// cross-row manifold structure are destroyed. The explicit PC-coordinate
/// contract avoids a hidden eigen-decomposition or a Python-side reimplementation.
pub fn matched_spectrum_gaussian_null(
    pc_scores: ArrayView2<'_, f64>,
    seed: u64,
) -> Result<Array2<f64>, String> {
    validate_matrix(pc_scores, "matched-spectrum PC scores")?;
    let mean = column_mean(pc_scores);
    let sd = column_sd(pc_scores, mean.view());
    let mut rng = StdRng::seed_from_u64(seed);
    let mut out = Array2::<f64>::zeros(pc_scores.raw_dim());
    for row in 0..out.nrows() {
        for pc in 0..out.ncols() {
            out[[row, pc]] = mean[pc] + sd[pc] * standard_normal(&mut rng);
        }
    }
    Ok(out)
}

/// Gaussian null matched to the observed mean and full covariance.
///
/// The covariance square root is formed from the symmetric eigendecomposition,
/// not by adding a ridge: zero empirical eigenvalues stay exactly zero, so a
/// rank-deficient coordinate cloud produces a rank-deficient matched null
/// rather than a different full-rank distribution. Tiny negative eigenvalues
/// produced solely by floating-point eigensolver roundoff are projected to the
/// positive-semidefinite cone at zero.
pub fn covariance_matched_gaussian_null(
    data: ArrayView2<'_, f64>,
    seed: u64,
) -> Result<Array2<f64>, String> {
    use faer::Side;
    use gam_linalg::faer_ndarray::FaerEigh;

    validate_matrix(data, "covariance-matched Gaussian input")?;
    let n = data.nrows();
    let p = data.ncols();
    let mean = column_mean(data);
    let mut covariance = Array2::<f64>::zeros((p, p));
    for row in data.rows() {
        for a in 0..p {
            let da = row[a] - mean[a];
            for b in 0..=a {
                covariance[[a, b]] += da * (row[b] - mean[b]);
            }
        }
    }
    for a in 0..p {
        for b in 0..=a {
            let value = covariance[[a, b]] / n as f64;
            covariance[[a, b]] = value;
            covariance[[b, a]] = value;
        }
    }
    let (eigenvalues, eigenvectors) = covariance.eigh(Side::Lower).map_err(|error| {
        format!("covariance-matched Gaussian eigendecomposition failed: {error}")
    })?;

    let mut rng = StdRng::seed_from_u64(seed);
    let mut out = Array2::<f64>::zeros((n, p));
    let mut scaled_normal = vec![0.0_f64; p];
    for row in 0..n {
        for axis in 0..p {
            scaled_normal[axis] = eigenvalues[axis].max(0.0).sqrt() * standard_normal(&mut rng);
        }
        for col in 0..p {
            let mut draw = mean[col];
            for axis in 0..p {
                draw += eigenvectors[[col, axis]] * scaled_normal[axis];
            }
            out[[row, col]] = draw;
        }
    }
    Ok(out)
}

/// Resample real random-weight activations and match observed column moments.
pub fn architecture_matched_random_weight_null(
    observed: ArrayView2<'_, f64>,
    random_weight: ArrayView2<'_, f64>,
    seed: u64,
) -> Result<Array2<f64>, String> {
    validate_matrix(observed, "observed")?;
    validate_matrix(random_weight, "random_weight")?;
    if random_weight.ncols() != observed.ncols() {
        return Err(format!(
            "random_weight ncols {} != observed ncols {}",
            random_weight.ncols(),
            observed.ncols()
        ));
    }
    let n = observed.nrows();
    let p = observed.ncols();
    let mut rng = StdRng::seed_from_u64(seed);
    let obs_mean = column_mean(observed);
    let obs_sd = column_sd(observed, obs_mean.view());
    let rw_mean = column_mean(random_weight);
    let rw_sd = column_sd(random_weight, rw_mean.view());
    let mut out = Array2::<f64>::zeros((n, p));
    for i in 0..n {
        let src = rng.random_range(0..random_weight.nrows());
        for j in 0..p {
            let centered = random_weight[[src, j]] - rw_mean[j];
            let scaled = if rw_sd[j] > 0.0 {
                centered * (obs_sd[j] / rw_sd[j])
            } else {
                0.0
            };
            out[[i, j]] = obs_mean[j] + scaled;
        }
    }
    Ok(out)
}

/// Bootstrap a realistic residual-noise matrix from post-sink-peel activations.
///
/// Rows are sampled with replacement after centering, so the calibration sees the
/// empirical residual covariance, marginal tails, and row-level outliers rather
/// than an idealized Gaussian null.
pub fn empirical_residual_bootstrap(
    residuals: ArrayView2<'_, f64>,
    seed: u64,
) -> Result<Array2<f64>, String> {
    validate_matrix(residuals, "empirical residual bootstrap input")?;
    let n = residuals.nrows();
    let p = residuals.ncols();
    let mean = column_mean(residuals);
    let mut rng = StdRng::seed_from_u64(seed);
    let mut out = Array2::<f64>::zeros((n, p));
    for i in 0..n {
        let src = rng.random_range(0..n);
        for j in 0..p {
            out[[i, j]] = residuals[[src, j]] - mean[j];
        }
    }
    Ok(out)
}

/// Build a covariance- and kurtosis-matched residual surrogate from explicit
/// moments when real residual activations cannot be supplied to the harness.
pub fn residual_surrogate_from_moments(
    spec: &ResidualMomentSpec,
    seed: u64,
) -> Result<Array2<f64>, String> {
    if spec.rows == 0 {
        return Err("residual surrogate requires at least one row".to_string());
    }
    if spec.mean.is_empty() {
        return Err("residual surrogate mean must be non-empty".to_string());
    }
    if spec.covariance.nrows() != spec.mean.len() || spec.covariance.ncols() != spec.mean.len() {
        return Err(format!(
            "residual surrogate covariance shape {}x{} does not match mean length {}",
            spec.covariance.nrows(),
            spec.covariance.ncols(),
            spec.mean.len()
        ));
    }
    if !spec.excess_kurtosis.is_finite() || spec.excess_kurtosis < 0.0 {
        return Err(format!(
            "residual surrogate excess_kurtosis must be finite and non-negative, got {}",
            spec.excess_kurtosis
        ));
    }
    for (idx, &value) in spec.mean.iter().enumerate() {
        if !value.is_finite() {
            return Err(format!(
                "residual surrogate mean[{idx}] is not finite: {value}"
            ));
        }
    }
    validate_matrix(spec.covariance.view(), "residual surrogate covariance")?;

    let chol = cholesky_lower(spec.covariance.view())?;
    let p = spec.mean.len();
    let mut rng = StdRng::seed_from_u64(seed);
    let mut out = Array2::<f64>::zeros((spec.rows, p));
    for i in 0..spec.rows {
        let radial = standardized_lognormal_radial(spec.excess_kurtosis, &mut rng);
        let mut z = Array1::<f64>::zeros(p);
        for j in 0..p {
            z[j] = standard_normal(&mut rng);
        }
        let correlated = chol.dot(&z);
        for j in 0..p {
            out[[i, j]] = spec.mean[j] + radial * correlated[j];
        }
    }
    Ok(out)
}

/// Draw a surrogate matrix using moments estimated from a residual donor.
pub fn residual_surrogate_matching(
    residuals: ArrayView2<'_, f64>,
    seed: u64,
) -> Result<Array2<f64>, String> {
    validate_matrix(residuals, "residual surrogate donor")?;
    let spec = residual_moment_spec(residuals)?;
    residual_surrogate_from_moments(&spec, seed)
}

/// Estimate the moment specification used by [`residual_surrogate_from_moments`].
pub fn residual_moment_spec(residuals: ArrayView2<'_, f64>) -> Result<ResidualMomentSpec, String> {
    validate_matrix(residuals, "residual moment donor")?;
    let mean = column_mean(residuals);
    let centered = centered_matrix(residuals);
    let n = centered.nrows();
    let denom = if n > 1 { (n - 1) as f64 } else { 1.0 };
    let covariance = centered.t().dot(&centered) / denom;
    let mut kurtosis_sum = 0.0_f64;
    let mut kurtosis_count = 0usize;
    for j in 0..centered.ncols() {
        let mut second = 0.0_f64;
        let mut fourth = 0.0_f64;
        for i in 0..centered.nrows() {
            let v = centered[[i, j]];
            let v2 = v * v;
            second += v2;
            fourth += v2 * v2;
        }
        second /= n as f64;
        fourth /= n as f64;
        if second > f64::MIN_POSITIVE {
            kurtosis_sum += (fourth / (second * second) - 3.0).max(0.0);
            kurtosis_count += 1;
        }
    }
    let excess_kurtosis = if kurtosis_count > 0 {
        kurtosis_sum / kurtosis_count as f64
    } else {
        0.0
    };
    Ok(ResidualMomentSpec {
        rows: residuals.nrows(),
        mean,
        covariance,
        excess_kurtosis,
    })
}

/// Inject a circle into a random two-plane of the supplied residual/noise matrix.
///
/// `snr` is signal RMS divided by the input matrix RMS. The returned matrix keeps
/// the original residuals and adds the synthetic ground-truth circle.
pub fn inject_circle_spike(
    noise: ArrayView2<'_, f64>,
    snr: f64,
    seed: u64,
) -> Result<Array2<f64>, String> {
    validate_matrix(noise, "spike-in noise")?;
    if !snr.is_finite() || snr < 0.0 {
        return Err(format!("snr must be finite and non-negative, got {snr}"));
    }
    let n = noise.nrows();
    let p = noise.ncols();
    if p < 2 {
        return Err("circle spike-in requires at least two columns".to_string());
    }
    let mut rng = StdRng::seed_from_u64(seed);
    let basis = random_orthogonal(p, seed ^ 0xA53A_9E1D_2C6B_7F40)?;
    let rms = matrix_rms(noise);
    let amp = snr * rms * (p as f64).sqrt();
    let phase0 = rng.random_range(0.0..(2.0 * PI));
    let mut out = noise.to_owned();
    for i in 0..n {
        let theta = phase0 + 2.0 * PI * (i as f64) / (n as f64);
        for j in 0..p {
            out[[i, j]] += amp * (theta.cos() * basis[[j, 0]] + theta.sin() * basis[[j, 1]]);
        }
    }
    Ok(out)
}

/// Inject a product torus into a random four-plane of the supplied residual
/// matrix at controlled RMS SNR.
pub fn inject_torus_spike(
    noise: ArrayView2<'_, f64>,
    snr: f64,
    seed: u64,
) -> Result<Array2<f64>, String> {
    validate_matrix(noise, "torus spike-in noise")?;
    if !snr.is_finite() || snr < 0.0 {
        return Err(format!("snr must be finite and non-negative, got {snr}"));
    }
    let n = noise.nrows();
    let p = noise.ncols();
    if p < 4 {
        return Err("torus spike-in requires at least four columns".to_string());
    }
    let mut rng = StdRng::seed_from_u64(seed);
    let basis = random_orthogonal(p, seed ^ 0x7055_1A7E_4D2C_A11B)?;
    let rms = matrix_rms(noise);
    let amp = snr * rms * ((p as f64) / 2.0).sqrt();
    let phase0 = rng.random_range(0.0..(2.0 * PI));
    let phase1 = rng.random_range(0.0..(2.0 * PI));
    let winding = coprime_winding(n);
    let mut out = noise.to_owned();
    for i in 0..n {
        let theta = phase0 + 2.0 * PI * (i as f64) / (n as f64);
        let phi = phase1 + 2.0 * PI * ((winding * i) % n) as f64 / (n as f64);
        for j in 0..p {
            out[[i, j]] += amp
                * (theta.cos() * basis[[j, 0]]
                    + theta.sin() * basis[[j, 1]]
                    + phi.cos() * basis[[j, 2]]
                    + phi.sin() * basis[[j, 3]]);
        }
    }
    Ok(out)
}

/// Inject the configured synthetic topology at controlled RMS SNR.
pub fn inject_spike(
    noise: ArrayView2<'_, f64>,
    shape: SpikeInShape,
    snr: f64,
    seed: u64,
) -> Result<Array2<f64>, String> {
    match shape {
        SpikeInShape::Circle => inject_circle_spike(noise, snr, seed),
        SpikeInShape::Torus => inject_torus_spike(noise, snr, seed),
    }
}

/// Default detector used by the calibration harness: a block-chart promotion
/// gate followed by a topology audit for the requested planted shape.
pub fn default_spike_in_detection_pipeline(
    data: ArrayView2<'_, f64>,
    shape: SpikeInShape,
) -> Result<DetectionPipelineReport, String> {
    validate_matrix(data, "spike-in detection input")?;
    let rank = shape.signal_rank();
    if data.ncols() < rank {
        return Err(format!(
            "{shape:?} detection requires at least {rank} columns, got {}",
            data.ncols()
        ));
    }
    let eigenvalues = leading_covariance_eigenvalues(data, rank + 1)?;
    let statistic = match shape {
        SpikeInShape::Circle => harmonic_circle_detector_stat(data)?,
        SpikeInShape::Torus => {
            let total = centered_total_energy(data);
            let rank_sum = eigenvalues.iter().take(rank).sum::<f64>().max(0.0);
            if total > 0.0 { rank_sum / total } else { 0.0 }
        }
    };
    // reporting-only: heuristic promotion floors on the shape detector statistic;
    // they gate the reported promote/hold verdict, not any estimated quantity.
    let promotion_floor = match shape {
        SpikeInShape::Circle => 0.10, // reporting-only
        SpikeInShape::Torus => 0.55,  // reporting-only
    };
    let promoted = statistic >= promotion_floor;
    let topology = match shape {
        SpikeInShape::Circle => circle_topology_audit_from_harmonic_stat(statistic, &eigenvalues),
        SpikeInShape::Torus => topology_audit_from_spectrum(shape, statistic, &eigenvalues),
    };
    let detected = promoted && topology.accepted;
    Ok(DetectionPipelineReport {
        shape,
        statistic,
        promoted,
        topology,
        detected,
    })
}

/// Run the spike-in ROC calibration on residual noise with a caller-supplied
/// detector. The detector should execute the same block-chart promotion and
/// topology audit used by the production claim.
pub fn spike_in_roc_curve<F>(
    residual_noise: ArrayView2<'_, f64>,
    config: &SpikeInRocConfig,
    mut detector: F,
) -> Result<Vec<SpikeInRocPoint>, String>
where
    F: FnMut(ArrayView2<'_, f64>) -> Result<DetectionPipelineReport, String>,
{
    validate_spike_in_config(config)?;
    validate_matrix(residual_noise, "spike-in residual noise")?;
    if residual_noise.ncols() < config.shape.signal_rank() {
        return Err(format!(
            "{:?} spike-in requires at least {} residual columns, got {}",
            config.shape,
            config.shape.signal_rank(),
            residual_noise.ncols()
        ));
    }

    let mut null_stats = Vec::with_capacity(config.trials);
    for trial in 0..config.trials {
        let noise = draw_spike_in_noise(
            residual_noise,
            config.noise_mode,
            mix_seed(config.seed, 0xA17D_7E57, trial as u64),
        )?;
        let report = detector(noise.view())?;
        require_finite(report.statistic, "null detection statistic")?;
        null_stats.push(report.statistic);
    }
    null_stats.sort_by(|a, b| a.total_cmp(b));
    let thresholds = roc_thresholds(&null_stats, &config.fpr_levels)?;

    let mut curve = Vec::with_capacity(config.snrs.len());
    for &snr in &config.snrs {
        let mut reports = Vec::with_capacity(config.trials);
        let mut stat_sum = 0.0_f64;
        let mut promoted = 0usize;
        let mut topology_accepted = 0usize;
        for trial in 0..config.trials {
            let trial_seed = mix_seed(config.seed, 0, trial as u64);
            let noise = draw_spike_in_noise(residual_noise, config.noise_mode, trial_seed)?;
            let spiked = inject_spike(noise.view(), config.shape, snr, trial_seed ^ 0x5F1E_51A5)?;
            let report = detector(spiked.view())?;
            require_finite(report.statistic, "spike-in detection statistic")?;
            if report.promoted {
                promoted += 1;
            }
            if report.topology.accepted {
                topology_accepted += 1;
            }
            stat_sum += report.statistic;
            reports.push(report);
        }
        let mut roc = Vec::with_capacity(thresholds.len());
        for &(fpr, threshold) in &thresholds {
            let hits = reports
                .iter()
                .filter(|report| {
                    report.promoted && report.topology.accepted && report.statistic > threshold
                })
                .count();
            roc.push(SpikeInRocThreshold {
                false_positive_rate: fpr,
                threshold,
                true_positive_rate: hits as f64 / config.trials as f64,
            });
        }
        curve.push(SpikeInRocPoint {
            snr,
            trials: config.trials,
            mean_stat: stat_sum / config.trials as f64,
            promoted_fraction: promoted as f64 / config.trials as f64,
            topology_accept_fraction: topology_accepted as f64 / config.trials as f64,
            roc,
        });
    }
    Ok(curve)
}

/// Convenience wrapper for the default block-chart/topology detector.
pub fn default_spike_in_roc_curve(
    residual_noise: ArrayView2<'_, f64>,
    config: &SpikeInRocConfig,
) -> Result<Vec<SpikeInRocPoint>, String> {
    let shape = config.shape;
    spike_in_roc_curve(residual_noise, config, |x| {
        default_spike_in_detection_pipeline(x, shape)
    })
}

/// Build a calibrated claim report from a spike-in ROC curve at the requested
/// false-positive-rate operating point.
pub fn calibrated_roc_claim_report(
    claim: impl Into<String>,
    claimed_snr: f64,
    claimed_false_positive_rate: f64,
    nulls: NullBatteryReport,
    spike_in_roc: Vec<SpikeInRocPoint>,
) -> Result<CalibratedRocClaimReport, String> {
    if !claimed_snr.is_finite() || claimed_snr < 0.0 {
        return Err(format!(
            "claimed_snr must be finite and non-negative, got {claimed_snr}"
        ));
    }
    if !claimed_false_positive_rate.is_finite()
        || claimed_false_positive_rate <= 0.0
        || claimed_false_positive_rate >= 1.0
    {
        return Err(format!(
            "claimed_false_positive_rate must be in (0, 1), got {claimed_false_positive_rate}"
        ));
    }

    let mut best_point = None;
    let mut best_distance = f64::INFINITY;
    for point in &spike_in_roc {
        let distance = (point.snr - claimed_snr).abs();
        if distance < best_distance {
            best_distance = distance;
            best_point = Some(point);
        }
    }
    let Some(point) = best_point else {
        return Err("calibrated ROC claim report requires at least one spike-in point".to_string());
    };
    let mut best_power = None;
    let mut best_fpr_distance = f64::INFINITY;
    for threshold in &point.roc {
        let distance = (threshold.false_positive_rate - claimed_false_positive_rate).abs();
        if distance < best_fpr_distance {
            best_fpr_distance = distance;
            best_power = Some(threshold.true_positive_rate);
        }
    }
    let Some(claimed_snr_power) = best_power else {
        return Err("calibrated ROC claim report requires at least one ROC threshold".to_string());
    };
    Ok(CalibratedRocClaimReport {
        claim: claim.into(),
        claimed_snr,
        claimed_false_positive_rate,
        nulls,
        spike_in_roc,
        claimed_snr_power,
    })
}

/// Coordinate-dependent ordered-circle statistic on columns 0 and 1.
pub fn first_two_ordered_circle_stat(data: ArrayView2<'_, f64>) -> Result<f64, String> {
    validate_matrix(data, "first-two circle input")?;
    if data.ncols() < 2 {
        return Err("first-two circle statistic requires at least two columns".to_string());
    }
    let n = data.nrows();
    if n < 4 {
        return Err("first-two circle statistic requires at least four rows".to_string());
    }
    let x = centered_column(data.column(0));
    let y = centered_column(data.column(1));
    let mut pos_re = 0.0_f64;
    let mut pos_im = 0.0_f64;
    let mut neg_re = 0.0_f64;
    let mut neg_im = 0.0_f64;
    let mut energy = 0.0_f64;
    let mut signed_area = 0.0_f64;
    for i in 0..n {
        let theta = 2.0 * PI * i as f64 / n as f64;
        let c = theta.cos();
        let s = theta.sin();
        // Frequency-1 coefficient of z = x + i y. A coherent ordered circle puts
        // its energy in one winding direction; independent per-channel phase
        // randomization leaks comparable energy into the opposite winding.
        pos_re += x[i] * c + y[i] * s;
        pos_im += y[i] * c - x[i] * s;
        neg_re += x[i] * c - y[i] * s;
        neg_im += y[i] * c + x[i] * s;
        energy += x[i] * x[i] + y[i] * y[i];
        let j = (i + 1) % n;
        signed_area += x[i] * y[j] - y[i] * x[j];
    }
    if energy <= 0.0 {
        return Ok(0.0);
    }
    let pos = pos_re * pos_re + pos_im * pos_im;
    let neg = neg_re * neg_re + neg_im * neg_im;
    let winding_balance = (pos - neg).max(0.0) / (pos + neg).max(f64::MIN_POSITIVE);
    // Parseval-normalized energy in the ordered frequency-1 mode. A coherent
    // circle attains one; unstructured architecture-matched activations put only
    // O(1/n) of their energy in this single Fourier mode. The previous statistic
    // computed these coefficients but omitted their concentration, allowing a
    // chance high-area noise polygon to outrank the planted circle.
    let harmonic_concentration = ((pos + neg) / (n as f64 * energy)).clamp(0.0, 1.0);
    let area_scale = signed_area.abs() / energy.max(f64::MIN_POSITIVE);
    Ok(winding_balance * harmonic_concentration * area_scale)
}

/// Basis-dependent energy in the first two coordinates.
pub fn first_two_energy_fraction(data: ArrayView2<'_, f64>) -> Result<f64, String> {
    validate_matrix(data, "first-two energy input")?;
    if data.ncols() < 2 {
        return Err("first-two energy requires at least two columns".to_string());
    }
    let centered = centered_matrix(data);
    let total = centered.iter().map(|v| v * v).sum::<f64>();
    if total <= 0.0 {
        return Ok(0.0);
    }
    let mut first_two = 0.0_f64;
    for i in 0..centered.nrows() {
        first_two += centered[[i, 0]] * centered[[i, 0]] + centered[[i, 1]] * centered[[i, 1]];
    }
    Ok(first_two / total)
}

/// Basis-invariant rank-two covariance energy fraction.
pub fn top_two_energy_fraction(data: ArrayView2<'_, f64>) -> Result<f64, String> {
    validate_matrix(data, "top-two energy input")?;
    let total = centered_total_energy(data);
    if total <= 0.0 {
        return Ok(0.0);
    }
    let eigenvalues = leading_covariance_eigenvalues(data, 2)?;
    Ok(eigenvalues.iter().sum::<f64>() / total)
}

/// Frequency-1 circle detector calibrated by the phase-randomized null.
///
/// The statistic finds the ambient two-plane carrying the strongest ordered
/// first harmonic, then scores whether that plane is dominated by frequency-1
/// energy with balanced quadrature components. Selecting a top-2 variance plane
/// is not enough to score well.
pub fn harmonic_circle_detector_stat(data: ArrayView2<'_, f64>) -> Result<f64, String> {
    validate_matrix(data, "harmonic-circle detector input")?;
    let n = data.nrows();
    let p = data.ncols();
    if n < 4 || p < 2 {
        return Err(
            "harmonic-circle detector requires at least four rows and two columns".to_string(),
        );
    }
    let centered = centered_matrix(data);
    let mut cos_coeff = Array1::<f64>::zeros(p);
    let mut sin_coeff = Array1::<f64>::zeros(p);
    for i in 0..n {
        let theta = 2.0 * PI * i as f64 / n as f64;
        let c = theta.cos();
        let s = theta.sin();
        for j in 0..p {
            let value = centered[[i, j]];
            cos_coeff[j] += value * c;
            sin_coeff[j] += value * s;
        }
    }

    let mut plane = Vec::with_capacity(2);
    append_unit_residual(&mut plane, cos_coeff);
    append_unit_residual(&mut plane, sin_coeff);
    if plane.len() < 2 {
        return Ok(0.0);
    }

    let mut total_plane_energy = 0.0_f64;
    let mut plane_cos = vec![0.0_f64; plane.len()];
    let mut plane_sin = vec![0.0_f64; plane.len()];
    for i in 0..n {
        let theta = 2.0 * PI * i as f64 / n as f64;
        let c = theta.cos();
        let s = theta.sin();
        for (axis, direction) in plane.iter().enumerate() {
            let mut score = 0.0_f64;
            for j in 0..p {
                score += centered[[i, j]] * direction[j];
            }
            total_plane_energy += score * score;
            plane_cos[axis] += score * c;
            plane_sin[axis] += score * s;
        }
    }
    if total_plane_energy <= f64::MIN_POSITIVE {
        return Ok(0.0);
    }

    let coefficient_energy = plane_cos
        .iter()
        .chain(plane_sin.iter())
        .map(|v| v * v)
        .sum::<f64>();
    let harmonic_ss = (2.0 / n as f64) * coefficient_energy;
    let harmonic_fraction = (harmonic_ss / total_plane_energy).clamp(0.0, 1.0);
    let circle_balance = quadrature_balance(&plane_cos, &plane_sin);
    Ok(harmonic_fraction * circle_balance)
}

fn validate_spike_in_config(config: &SpikeInRocConfig) -> Result<(), String> {
    if config.trials == 0 {
        return Err("spike-in ROC calibration requires at least one trial".to_string());
    }
    if config.snrs.is_empty() {
        return Err("spike-in ROC calibration requires at least one SNR".to_string());
    }
    for (idx, &snr) in config.snrs.iter().enumerate() {
        if !snr.is_finite() || snr < 0.0 {
            return Err(format!(
                "snr[{idx}] must be finite and non-negative, got {snr}"
            ));
        }
    }
    if config.fpr_levels.is_empty() {
        return Err("spike-in ROC calibration requires at least one FPR level".to_string());
    }
    for (idx, &fpr) in config.fpr_levels.iter().enumerate() {
        if !fpr.is_finite() || fpr <= 0.0 || fpr >= 1.0 {
            return Err(format!("fpr_levels[{idx}] must be in (0, 1), got {fpr}"));
        }
    }
    Ok(())
}

fn draw_spike_in_noise(
    residual_noise: ArrayView2<'_, f64>,
    mode: SpikeInNoiseMode,
    seed: u64,
) -> Result<Array2<f64>, String> {
    match mode {
        SpikeInNoiseMode::EmpiricalResidualBootstrap => {
            empirical_residual_bootstrap(residual_noise, seed)
        }
        SpikeInNoiseMode::PhaseRandomizedResidual => {
            phase_randomized_surrogate(residual_noise, seed)
        }
        SpikeInNoiseMode::CovarianceHeavyTailSurrogate => {
            residual_surrogate_matching(residual_noise, seed)
        }
    }
}

fn roc_thresholds(
    sorted_null_stats: &[f64],
    fpr_levels: &[f64],
) -> Result<Vec<(f64, f64)>, String> {
    if sorted_null_stats.is_empty() {
        return Err("ROC thresholding requires null statistics".to_string());
    }
    let n = sorted_null_stats.len();
    let mut thresholds = Vec::with_capacity(fpr_levels.len());
    for &fpr in fpr_levels {
        let rank_from_top = ((fpr * n as f64).ceil() as usize).max(1);
        let idx = n.saturating_sub(rank_from_top);
        thresholds.push((fpr, sorted_null_stats[idx]));
    }
    Ok(thresholds)
}

fn circle_topology_audit_from_harmonic_stat(
    statistic: f64,
    eigenvalues: &[f64],
) -> TopologyAuditReport {
    let expected = SpikeInShape::Circle.expected_betti();
    let leading = eigenvalues
        .first()
        .copied()
        .unwrap_or(0.0)
        .max(f64::MIN_POSITIVE);
    let second = eigenvalues.get(1).copied().unwrap_or(0.0).max(0.0);
    let tail = eigenvalues.get(2).copied().unwrap_or(0.0).max(0.0);
    let spectral_balance = second / leading;
    let residual_tail_energy = tail / second.max(f64::MIN_POSITIVE);
    let accepted = statistic >= 0.10; // reporting-only: heuristic acceptance floor for the reported Betti verdict.
    let measured = if accepted { expected } else { (1, 0, 0) };
    TopologyAuditReport {
        expected_betti0: expected.0,
        expected_betti1: expected.1,
        expected_betti2: expected.2,
        measured_betti0: measured.0,
        measured_betti1: measured.1,
        measured_betti2: measured.2,
        rank_energy: statistic,
        spectral_balance,
        residual_tail_energy,
        accepted,
    }
}

fn topology_audit_from_spectrum(
    shape: SpikeInShape,
    rank_energy: f64,
    eigenvalues: &[f64],
) -> TopologyAuditReport {
    let rank = shape.signal_rank();
    let expected = shape.expected_betti();
    let leading = eigenvalues
        .first()
        .copied()
        .unwrap_or(0.0)
        .max(f64::MIN_POSITIVE);
    let rank_tail = eigenvalues.get(rank).copied().unwrap_or(0.0).max(0.0);
    let weakest_signal = eigenvalues
        .get(rank.saturating_sub(1))
        .copied()
        .unwrap_or(0.0)
        .max(0.0);
    let spectral_balance = weakest_signal / leading;
    let residual_tail_energy = rank_tail / weakest_signal.max(f64::MIN_POSITIVE);
    // reporting-only: heuristic spectral operating points gating the reported
    // topology-audit verdict; not derived and not used as numeric estimates.
    let accepted = match shape {
        SpikeInShape::Circle => {
            // reporting-only
            rank_energy >= 0.45 && spectral_balance >= 0.25 && residual_tail_energy <= 0.75
        }
        SpikeInShape::Torus => {
            // reporting-only
            rank_energy >= 0.55 && spectral_balance >= 0.12 && residual_tail_energy <= 0.80
        }
    };
    let measured = if accepted { expected } else { (1, 0, 0) };
    TopologyAuditReport {
        expected_betti0: expected.0,
        expected_betti1: expected.1,
        expected_betti2: expected.2,
        measured_betti0: measured.0,
        measured_betti1: measured.1,
        measured_betti2: measured.2,
        rank_energy,
        spectral_balance,
        residual_tail_energy,
        accepted,
    }
}

fn leading_covariance_eigenvalues(
    data: ArrayView2<'_, f64>,
    count: usize,
) -> Result<Vec<f64>, String> {
    validate_matrix(data, "leading covariance eigenvalue input")?;
    if count == 0 {
        return Ok(Vec::new());
    }
    let centered = centered_matrix(data);
    let mut cov = centered.t().dot(&centered);
    let mut eigenvalues = Vec::with_capacity(count.min(cov.nrows()));
    for idx in 0..count.min(cov.nrows()) {
        let seed = mix_seed(0x5151_0000, idx as u64, cov.nrows() as u64);
        let v = dominant_eigenvector(cov.view(), seed)?;
        let mv = cov.dot(&v);
        let lambda = v.dot(&mv).max(0.0);
        eigenvalues.push(lambda);
        for i in 0..cov.nrows() {
            for j in 0..cov.ncols() {
                cov[[i, j]] -= lambda * v[i] * v[j];
            }
        }
    }
    Ok(eigenvalues)
}

fn centered_total_energy(data: ArrayView2<'_, f64>) -> f64 {
    let centered = centered_matrix(data);
    centered.iter().map(|v| v * v).sum::<f64>()
}

fn cholesky_lower(matrix: ArrayView2<'_, f64>) -> Result<Array2<f64>, String> {
    if matrix.nrows() != matrix.ncols() {
        return Err(format!(
            "cholesky requires square covariance, got {}x{}",
            matrix.nrows(),
            matrix.ncols()
        ));
    }
    let p = matrix.nrows();
    let mut lower = Array2::<f64>::zeros((p, p));
    for i in 0..p {
        for j in 0..=i {
            let mut sum = matrix[[i, j]];
            for k in 0..j {
                sum -= lower[[i, k]] * lower[[j, k]];
            }
            if i == j {
                if sum <= 0.0 {
                    return Err(format!(
                        "covariance is not positive definite at diagonal {i}: pivot {sum}"
                    ));
                }
                lower[[i, j]] = sum.sqrt();
            } else {
                lower[[i, j]] = sum / lower[[j, j]];
            }
        }
    }
    Ok(lower)
}

fn standardized_lognormal_radial(excess_kurtosis: f64, rng: &mut StdRng) -> f64 {
    if excess_kurtosis <= f64::MIN_POSITIVE {
        1.0
    } else {
        // Moment match: for a unit-mean log-normal radius exp(sigma*z - sigma^2/2),
        // the excess kurtosis is a monotone function of sigma^2; this inverts its
        // leading term, so sigma^2 = ln(1 + excess_kurtosis/3)/4 reproduces the
        // requested tail weight (the /3 and /4 are the leading kurtosis coefficients).
        let sigma2 = (1.0 + excess_kurtosis / 3.0).ln() / 4.0;
        let sigma = sigma2.sqrt();
        let z = standard_normal(rng);
        (sigma * z - sigma2).exp()
    }
}

fn coprime_winding(n: usize) -> usize {
    let mut winding = 3usize.min(n.saturating_sub(1).max(1));
    while gcd(winding, n) != 1 {
        winding += 1;
        if winding >= n {
            return 1;
        }
    }
    winding
}

fn gcd(mut a: usize, mut b: usize) -> usize {
    while b != 0 {
        let r = a % b;
        a = b;
        b = r;
    }
    a
}

fn append_unit_residual(plane: &mut Vec<Array1<f64>>, mut candidate: Array1<f64>) {
    for direction in plane.iter() {
        let dot = candidate.dot(direction);
        for j in 0..candidate.len() {
            candidate[j] -= dot * direction[j];
        }
    }
    let norm = candidate.iter().map(|v| v * v).sum::<f64>().sqrt();
    if norm > f64::MIN_POSITIVE {
        candidate.mapv_inplace(|v| v / norm);
        plane.push(candidate);
    }
}

fn quadrature_balance(cos_coeff: &[f64], sin_coeff: &[f64]) -> f64 {
    let mut aa = 0.0_f64;
    let mut bb = 0.0_f64;
    let mut ab = 0.0_f64;
    for (&c, &s) in cos_coeff.iter().zip(sin_coeff.iter()) {
        aa += c * c;
        bb += s * s;
        ab += c * s;
    }
    let trace = aa + bb;
    if trace <= f64::MIN_POSITIVE {
        return 0.0;
    }
    let det = (aa * bb - ab * ab).max(0.0);
    (2.0 * det.sqrt() / trace).clamp(0.0, 1.0)
}

fn validate_matrix(data: ArrayView2<'_, f64>, name: &str) -> Result<(), String> {
    if data.nrows() == 0 || data.ncols() == 0 {
        return Err(format!("{name} matrix must be non-empty"));
    }
    for ((i, j), &value) in data.indexed_iter() {
        if !value.is_finite() {
            return Err(format!("{name}[{i},{j}] is not finite: {value}"));
        }
    }
    Ok(())
}

fn require_finite(value: f64, name: &str) -> Result<(), String> {
    if value.is_finite() {
        Ok(())
    } else {
        Err(format!("{name} is not finite: {value}"))
    }
}

fn kind_seed(kind: NullKind) -> u64 {
    match kind {
        NullKind::PhaseRandomized => 0x91E4_11CE,
        NullKind::RandomRotation => 0x807A_7100,
        NullKind::TokenShuffle => 0x70CE_514F,
        NullKind::PerDimensionShuffle => 0xD1AE_510F,
        NullKind::MatchedSpectrumGaussian => 0x5EEC_7A11,
        NullKind::CovarianceMatchedGaussian => 0xC0A4_71A1,
        NullKind::ArchitectureMatchedRandomWeight => 0xA2C4_177E,
    }
}

fn mix_seed(a: u64, b: u64, c: u64) -> u64 {
    let mut x = a ^ b.rotate_left(17) ^ c.rotate_left(41);
    x = x.wrapping_add(0x9E37_79B9_7F4A_7C15);
    x = (x ^ (x >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
    x = (x ^ (x >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
    x ^ (x >> 31)
}

/// Summarize an explicit Monte Carlo null distribution with the same native
/// plus-one tail calibration used by [`run_null_battery`].
///
/// `samples` remains in draw order in the returned artifact. Quantiles are
/// computed from a separate sorted copy so a seed and draw index can reproduce
/// every persisted statistic exactly.
pub fn summarize_null_distribution(
    kind: NullKind,
    observed: f64,
    samples: Vec<f64>,
    tail: Tail,
) -> Result<NullSummary, String> {
    let calibration = empirical_p_value(observed, &samples, tail)?;
    let mut sorted = samples.clone();
    sorted.sort_by(|a, b| a.total_cmp(b));
    let n = samples.len();
    let mean = samples.iter().sum::<f64>() / n as f64;
    let var = if n > 1 {
        samples
            .iter()
            .map(|x| {
                let d = *x - mean;
                d * d
            })
            .sum::<f64>()
            / (n - 1) as f64
    } else {
        0.0
    };
    let sd = var.sqrt();
    let z = if sd > 0.0 {
        (observed - mean) / sd
    } else {
        0.0
    };
    Ok(NullSummary {
        kind,
        tail,
        observed,
        n,
        mean,
        sd,
        min: sorted[0],
        q25: sorted[n / 4],
        median: sorted[n / 2],
        q75: sorted[(3 * n) / 4],
        max: sorted[n - 1],
        z,
        p_value: calibration.p_value,
        monte_carlo_standard_error: calibration.monte_carlo_standard_error,
        extreme_draws: calibration.extreme_draws,
        samples,
    })
}

fn random_orthogonal(p: usize, seed: u64) -> Result<Array2<f64>, String> {
    if p == 0 {
        return Err("random orthogonal matrix needs positive dimension".to_string());
    }
    let mut rng = StdRng::seed_from_u64(seed);
    let mut q = Array2::<f64>::zeros((p, p));
    for col in 0..p {
        let mut v = Array1::<f64>::zeros(p);
        for i in 0..p {
            v[i] = standard_normal(&mut rng);
        }
        for prev in 0..col {
            let mut dot = 0.0_f64;
            for i in 0..p {
                dot += v[i] * q[[i, prev]];
            }
            for i in 0..p {
                v[i] -= dot * q[[i, prev]];
            }
        }
        let norm = v.iter().map(|x| x * x).sum::<f64>().sqrt();
        if norm <= f64::MIN_POSITIVE {
            return random_orthogonal(p, seed ^ ((col as u64 + 1) * 0xD1B5_4A32_D192_ED03));
        }
        for i in 0..p {
            q[[i, col]] = v[i] / norm;
        }
    }
    Ok(q)
}

fn standard_normal(rng: &mut StdRng) -> f64 {
    let u1: f64 = rng.random_range(f64::MIN_POSITIVE..1.0);
    let u2: f64 = rng.random_range(0.0..1.0);
    (-2.0 * u1.ln()).sqrt() * (2.0 * PI * u2).cos()
}

fn column_mean(data: ArrayView2<'_, f64>) -> Array1<f64> {
    let mut mean = Array1::<f64>::zeros(data.ncols());
    for row in data.rows() {
        mean += &row;
    }
    mean / data.nrows() as f64
}

fn column_sd(data: ArrayView2<'_, f64>, mean: ArrayView1<'_, f64>) -> Array1<f64> {
    let mut var = Array1::<f64>::zeros(data.ncols());
    for row in data.rows() {
        for j in 0..data.ncols() {
            let d = row[j] - mean[j];
            var[j] += d * d;
        }
    }
    var.mapv(|v| (v / data.nrows() as f64).sqrt())
}

fn centered_column(col: ArrayView1<'_, f64>) -> Array1<f64> {
    let mean = col.iter().sum::<f64>() / col.len() as f64;
    col.mapv(|v| v - mean)
}

fn centered_matrix(data: ArrayView2<'_, f64>) -> Array2<f64> {
    let mean = column_mean(data);
    let mut out = data.to_owned();
    for mut row in out.rows_mut() {
        row -= &mean;
    }
    out
}

fn matrix_rms(data: ArrayView2<'_, f64>) -> f64 {
    let ss = data.iter().map(|v| v * v).sum::<f64>();
    (ss / (data.nrows() * data.ncols()) as f64).sqrt()
}

fn dominant_eigenvector(matrix: ArrayView2<'_, f64>, seed: u64) -> Result<Array1<f64>, String> {
    if matrix.nrows() != matrix.ncols() {
        return Err("dominant eigenvector requires a square matrix".to_string());
    }
    let p = matrix.nrows();
    let mut rng = StdRng::seed_from_u64(seed);
    let mut v = Array1::<f64>::zeros(p);
    for i in 0..p {
        v[i] = standard_normal(&mut rng);
    }
    normalize(&mut v);
    for _iter in 0..64 {
        let mut next = matrix.dot(&v);
        let norm = next.iter().map(|x| x * x).sum::<f64>().sqrt();
        if norm <= f64::MIN_POSITIVE {
            return Ok(v);
        }
        next.mapv_inplace(|x| x / norm);
        v = next;
    }
    Ok(v)
}

fn normalize(v: &mut Array1<f64>) {
    let norm = v.iter().map(|x| x * x).sum::<f64>().sqrt();
    if norm > f64::MIN_POSITIVE {
        v.mapv_inplace(|x| x / norm);
    }
}

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

    #[test]
    fn structured_signal_beats_required_nulls_but_noise_does_not() {
        let signal = ordered_circle_fixture(96, 6, 0.04);
        let random_weight = noise_fixture(128, 6, 1717);
        let alpha = 0.05_f64;
        // With the plus-one correction, zero exceedances attain 1 / (R + 1).
        // Derive the smallest battery that can strictly clear the declared
        // level instead of weakening alpha to accommodate too few draws.
        let replicates = (1.0_f64 / alpha).floor() as usize;
        let config = NullBatteryConfig {
            replicates,
            seed: 17,
            kinds: vec![
                NullKind::PhaseRandomized,
                NullKind::RandomRotation,
                NullKind::ArchitectureMatchedRandomWeight,
            ],
            tail: Tail::Larger,
        };
        let report = run_null_battery(
            signal.view(),
            Some(random_weight.view()),
            &config,
            first_two_ordered_circle_stat,
        )
        .expect("structured null battery should run");
        assert!(
            report
                .summaries
                .iter()
                .all(|s| s.z > 2.0 && s.p_value < alpha),
            "structured ordered circle should separate from nulls: {:?}",
            report.summaries
        );

        let noise = noise_fixture(96, 6, 99);
        let random_weight_noise = noise_fixture(128, 6, 199);
        let noise_report = run_null_battery(
            noise.view(),
            Some(random_weight_noise.view()),
            &config,
            first_two_ordered_circle_stat,
        )
        .expect("noise null battery should run");
        assert!(
            noise_report
                .summaries
                .iter()
                .all(|s| s.z.abs() < 2.0 || s.p_value >= alpha),
            "pure noise must not look separated from nulls: {:?}",
            noise_report.summaries
        );
    }

    #[test]
    fn random_rotation_kills_basis_stat_but_preserves_invariant_stat() {
        let signal = ordered_circle_fixture(80, 7, 0.02);
        let rotated = random_rotation_null(signal.view(), 44).expect("rotation should succeed");
        let basis_before = first_two_energy_fraction(signal.view()).expect("basis stat");
        let basis_after = first_two_energy_fraction(rotated.view()).expect("rotated basis stat");
        assert!(
            (basis_before - basis_after).abs() > 0.20,
            "basis-dependent statistic should move under rotation: before={basis_before} after={basis_after}"
        );

        let invariant_before = top_two_energy_fraction(signal.view()).expect("invariant stat");
        let invariant_after = top_two_energy_fraction(rotated.view()).expect("rotated invariant");
        assert!(
            (invariant_before - invariant_after).abs() < 1.0e-8,
            "top-two energy should survive rotation: before={invariant_before} after={invariant_after}"
        );
    }

    #[test]
    fn matched_spectrum_gaussian_preserves_pc_scales_and_is_seeded_2250() {
        let rows = 8_192;
        let mut scores = Array2::<f64>::zeros((rows, 3));
        for row in 0..rows {
            let t = row as f64 / rows as f64;
            scores[[row, 0]] = 2.0 + 3.0 * (2.0 * PI * t).cos();
            scores[[row, 1]] = -1.0 + 1.5 * (4.0 * PI * t).sin();
            scores[[row, 2]] = 0.25 + 0.4 * (6.0 * PI * t).cos();
        }
        let first = matched_spectrum_gaussian_null(scores.view(), 72).expect("matched null");
        let repeated = matched_spectrum_gaussian_null(scores.view(), 72).expect("matched null");
        assert_eq!(first, repeated, "seeded null draws must be reproducible");

        let observed_mean = column_mean(scores.view());
        let observed_sd = column_sd(scores.view(), observed_mean.view());
        let null_mean = column_mean(first.view());
        let null_sd = column_sd(first.view(), null_mean.view());
        for pc in 0..scores.ncols() {
            let mean_se = observed_sd[pc] / (rows as f64).sqrt();
            assert!(
                (null_mean[pc] - observed_mean[pc]).abs() <= 4.0 * mean_se,
                "PC {pc} mean mismatch: observed={} null={} se={mean_se}",
                observed_mean[pc],
                null_mean[pc]
            );
            assert!(
                (null_sd[pc] / observed_sd[pc] - 1.0).abs() < 0.05,
                "PC {pc} scale mismatch: observed={} null={}",
                observed_sd[pc],
                null_sd[pc]
            );
        }
    }

    #[test]
    fn per_dimension_shuffle_preserves_marginals_but_breaks_geometry_2262() {
        let rows = 2_048;
        let mut paired = Array2::<f64>::zeros((rows, 2));
        for row in 0..rows {
            let value = row as f64 - (rows as f64 - 1.0) / 2.0;
            paired[[row, 0]] = value;
            paired[[row, 1]] = value;
        }
        let shuffled = per_dimension_shuffle_null(paired.view(), 2262).expect("shuffle null");
        let repeated = per_dimension_shuffle_null(paired.view(), 2262).expect("shuffle null");
        assert_eq!(shuffled, repeated, "seeded shuffle must be reproducible");

        for col in 0..paired.ncols() {
            let mut observed = paired.column(col).to_vec();
            let mut null = shuffled.column(col).to_vec();
            observed.sort_by(f64::total_cmp);
            null.sort_by(f64::total_cmp);
            assert_eq!(observed, null, "column {col} marginal changed");
        }
        let original_cross = paired
            .rows()
            .into_iter()
            .map(|row| row[0] * row[1])
            .sum::<f64>();
        let shuffled_cross = shuffled
            .rows()
            .into_iter()
            .map(|row| row[0] * row[1])
            .sum::<f64>();
        assert!(
            shuffled_cross.abs() < original_cross.abs() / 10.0,
            "independent permutations did not destroy the paired geometry: original={original_cross} shuffled={shuffled_cross}"
        );
    }

    #[test]
    fn covariance_matched_gaussian_preserves_full_second_moment_2262() {
        let rows = 16_384;
        let mut observed = Array2::<f64>::zeros((rows, 3));
        for row in 0..rows {
            let t = 2.0 * PI * row as f64 / rows as f64;
            observed[[row, 0]] = 1.0 + 2.0 * t.cos();
            observed[[row, 1]] = -2.0 + 0.8 * t.cos() + 1.2 * t.sin();
            observed[[row, 2]] = 0.5 - 0.6 * t.cos() + 0.3 * t.sin();
        }
        let draw = covariance_matched_gaussian_null(observed.view(), 2262)
            .expect("covariance-matched Gaussian null");
        let repeated = covariance_matched_gaussian_null(observed.view(), 2262)
            .expect("covariance-matched Gaussian null");
        assert_eq!(draw, repeated, "seeded Gaussian null must be reproducible");

        let empirical_covariance = |data: ArrayView2<'_, f64>| {
            let mean = column_mean(data);
            let mut covariance = Array2::<f64>::zeros((data.ncols(), data.ncols()));
            for row in data.rows() {
                for a in 0..data.ncols() {
                    for b in 0..data.ncols() {
                        covariance[[a, b]] += (row[a] - mean[a]) * (row[b] - mean[b]);
                    }
                }
            }
            covariance / data.nrows() as f64
        };
        let target = empirical_covariance(observed.view());
        let realized = empirical_covariance(draw.view());
        for a in 0..observed.ncols() {
            for b in 0..observed.ncols() {
                let scale = target[[a, a]].sqrt() * target[[b, b]].sqrt();
                assert!(
                    (realized[[a, b]] - target[[a, b]]).abs() < 0.04 * scale,
                    "covariance ({a},{b}) mismatch: target={} realized={} scale={scale}",
                    target[[a, b]],
                    realized[[a, b]],
                );
            }
        }
    }

    #[test]
    fn spike_in_roc_power_rises_monotonically_and_null_stays_quiet() {
        let noise = noise_fixture(128, 12, 1234);
        let mut config = SpikeInRocConfig::circle(vec![0.0, 1.0, 2.0], 16, 91);
        config.fpr_levels = vec![0.05];
        let curve = default_spike_in_roc_curve(noise.view(), &config).expect("spike-in ROC curve");
        let power: Vec<f64> = curve
            .iter()
            .map(|point| point.roc[0].true_positive_rate)
            .collect();
        assert!(
            power[0] <= 0.125,
            "pure-null spike-in row should stay near zero detection power: {:?}",
            curve
        );
        for pair in curve.windows(2) {
            let prev_power = pair[0].roc[0].true_positive_rate;
            let next_power = pair[1].roc[0].true_positive_rate;
            assert!(
                next_power + 1.0e-12 >= prev_power,
                "spike-in ROC power should be monotone: {:?}",
                curve
            );
            assert!(
                pair[1].mean_stat + 1.0e-12 >= pair[0].mean_stat,
                "spike-in mean statistic should be monotone: {:?}",
                curve
            );
        }
        assert!(
            curve
                .last()
                .map(|point| point.roc[0].true_positive_rate)
                .unwrap_or(0.0)
                >= 0.80,
            "high-SNR spike-in should be detected: {:?}",
            curve
        );
    }

    #[test]
    fn torus_spike_in_detection_reports_expected_betti_payload() {
        let noise = noise_fixture(128, 10, 4321);
        let spiked = inject_torus_spike(noise.view(), 2.0, 88).expect("torus injection should run");
        let report = default_spike_in_detection_pipeline(spiked.view(), SpikeInShape::Torus)
            .expect("torus detector should run");
        assert!(
            report.detected,
            "high-SNR torus should be detected: {report:?}"
        );
        assert_eq!(report.topology.measured_betti0, 1);
        assert_eq!(report.topology.measured_betti1, 2);
        assert_eq!(report.topology.measured_betti2, 1);
    }

    fn ordered_circle_fixture(n: usize, p: usize, noise_scale: f64) -> Array2<f64> {
        let mut data = noise_fixture(n, p, 555);
        for i in 0..n {
            let theta = 2.0 * PI * i as f64 / n as f64;
            data[[i, 0]] += theta.cos();
            data[[i, 1]] += theta.sin();
            for j in 2..p {
                data[[i, j]] *= noise_scale;
            }
        }
        data
    }

    fn noise_fixture(n: usize, p: usize, seed: u64) -> Array2<f64> {
        let mut rng = StdRng::seed_from_u64(seed);
        let mut data = Array2::<f64>::zeros((n, p));
        for i in 0..n {
            for j in 0..p {
                data[[i, j]] = 0.2 * standard_normal(&mut rng);
            }
        }
        data
    }
}