Aura/crates/aura-backtest/src/mc.rs

//! Monte-Carlo orchestration family (C12 axis 4): Monte-Carlo as a sweep over
//! seeds. `monte_carlo(base_point, seeds, run_one)` runs a fixed base point over
//! a seed set — each seed a disjoint C1 realization — and collects an
//! [`McFamily`]: the per-seed [`McDraw`]s in seed-input order plus a stored
//! [`McAggregate`] (mean + quantiles of all three run metrics across the
//! realizations). It reuses the disjoint-parallel [`run_indexed`](aura_engine::run_indexed)
//! core `sweep` drives — the varying dimension is the *seed*, not a tuning param.
//! Eager-agnostic (C12/#71): the API takes seeds + a per-draw closure, never a
//! materialized stream `Vec`; the seed -> `Source` construction is a closure-body
//! concern. Moved verbatim from `aura-engine::mc` (#291, C28 phase 2): the
//! engine stays metric-agnostic, so the concrete `RunMetrics` instantiation of
//! `RunReport` lives here, the outer C28 rung.

// MetricStats/resample_block/SplitMix64 are foundation-grade (aura-analysis),
// but reached here through aura-engine's re-export rather than a direct
// aura-analysis dependency — the C28 ladder (c28_layering) permits
// aura-backtest -> {aura-core, aura-engine} only, not a same-rung-skipping
// direct edge to the foundation crate.
use aura_core::Scalar;
use aura_engine::{run_indexed, resample_block, MetricStats, SplitMix64};
use crate::{RunMetrics, RunReport};

/// One Monte-Carlo realization: the seed that drove it and the full `RunReport`.
/// Self-describing, analog to [`SweepPoint`](aura_engine::SweepPoint).
#[derive(Clone, Debug, PartialEq)]
pub struct McDraw {
    pub seed: u64,
    pub report: RunReport,
}

/// The result family of a Monte-Carlo run — one [`McDraw`] per seed, in
/// seed-**input** order (independent of thread completion), plus the stored
/// [`McAggregate`]. Analog to [`SweepFamily`](aura_engine::SweepFamily).
#[derive(Clone, Debug, PartialEq)]
pub struct McFamily {
    pub draws: Vec<McDraw>,
    pub aggregate: McAggregate,
}

/// Distribution summary of all three run metrics across the realizations: covers
/// every metric (not a single "chosen" one). A pure post-run reduction over the
/// draws — stored for the common robustness case; custom statistics read the raw
/// draws directly.
#[derive(Clone, Debug, PartialEq)]
pub struct McAggregate {
    pub total_pips: MetricStats,
    pub max_drawdown: MetricStats,
    pub bias_sign_flips: MetricStats,
}

impl McAggregate {
    /// Pure reduction over the realizations — recomputable from `draws` alone (it
    /// is exactly what [`monte_carlo`] stored). `draws` must be non-empty (a
    /// Monte-Carlo over zero realizations has no defined mean/quantile); on a
    /// non-empty slice every field is finite.
    pub fn from_draws(draws: &[McDraw]) -> McAggregate {
        let pick = |f: fn(&RunMetrics) -> f64| -> MetricStats {
            let xs: Vec<f64> = draws.iter().map(|d| f(&d.report.metrics)).collect();
            MetricStats::from_values(&xs)
        };
        McAggregate {
            total_pips: pick(|m| m.total_pips),
            max_drawdown: pick(|m| m.max_drawdown),
            bias_sign_flips: pick(|m| m.bias_sign_flips as f64),
        }
    }
}

/// The single assembly path behind both [`monte_carlo`] and the thread-count-
/// explicit test wrapper: run `run_one` over every seed via `run_indexed` (on
/// whichever pool is ambient when called), build the per-seed draws in
/// seed-input order, and reduce them into the aggregate. One source of truth so
/// the two callers can never drift on assembly logic while their pool selection
/// differs.
fn assemble_mc<F>(base_point: &[Scalar], seeds: &[u64], run_one: &F) -> McFamily
where
    F: Fn(u64, &[Scalar]) -> RunReport + Sync,
{
    let reports = run_indexed(seeds.len(), |i| run_one(seeds[i], base_point));
    let draws: Vec<McDraw> = seeds
        .iter()
        .zip(reports)
        .map(|(&seed, report)| McDraw { seed, report })
        .collect();
    let aggregate = McAggregate::from_draws(&draws);
    McFamily { draws, aggregate }
}

