Aura/crates/aura-analysis/src/lib.rs

//! Domain-free post-run statistics and selection provenance (C1/C12): the
//! multiple-comparison hurdle math (`inv_norm_cdf`, `expected_max_of_normals`)
//! and the selection-provenance record a sweep winner carries
//! (`FamilySelection` / `SelectionMode`, embedded by the engine's
//! `RunManifest.selection`). No trading vocabulary lives here: the backtest
//! reductions (R-metrics, the position-event table) moved to
//! `aura-backtest::metrics` (issue #291, C28 phase 5). Foundation-grade — any
//! ladder rung may depend on this crate without violating the C28 import
//! direction. Every type's serde shape is byte-pinned by C18 goldens; bodies
//! moved verbatim.

/// Which selection objective produced the record (additive provenance, C23).
/// `Argmax` is the bare-best pick (cycle 0076), deflated for the number of trials.
/// `PlateauMean` / `PlateauWorst` (cycle 0077) argmax the neighbourhood-smoothed
/// surface instead, so peak-vs-plateau runs stay distinguishable on this field.
#[derive(Clone, Copy, Debug, PartialEq, serde::Serialize, serde::Deserialize)]
pub enum SelectionMode {
    Argmax,
    PlateauMean,
    PlateauWorst,
}

/// Selection-provenance for a sweep winner. The selection RULE (`mode`) and its
/// ANNOTATION are orthogonal: `Argmax` carries the trials-deflation annotation
/// (`deflated_score` / `overfit_probability` / `n_resamples` / `block_len` /
/// `seed`); `Plateau*` carries the smoothing annotation (`neighbourhood_score` /
/// `n_neighbours`). Each annotation is `Option`, present iff its rule produced the
/// record; all are additive — recorded, never re-ranking (C23). A legacy line
/// (pre-0077) carries the deflation scalars as bare values that deserialize to
/// `Some` (serde default), so it still loads (C14/C18).
#[derive(Clone, Debug, PartialEq, serde::Serialize, serde::Deserialize)]
pub struct FamilySelection {
    pub selection_metric: String,
    pub n_trials: usize,
    pub raw_winner_metric: f64,
    pub mode: SelectionMode,
    // deflation annotation (present iff mode == Argmax with a deflation run)
    #[serde(default, skip_serializing_if = "Option::is_none")]
    pub deflated_score: Option<f64>,
    #[serde(default, skip_serializing_if = "Option::is_none")]
    pub overfit_probability: Option<f64>,
    #[serde(default, skip_serializing_if = "Option::is_none")]
    pub n_resamples: Option<usize>,
    #[serde(default, skip_serializing_if = "Option::is_none")]
    pub block_len: Option<usize>,
    #[serde(default, skip_serializing_if = "Option::is_none")]
    pub seed: Option<u64>,
    // plateau annotation (present iff mode is Plateau*)
    #[serde(default, skip_serializing_if = "Option::is_none")]
    pub neighbourhood_score: Option<f64>,
    #[serde(default, skip_serializing_if = "Option::is_none")]
    pub n_neighbours: Option<usize>,
}

