Aura/crates/aura-engine/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`](crate::sweep)
//! 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.

use crate::sweep::run_indexed;
use crate::{RunMetrics, RunReport, Scalar};

/// One Monte-Carlo realization: the seed that drove it and the full `RunReport`.
/// Self-describing, analog to [`SweepPoint`](crate::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`](crate::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 exposure_sign_flips: MetricStats,
}

/// 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)]
pub struct MetricStats {
    pub mean: f64,
    pub p5: f64,
    pub p25: f64,
    pub p50: f64, // == median
    pub p75: f64,
    pub p95: f64,
}

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 mut xs: Vec<f64> = draws.iter().map(|d| f(&d.report.metrics)).collect();
            xs.sort_by(|a, b| a.partial_cmp(b).expect("run metrics 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),
            }
        };
        McAggregate {
            total_pips: pick(|m| m.total_pips),
            max_drawdown: pick(|m| m.max_drawdown),
            exposure_sign_flips: pick(|m| m.exposure_sign_flips as f64),
        }
    }
}

/// 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.
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]
    }
}

/// 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`](crate::sweep).
pub fn monte_carlo<F>(base_point: &[Scalar], seeds: &[u64], run_one: F) -> McFamily
where
    F: Fn(u64, &[Scalar]) -> RunReport + Sync,
{
    let nthreads = std::thread::available_parallelism().map(|n| n.get()).unwrap_or(1);
    monte_carlo_with_threads(base_point, seeds, nthreads, run_one)
}

/// The thread-count-explicit core of [`monte_carlo`]. Module-private: the public
/// `monte_carlo` derives the count, while the tests drive it at 1 and at N to pin
/// determinism under concurrency (C1). A thin adapter over
/// [`run_indexed`](crate::sweep): it runs each seed disjointly and zips the
/// reports back onto their seeds in seed-input order.
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 reports = run_indexed(seeds.len(), nthreads, |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 }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::test_fixtures::composite_sma_cross_harness;
    use crate::{f64_field, summarize, RunManifest, SyntheticSpec, Timestamp};
    use aura_core::Scalar;

    /// 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(),
                window,
                seed,
                broker: "test".to_string(),
            },
            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(),
                        window: (Timestamp(0), Timestamp(0)),
                        seed: i as u64,
                        broker: "t".to_string(),
                    },
                    metrics: RunMetrics {
                        total_pips: v,
                        max_drawdown: v,
                        exposure_sign_flips: v as u64,
                    },
                },
            })
            .collect()
    }

    #[test]
    fn monte_carlo_runs_one_draw_per_seed_in_input_order() {
        // spec §Testing 1: 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() {
        // spec §Testing 2: 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() {
        // spec §Testing 3: 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() {
        // spec §Testing 4: 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() {
        // spec §Testing 5: 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() {
        // spec §Testing 6: 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 quantile_endpoints_and_singleton() {
        // spec §Testing 7: 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);
    }
}