/// Run `run_one(seed, base_point)` over every seed, disjointly in parallel (C1),
/// and collect the family in seed-input order plus the aggregate. The varying
/// dimension is the *seed* (C12 axis 4: MC = sweep over seeds), not a tuning
/// param; `base_point` is constant across draws. Eager-agnostic: the seed ->
/// `Source` construction lives inside `run_one`, never a materialized stream
/// `Vec` in this API (#71).
///
/// Precondition: `seeds` is non-empty (a Monte-Carlo over zero realizations has
/// no defined aggregate). A `debug_assert!` guards it; the `-> McFamily`
/// signature is preserved (no `Result`), matching [`sweep`](aura_engine::sweep).
pub fn monte_carlo<F>(base_point: &[Scalar], seeds: &[u64], run_one: F) -> McFamily
where
    F: Fn(u64, &[Scalar]) -> RunReport + Sync,
{
    debug_assert!(!seeds.is_empty(), "monte_carlo requires a non-empty seed set");
    assemble_mc(base_point, seeds, &run_one)
}

/// The thread-count-explicit wrapper of [`monte_carlo`]. Module-private: the
/// public `monte_carlo` runs on the ambient pool; the tests drive this at 1 and
/// at N to pin determinism (C1). `assemble_mc` (run_indexed + draw assembly +
/// aggregate) runs inside a local rayon pool of `nthreads` workers, sharing the
/// exact assembly path `monte_carlo` uses — no second copy to drift. Bounds
/// match `monte_carlo` (`F: Sync`): the install closure only borrows `run_one`,
/// never moves it.
#[cfg(test)]
fn monte_carlo_with_threads<F>(
    base_point: &[Scalar],
    seeds: &[u64],
    nthreads: usize,
    run_one: F,
) -> McFamily
where
    F: Fn(u64, &[Scalar]) -> RunReport + Sync,
{
    debug_assert!(!seeds.is_empty(), "monte_carlo requires a non-empty seed set");
    let pool = rayon::ThreadPoolBuilder::new()
        .num_threads(nthreads)
        .build()
        .expect("rayon thread pool");
    pool.install(|| assemble_mc(base_point, seeds, &run_one))
}

/// Distribution of `E[R]` under a moving-block bootstrap of an OOS per-trade R series
/// (#139). `block_len == 1` is the i.i.d. trade-shuffle; `block_len > 1` resamples
/// contiguous runs, preserving the serial correlation of sequential trades a pure
/// shuffle would erase. Deterministic (C1) given `seed`.
#[derive(Clone, Debug, PartialEq, serde::Serialize, serde::Deserialize)]
pub struct RBootstrap {
    pub e_r: MetricStats,   // mean + p5/p25/p50/p75/p95 of the resampled E[R]
    pub prob_le_zero: f64,  // fraction of resamples whose mean R <= 0
    pub n_trades: usize,
    pub block_len: usize,
    pub n_resamples: usize,
}

/// Moving-block bootstrap of `rs` (non-circular; the final block of each resample is
/// truncated so the resample has exactly `n` values). `block_len` is clamped to
/// `[1, n]`. Empty `rs` or zero resamples -> an all-zero `RBootstrap`. Pure given
/// `seed` (drives the existing `SplitMix64`).
pub fn r_bootstrap(rs: &[f64], n_resamples: usize, block_len: usize, seed: u64) -> RBootstrap {
    let n = rs.len();
    if n == 0 || n_resamples == 0 {
        return RBootstrap {
            e_r: MetricStats { mean: 0.0, p5: 0.0, p25: 0.0, p50: 0.0, p75: 0.0, p95: 0.0 },
            prob_le_zero: 0.0,
            n_trades: n,
            block_len: block_len.clamp(1, n.max(1)),
            n_resamples,
        };
    }
    let block_len = block_len.clamp(1, n);
    let mut rng = SplitMix64::new(seed);
    let mut means: Vec<f64> = Vec::with_capacity(n_resamples);
    for _ in 0..n_resamples {
        let sample = resample_block(rs, block_len, &mut rng);
        means.push(sample.iter().sum::<f64>() / n as f64);
    }
    let e_r = MetricStats::from_values(&means);
    let prob_le_zero = means.iter().filter(|&&m| m <= 0.0).count() as f64 / n_resamples as f64;
    RBootstrap { e_r, prob_le_zero, n_trades: n, block_len, n_resamples }
}