/// Inverse standard-normal CDF (quantile function), Acklam's rational
/// approximation — absolute error < ~1.2e-9 over `p ∈ (0,1)`. Pure (C1).
/// Callers pass strictly-interior `p` (`k >= 2` keeps the argument off the
/// boundaries); the boundary branches return ±inf as a defined limit.
// Acklam's published coefficients verbatim; some carry more decimals than an
// f64 holds (clippy::excessive_precision) — kept literal as reference constants.
#[allow(clippy::excessive_precision)]
pub fn inv_norm_cdf(p: f64) -> f64 {
    const A: [f64; 6] = [-3.969683028665376e+01, 2.209460984245205e+02, -2.759285104469687e+02,
                          1.383577518672690e+02, -3.066479806614716e+01, 2.506628277459239e+00];
    const B: [f64; 5] = [-5.447609879822406e+01, 1.615858368580409e+02, -1.556989798598866e+02,
                          6.680131188771972e+01, -1.328068155288572e+01];
    const C: [f64; 6] = [-7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e+00,
                         -2.549732539343734e+00, 4.374664141464968e+00, 2.938163982698783e+00];
    const D: [f64; 4] = [7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e+00,
                         3.754408661907416e+00];
    const P_LOW: f64 = 0.02425;
    if p <= 0.0 { return f64::NEG_INFINITY; }
    if p >= 1.0 { return f64::INFINITY; }
    if p < P_LOW {
        let q = (-2.0 * p.ln()).sqrt();
        (((((C[0]*q + C[1])*q + C[2])*q + C[3])*q + C[4])*q + C[5])
            / ((((D[0]*q + D[1])*q + D[2])*q + D[3])*q + 1.0)
    } else if p <= 1.0 - P_LOW {
        let q = p - 0.5;
        let r = q * q;
        (((((A[0]*r + A[1])*r + A[2])*r + A[3])*r + A[4])*r + A[5]) * q
            / (((((B[0]*r + B[1])*r + B[2])*r + B[3])*r + B[4])*r + 1.0)
    } else {
        let q = (-2.0 * (1.0 - p).ln()).sqrt();
        -(((((C[0]*q + C[1])*q + C[2])*q + C[3])*q + C[4])*q + C[5])
            / ((((D[0]*q + D[1])*q + D[2])*q + D[3])*q + 1.0)
    }
}

/// Expected maximum of `k` i.i.d. standard normals — the hurdle a search of `k`
/// configurations clears by chance alone. `k <= 1` -> `0.0` (a single trial has no
/// multiple-comparison inflation, and guards the `Φ⁻¹(1 - 1/k) -> Φ⁻¹(0)`
/// divergence). For `k >= 2`:
/// `(1 - γ)·Φ⁻¹(1 - 1/k) + γ·Φ⁻¹(1 - 1/(k·e))`, `γ = 0.5772156649…`. Pure (C1).
pub fn expected_max_of_normals(k: usize) -> f64 {
    if k <= 1 {
        return 0.0;
    }
    const GAMMA: f64 = 0.577_215_664_901_532_9;
    let kf = k as f64;
    (1.0 - GAMMA) * inv_norm_cdf(1.0 - 1.0 / kf)
        + GAMMA * inv_norm_cdf(1.0 - 1.0 / (kf * std::f64::consts::E))
}

/// A tiny, fully-deterministic, dependency-free PRNG (SplitMix64). The seed
/// completely determines the sequence; no external entropy, no global state.
/// Bit-stable across toolchains and crate versions — the property C1 needs for
/// seed-as-input reproducibility (C12).
pub struct SplitMix64 {
    state: u64,
}

impl SplitMix64 {
    pub fn new(seed: u64) -> Self {
        Self { state: seed }
    }
    pub fn next_u64(&mut self) -> u64 {
        self.state = self.state.wrapping_add(0x9E37_79B9_7F4A_7C15);
        let mut z = self.state;
        z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
        z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
        z ^ (z >> 31)
    }
    /// A uniform `f64` in `[0, 1)` from the top 53 bits.
    pub fn next_f64(&mut self) -> f64 {
        (self.next_u64() >> 11) as f64 / ((1u64 << 53) as f64)
    }
}

/// Mean + a fixed quantile set of one metric across the realizations. `p50` is
/// the median (the two coincide by definition), so "mean/median/quantiles" is
/// `mean` + `p50` + the surrounding quantiles, with no redundant field.
#[derive(Clone, Debug, PartialEq, serde::Serialize, serde::Deserialize)]
pub struct MetricStats {
    pub mean: f64,
    pub p5: f64,
    pub p25: f64,
    pub p50: f64, // == median
    pub p75: f64,
    pub p95: f64,
}

impl MetricStats {
    /// Mean + the fixed type-7 quantile set over a value set. Sorts a copy;
    /// `values` must be finite and non-empty. The shared reduction behind both the
    /// MC aggregate (a metric across draws) and the walk-forward param-stability
    /// summary (a param across windows; `aura-engine`'s `param_stability`).
    pub fn from_values(values: &[f64]) -> MetricStats {
        let mut xs = values.to_vec();
        xs.sort_by(|a, b| a.partial_cmp(b).expect("values are finite"));
        MetricStats {
            mean: xs.iter().sum::<f64>() / xs.len() as f64,
            p5: quantile(&xs, 0.05),
            p25: quantile(&xs, 0.25),
            p50: quantile(&xs, 0.50),
            p75: quantile(&xs, 0.75),
            p95: quantile(&xs, 0.95),
        }
    }
}