#[cfg(test)]
mod tests {
    use super::*;
    use aura_engine::{
        f64_field, BlueprintNode, Composite, Edge, OutField, Role, RunManifest, SyntheticSpec,
        Target, Timestamp,
    };
    use aura_core::{Firing, Scalar, ScalarKind};
    use aura_std::{Recorder, Sma, Sub};
    use aura_strategy::Bias;
    use crate::{summarize, SimBroker};
    use std::sync::mpsc;

    // The SMA-cross signal-quality harness fixture, reconstructed here through
    // the PUBLIC aura-engine/aura-std/aura-strategy/aura-backtest surface —
    // `aura-engine::test_fixtures` is `#[cfg(test)]`-private to the engine crate
    // itself and unreachable from here (the same reason `random_sweep_e2e.rs` /
    // `list_sweep_e2e.rs` carry their own copy). Bodies verbatim from
    // `aura-engine/src/test_fixtures.rs`'s `sma_cross`/`composite_sma_cross_harness`.

    /// The SMA-cross signal as a reusable value-empty composite: one input role
    /// (price), one output (fast-minus-slow spread). Lengths injected at compile.
    fn sma_cross() -> Composite {
        Composite::new(
            "sma_cross",
            vec![
                Sma::builder().named("fast").into(),
                Sma::builder().named("slow").into(),
                Sub::builder().into(),
            ],
            vec![
                Edge { from: 0, to: 2, slot: 0, from_field: 0 },
                Edge { from: 1, to: 2, slot: 1, from_field: 0 },
            ],
            vec![Role {
                name: "price".into(),
                targets: vec![Target { node: 0, slot: 0 }, Target { node: 1, slot: 0 }],
                source: None,
            }],
            vec![OutField { node: 2, field: 0, name: "out".into() }],
        )
    }

    /// The signal-quality harness as a value-empty composite blueprint with two
    /// recording sinks (equity, exposure).
    #[allow(clippy::type_complexity)]
    fn composite_sma_cross_harness() -> (
        Composite,
        mpsc::Receiver<(Timestamp, Vec<Scalar>)>,
        mpsc::Receiver<(Timestamp, Vec<Scalar>)>,
    ) {
        let (tx_eq, rx_eq) = mpsc::channel();
        let (tx_ex, rx_ex) = mpsc::channel();
        let bp = Composite::new(
            "root",
            vec![
                BlueprintNode::Composite(sma_cross()),
                Bias::builder().named("bias").into(),
                SimBroker::builder(0.0001).into(),
                Recorder::builder(vec![ScalarKind::F64], Firing::Any, tx_eq).into(),
                Recorder::builder(vec![ScalarKind::F64], Firing::Any, tx_ex).into(),
            ],
            vec![
                Edge { from: 0, to: 1, slot: 0, from_field: 0 }, // composite out -> Bias
                Edge { from: 1, to: 2, slot: 0, from_field: 0 }, // exposure -> broker slot 0
                Edge { from: 2, to: 3, slot: 0, from_field: 0 }, // equity -> sink
                Edge { from: 1, to: 4, slot: 0, from_field: 0 }, // exposure -> sink
            ],
            vec![Role {
                name: "src".into(),
                targets: vec![
                    Target { node: 0, slot: 0 }, // price -> sma_cross role 0
                    Target { node: 2, slot: 1 }, // price -> SimBroker price slot
                ],
                source: Some(ScalarKind::F64),
            }],
            vec![], // output: the root ends in sinks
        );
        (bp, rx_eq, rx_ex)
    }

    /// Build + bootstrap + run + drain + summarize one (seed, point) into a
    /// `RunReport`. A free `fn` (Copy + Sync) so it serves as the `monte_carlo`
    /// closure AND a direct reference for the "draw == independent run"
    /// comparison. The stream is generated from `seed` (seed-as-input, #66), and
    /// `seed` is recorded into the manifest.
    fn run_draw(seed: u64, point: &[Scalar]) -> RunReport {
        let (bp, rx_eq, rx_ex) = composite_sma_cross_harness();
        let mut h = bp
            .bootstrap_with_params(point.to_vec())
            .expect("base point is kind-checked against param_space");
        let spec = SyntheticSpec { start: 1.0, len: 64, step: 1 };
        let window = (Timestamp(1), Timestamp((spec.len as i64 - 1) * spec.step + 1));
        h.run(vec![Box::new(spec.source(seed))]);
        let equity = f64_field(&rx_eq.try_iter().collect::<Vec<_>>(), 0);
        let exposure = f64_field(&rx_ex.try_iter().collect::<Vec<_>>(), 0);
        RunReport {
            manifest: RunManifest {
                commit: "test".to_string(),
                params: Vec::new(),
                defaults: Vec::new(),
                window,
                seed,
                broker: "test".to_string(),
                selection: None,
                instrument: None,
                topology_hash: None,
                project: None,
            },
            metrics: summarize(&equity, &exposure),
        }
    }

    fn base_point() -> Vec<Scalar> {
        // matches the composite_sma_cross param_space order: fast, slow, scale
        vec![Scalar::i64(2), Scalar::i64(4), Scalar::f64(0.5)]
    }

    fn draws_with_metrics(values: &[f64]) -> Vec<McDraw> {
        values
            .iter()
            .enumerate()
            .map(|(i, &v)| McDraw {
                seed: i as u64,
                report: RunReport {
                    manifest: RunManifest {
                        commit: "t".to_string(),
                        params: Vec::new(),
                        defaults: Vec::new(),
                        window: (Timestamp(0), Timestamp(0)),
                        seed: i as u64,
                        broker: "t".to_string(),
                        selection: None,
                        instrument: None,
                        topology_hash: None,
                        project: None,
                    },
                    metrics: RunMetrics {
                        total_pips: v,
                        max_drawdown: v,
                        bias_sign_flips: v as u64,
                        r: None,
                    },
                },
            })
            .collect()
    }

    #[test]
    fn monte_carlo_runs_one_draw_per_seed_in_input_order() {
        // N seeds -> N draws, seeds carried in INPUT order.
        let family = monte_carlo(&base_point(), &[1, 2, 3], run_draw);
        assert_eq!(family.draws.len(), 3);
        assert_eq!(
            family.draws.iter().map(|d| d.seed).collect::<Vec<_>>(),
            vec![1, 2, 3],
        );
    }

    #[test]
    fn family_is_deterministic_across_thread_counts() {
        // order = seed input, not completion (C1).
        let point = base_point();
        let seeds: Vec<u64> = (1..=8).collect();
        let one = monte_carlo_with_threads(&point, &seeds, 1, run_draw);
        let many = monte_carlo_with_threads(&point, &seeds, 8, run_draw);
        assert_eq!(one, many);
        assert_eq!(one, monte_carlo(&point, &seeds, run_draw));
    }

    #[test]
    fn draw_equals_independent_seeded_run() {
        // each draw == an independent run of that (seed, point)
        // — MC adds enumeration + execution, never a metrics change (C1).
        let point = base_point();
        let family = monte_carlo(&point, &[10, 11, 12], run_draw);
        for draw in &family.draws {
            assert_eq!(draw.report, run_draw(draw.seed, &point));
            assert!(draw.report.metrics.total_pips.is_finite());
        }
    }

    #[test]
    fn distinct_seeds_produce_distinct_draws() {
        // the seed actually perturbs each run — the family does
        // not collapse to one metric.
        let family = monte_carlo(&base_point(), &[1, 2, 3, 4, 5], run_draw);
        let first = family.draws[0].report.metrics.total_pips;
        assert!(
            family.draws.iter().any(|d| d.report.metrics.total_pips != first),
            "a multi-seed family must not collapse to one metric",
        );
    }

    #[test]
    fn aggregate_is_a_pure_reduction() {
        // the stored aggregate is recomputable from the draws.
        let family = monte_carlo(&base_point(), &[1, 2, 3, 4], run_draw);
        assert_eq!(McAggregate::from_draws(&family.draws), family.aggregate);
    }

    #[test]
    fn aggregate_stats_on_known_fixture() {
        // type-7 quantile + mean over known metric values
        // [0,1,2,3,4]: mean=2.0, p50=2.0, p5≈0.2, p95≈3.8.
        let agg = McAggregate::from_draws(&draws_with_metrics(&[0.0, 1.0, 2.0, 3.0, 4.0]));
        assert_eq!(agg.total_pips.mean, 2.0);
        assert_eq!(agg.total_pips.p50, 2.0);
        assert!((agg.total_pips.p5 - 0.2).abs() < 1e-9, "p5 = {}", agg.total_pips.p5);
        assert!((agg.total_pips.p95 - 3.8).abs() < 1e-9, "p95 = {}", agg.total_pips.p95);
    }

    #[test]
    fn r_bootstrap_empty_series_is_all_zero() {
        let b = r_bootstrap(&[], 100, 1, 7);
        assert_eq!(b.n_trades, 0);
        assert_eq!(b.e_r.mean, 0.0);
        assert_eq!(b.prob_le_zero, 0.0);
    }