/// Linear-interpolation quantile (the numpy/Excel "type 7" default) over a
/// pre-sorted, non-empty slice. `p` in `[0, 1]`. `rank = p * (n - 1)`;
/// interpolate between the two bracketing order statistics. `n == 1` returns the
/// sole value.
pub fn quantile(sorted: &[f64], p: f64) -> f64 {
    let n = sorted.len();
    if n == 1 {
        return sorted[0];
    }
    let rank = p * (n - 1) as f64;
    let lo = rank.floor() as usize;
    let frac = rank - lo as f64;
    if lo + 1 < n {
        sorted[lo] + frac * (sorted[lo + 1] - sorted[lo])
    } else {
        sorted[n - 1]
    }
}

/// One moving-block resample of `rs` to length `rs.len()` (non-circular; the final
/// block is truncated so the resample has exactly `n` values). `block_len` must be
/// pre-clamped to `[1, rs.len()]` by the caller. Pure given the `rng` state — the
/// shared kernel behind `r_bootstrap` and the trials-deflation reality-check (C1).
pub fn resample_block(rs: &[f64], block_len: usize, rng: &mut SplitMix64) -> Vec<f64> {
    let n = rs.len();
    debug_assert!(
        (1..=n).contains(&block_len),
        "resample_block requires block_len in [1, rs.len()] (caller pre-clamps); \
         got block_len={block_len}, rs.len()={n}",
    );
    let mut sample: Vec<f64> = Vec::with_capacity(n);
    while sample.len() < n {
        let start = (rng.next_u64() % (n - block_len + 1) as u64) as usize;
        let take = block_len.min(n - sample.len());
        sample.extend_from_slice(&rs[start..start + take]);
    }
    sample
}

/// Pearson product-moment correlation of two equal-length finite series.
/// `n < 2`, unequal lengths, or zero variance on either side → `0.0`: no linear
/// relationship is defined, and `0.0` is the honest "no correlation" value (mirrors
/// how the R bootstrap floors a degenerate series rather than emitting `NaN`).
pub fn pearson_corr(xs: &[f64], ys: &[f64]) -> f64 {
    let n = xs.len();
    if n < 2 || ys.len() != n {
        return 0.0;
    }
    let nf = n as f64;
    let mx = xs.iter().sum::<f64>() / nf;
    let my = ys.iter().sum::<f64>() / nf;
    let (mut sxy, mut sxx, mut syy) = (0.0, 0.0, 0.0);
    for i in 0..n {
        let (dx, dy) = (xs[i] - mx, ys[i] - my);
        sxy += dx * dy;
        sxx += dx * dx;
        syy += dy * dy;
    }
    if sxx <= 0.0 || syy <= 0.0 {
        return 0.0;
    }
    sxy / (sxx.sqrt() * syy.sqrt())
}

/// In-place Fisher-Yates shuffle driven by `SplitMix64` — a uniform permutation,
/// sampling WITHOUT replacement (the permutation null). Distinct from
/// `resample_block`, which resamples contiguous blocks WITH replacement. Pure given
/// the `rng` state (C1).
pub fn permute<T>(xs: &mut [T], rng: &mut SplitMix64) {
    for i in (1..xs.len()).rev() {
        let j = (rng.next_u64() % (i as u64 + 1)) as usize;
        xs.swap(i, j);
    }
}