    #[test]
    fn r_bootstrap_single_block_equals_full_series_mean() {
        // block_len == n: the only valid start is 0, so every resample IS the full
        // series -> every resample mean == the series mean -> zero spread.
        let rs = [1.0, -2.0, 3.0]; // mean = 2/3
        let b = r_bootstrap(&rs, 500, rs.len(), 7);
        let mean = 2.0 / 3.0;
        assert!((b.e_r.mean - mean).abs() < 1e-12);
        assert!((b.e_r.p5 - mean).abs() < 1e-12);
        assert!((b.e_r.p95 - mean).abs() < 1e-12);
        assert_eq!(b.prob_le_zero, 0.0); // mean > 0
        assert_eq!(b.n_trades, 3);
        assert_eq!(b.block_len, 3);
        assert_eq!(b.n_resamples, 500);
    }

    #[test]
    fn r_bootstrap_is_deterministic_given_seed() {
        let rs = [0.5, -1.0, 2.0, -0.5, 1.5, -2.0, 0.25];
        let a = r_bootstrap(&rs, 1000, 1, 42);
        let b = r_bootstrap(&rs, 1000, 1, 42);
        assert_eq!(a, b, "same rs+seed+params -> identical RBootstrap (C1)");
        let c = r_bootstrap(&rs, 1000, 1, 43);
        assert_ne!(a.e_r.p5, c.e_r.p5, "a different seed reshuffles differently");
    }

    #[test]
    fn r_bootstrap_block_len_is_clamped_to_series_len() {
        // block_len > n clamps to n -> single-block behaviour (no panic, no OOB).
        let rs = [1.0, 2.0];
        let b = r_bootstrap(&rs, 10, 99, 1);
        assert_eq!(b.block_len, 2);
        assert!((b.e_r.mean - 1.5).abs() < 1e-12);
    }

    #[test]
    fn r_bootstrap_moving_block_resamples_are_contiguous_truncated_runs() {
        // The serial-correlation-preserving branch: 1 < block_len < n with n NOT a
        // multiple of block_len, so the final block of every resample is truncated
        // (`take = block_len.min(n - sample.len())`). Every resample is a sequence of
        // contiguous runs of `rs` (the last one cut short to reach exactly n values),
        // so each resample mean must lie in the finite set of legal block-composition
        // means. Distinct powers of ten make every distinct multiset of contributions
        // a distinct sum -> a faithful membership check. n=5, block_len=2 -> blocks 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;

        // Enumerate every legal resample sum: pick a start in [0, n-block_len] for each
        // of the three blocks (lengths 2, 2, 1), summing min(block_len, n-filled) values
        // from that start. Membership-by-sum, mirroring r_bootstrap's own construction.
        let starts: Vec<usize> = (0..=n - block_len).collect();
        let mut legal_means: Vec<f64> = Vec::new();
        for &s0 in &starts {
            for &s1 in &starts {
                for &s2 in &starts {
                    let mut sum = 0.0;
                    let mut filled = 0usize;
                    for &start in &[s0, s1, s2] {
                        let take = block_len.min(n - filled);
                        if take == 0 {
                            break;
                        }
                        sum += rs[start..start + take].iter().sum::<f64>();
                        filled += take;
                    }
                    legal_means.push(sum / n as f64);
                }
            }
        }

        // Drive enough resamples that the truncated final block is exercised many
        // times; assert every produced mean is a legal block-composition mean.
        let b = r_bootstrap(&rs, 2000, block_len, 7);
        assert_eq!(b.block_len, 2);
        assert_eq!(b.n_trades, 5);
        // The mean of resample means is a weighted average of legal means, so it is
        // bounded by their min and max (not itself a single legal mean).
        let lo = legal_means.iter().cloned().fold(f64::INFINITY, f64::min);
        let hi = legal_means.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
        assert!(lo - 1e-6 <= b.e_r.mean && b.e_r.mean <= hi + 1e-6);
        // Each quantile of the resampled E[R] is an order statistic of the resample
        // means, hence itself a legal contiguous-block-composition mean.
        for &m in &[b.e_r.p5, b.e_r.p25, b.e_r.p50, b.e_r.p75, b.e_r.p95] {
            assert!(
                legal_means.iter().any(|&lm| (lm - m).abs() < 1e-6),
                "resampled quantile {m} is not a legal contiguous-block-composition mean",
            );
        }
    }
}