/// Per-metric vocabulary a report payload supplies to the registry's
/// ranking/deflation machinery (#147). Narrow by design (#136): exactly the
/// surface the two production implementors — the backtest metrics
/// (`RunMetrics`, centred moving-block null) and the measurement IC
/// (within-run permutation null) — demonstrably share. The registry owns the
/// family fold — best-of-K, p95, Laplace, dispersion floor — the payload
/// owns what each metric IS: its name, direction, value read, and what one
/// draw of its null statistic means.
pub trait MetricVocabulary: Sized {
    /// Resolved metric identity — cheap to copy, resolved once per call.
    type Key: Copy;
    /// Name → key; `None` for a name outside this vocabulary.
    fn resolve(name: &str) -> Option<Self::Key>;
    /// The canonical accepted-name roster, in refusal-prose order.
    fn known() -> &'static [&'static str];
    /// The key's optimisation direction: `true` when larger is better.
    fn higher_is_better(key: Self::Key) -> bool;
    /// The member's scalar under the key. A missing optional block reads
    /// `f64::NEG_INFINITY` (the worst rank for higher-is-better keys).
    fn value(&self, key: Self::Key) -> f64;
    /// Whether the key deflates via a resampled null (the best-of-K arm) or
    /// the dispersion floor.
    fn has_resampling_null(key: Self::Key) -> bool;
    /// ONE draw of this member's null statistic under H0, advancing `rng`.
    /// `None` means no usable null input for this member, and consumes no
    /// rng. `block_len` parameterises resampling vocabularies; permutation
    /// vocabularies ignore it.
    fn null_draw(&self, key: Self::Key, block_len: usize, rng: &mut SplitMix64) -> Option<f64>;
}

/// One-sided Laplace-smoothed tail estimate `(ge + 1) / (n + 1)`: the
/// probability mass at-or-above an observed statistic given that `ge` of `n`
/// null draws reached it. The single source of the formula: the registry's
/// R deflation arm and the CLI's IC reduction both call it.
pub fn one_sided_p_laplace(ge: usize, n: usize) -> f64 {
    (ge + 1) as f64 / (n + 1) as f64
}

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

    #[test]
    fn pearson_corr_known_values() {
        // identical series → +1; reversed monotone → −1; constant side → 0 (zero variance)
        assert!((pearson_corr(&[1.0, 2.0, 3.0, 4.0], &[1.0, 2.0, 3.0, 4.0]) - 1.0).abs() < 1e-12);
        assert!((pearson_corr(&[1.0, 2.0, 3.0, 4.0], &[4.0, 3.0, 2.0, 1.0]) + 1.0).abs() < 1e-12);
        assert_eq!(pearson_corr(&[1.0, 2.0, 3.0, 4.0], &[7.0, 7.0, 7.0, 7.0]), 0.0);
        // degenerate: fewer than two points, or unequal lengths → 0.0
        assert_eq!(pearson_corr(&[1.0], &[1.0]), 0.0);
        assert_eq!(pearson_corr(&[1.0, 2.0], &[1.0]), 0.0);
    }

    #[test]
    fn permute_is_a_permutation_and_deterministic() {
        let mut a = [1u32, 2, 3, 4, 5, 6, 7, 8];
        let mut b = a;
        let mut ra = SplitMix64::new(42);
        let mut rb = SplitMix64::new(42);
        permute(&mut a, &mut ra);
        permute(&mut b, &mut rb);
        assert_eq!(a, b, "same seed → same permutation (C1)");
        let mut sorted = a;
        sorted.sort_unstable();
        assert_eq!(sorted, [1, 2, 3, 4, 5, 6, 7, 8], "a permutation preserves the multiset");
        // a different seed generally yields a different order (not a hard guarantee, but true here)
        let mut c = [1u32, 2, 3, 4, 5, 6, 7, 8];
        let mut rc = SplitMix64::new(43);
        permute(&mut c, &mut rc);
        assert_ne!(a, c);
    }

    #[test]
    fn inv_norm_cdf_matches_known_quantiles() {
        assert!((inv_norm_cdf(0.975) - 1.959964).abs() < 1e-3);
        assert!(inv_norm_cdf(0.5).abs() < 1e-9);
        assert!((inv_norm_cdf(0.1) + inv_norm_cdf(0.9)).abs() < 1e-9); // symmetry
    }

    #[test]
    fn expected_max_of_normals_is_zero_at_one_and_monotone() {
        assert_eq!(expected_max_of_normals(1), 0.0);
        assert_eq!(expected_max_of_normals(0), 0.0);
        let (e2, e4, e10, e100) = (
            expected_max_of_normals(2), expected_max_of_normals(4),
            expected_max_of_normals(10), expected_max_of_normals(100),
        );
        assert!(e2 < e4 && e4 < e10 && e10 < e100); // strictly increasing
        assert!((1.0..1.1).contains(&e4));           // ~1.05 (not the √(2 ln 4)=1.66 asymptote)
        assert!((2.4..2.7).contains(&e100));         // ~2.53
    }

    #[test]
    fn quantile_endpoints_and_singleton() {
        // singleton returns the sole value; p=0 -> min, p=1 -> max.
        assert_eq!(quantile(&[42.0], 0.0), 42.0);
        assert_eq!(quantile(&[42.0], 0.5), 42.0);
        assert_eq!(quantile(&[42.0], 1.0), 42.0);
        let xs = [1.0, 2.0, 3.0, 4.0, 5.0];
        assert_eq!(quantile(&xs, 0.0), 1.0);
        assert_eq!(quantile(&xs, 1.0), 5.0);
    }

    #[test]
    fn metric_stats_from_values_matches_known_fixture() {
        // type-7 quantile + mean over [0,1,2,3,4], directly on the
        // extracted reduction: mean=2.0, p50=2.0, p5≈0.2, p95≈3.8 (same numbers the
        // aggregate fixture pins, now on from_values).
        let s = MetricStats::from_values(&[0.0, 1.0, 2.0, 3.0, 4.0]);
        assert_eq!(s.mean, 2.0);
        assert_eq!(s.p50, 2.0);
        assert!((s.p5 - 0.2).abs() < 1e-9, "p5 = {}", s.p5);
        assert!((s.p95 - 3.8).abs() < 1e-9, "p95 = {}", s.p95);
    }

    #[test]
    fn metric_stats_serde_round_trips() {
        // MetricStats gains serde (consistent with the
        // report types) so the CLI summary renders it and #70 lineage can persist.
        let s = MetricStats::from_values(&[1.0, 2.0, 3.0]);
        let json = serde_json::to_string(&s).expect("serialize MetricStats");
        let back: MetricStats = serde_json::from_str(&json).expect("deserialize MetricStats");
        assert_eq!(back, s);
    }

    #[test]
    fn resample_block_returns_n_values_from_contiguous_runs() {
        // The extracted kernel's standalone contract (now a public entry point):
        // given a pre-clamped block_len in [1, n], it yields exactly `n` values, each
        // appearing in `rs`, arranged as contiguous runs of `rs` (the final run
        // truncated to reach exactly n). Distinct powers of ten make membership a sum
        // check. n=5, block_len=2 -> runs of length 2, 2, 1.
        let rs = [1.0, 10.0, 100.0, 1000.0, 10000.0];
        let n = rs.len();
        let block_len = 2;
        let mut rng = SplitMix64::new(7);
        let sample = resample_block(&rs, block_len, &mut rng);
        assert_eq!(sample.len(), n, "resample has exactly n values");
        assert!(
            sample.iter().all(|v| rs.contains(v)),
            "every resampled value is drawn from rs",
        );
        // Reconstruct the contiguous-run structure: walk the sample in blocks of
        // `block_len` (last truncated) and assert each block is a contiguous slice of rs.
        let mut filled = 0usize;
        while filled < n {
            let take = block_len.min(n - filled);
            let run = &sample[filled..filled + take];
            let start = rs
                .iter()
                .position(|v| v == &run[0])
                .expect("run head is in rs");
            assert_eq!(
                run,
                &rs[start..start + take],
                "each run is a contiguous slice of rs",
            );
            filled += take;
        }
    }

    #[test]
    fn one_sided_p_laplace_known_values() {
        assert_eq!(one_sided_p_laplace(0, 3), 0.25);
        assert_eq!(one_sided_p_laplace(3, 3), 1.0);
        assert_eq!(one_sided_p_laplace(1, 999), 2.0 / 1000.0);
        // the registry's R-arm shape: count-of-ge over usable draws
        assert_eq!(one_sided_p_laplace(0, 999), 1.0 / 1000.0);
    }
}