diff --git a/crates/aura-engine/src/lib.rs b/crates/aura-engine/src/lib.rs index 0d7e313..34cafb7 100644 --- a/crates/aura-engine/src/lib.rs +++ b/crates/aura-engine/src/lib.rs @@ -37,6 +37,7 @@ mod blueprint; mod builder; mod graph_model; mod harness; +mod mc; mod report; mod sweep; @@ -51,6 +52,7 @@ pub use harness::{ }; pub use report::{f64_field, summarize, RunManifest, RunMetrics, RunReport}; pub use sweep::{sweep, GridSpace, SweepError, SweepFamily, SweepPoint}; +pub use mc::{monte_carlo, McAggregate, McDraw, McFamily, MetricStats}; // #29: re-export the core scalar vocabulary a Blueprint builder needs // (SourceSpec.kind is a ScalarKind; sources/Recorder columns are Scalar / // Firing / Timestamp) so a graph builder has one import surface, not two. diff --git a/crates/aura-engine/src/mc.rs b/crates/aura-engine/src/mc.rs new file mode 100644 index 0000000..73a07a9 --- /dev/null +++ b/crates/aura-engine/src/mc.rs @@ -0,0 +1,281 @@ +//! 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, + 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 = 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::() / 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(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( + 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 = 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::>(), 0); + let exposure = f64_field(&rx_ex.try_iter().collect::>(), 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 { + // 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 { + 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![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 = (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); + } +} diff --git a/crates/aura-engine/src/sweep.rs b/crates/aura-engine/src/sweep.rs index b384a65..5918312 100644 --- a/crates/aura-engine/src/sweep.rs +++ b/crates/aura-engine/src/sweep.rs @@ -121,34 +121,34 @@ where sweep_with_threads(space, nthreads, run_one) } -/// The thread-count-explicit core of [`sweep`]. Module-private: the public -/// `sweep` derives the count, while the tests drive it at 1 and at N to pin -/// determinism under concurrency (C1). Workers pull point indices from a shared -/// atomic cursor (work-stealing load-balances uneven per-point cost); each tags -/// its result with the point index, and the family is assembled by sorting on -/// that index after the scope joins — so the order is the enumeration order, not -/// the completion order. Only the cursor is shared; the results side is lock-free. -fn sweep_with_threads(space: &GridSpace, nthreads: usize, run_one: F) -> SweepFamily +/// Run `n` disjoint jobs in parallel and collect their results in **job-index +/// order** — the shared disjoint-parallel core both `sweep` (over grid points) +/// and `monte_carlo` (over seeds) drive (C1: order is the input order, not the +/// completion order). Workers pull job indices from a shared atomic cursor +/// (work-stealing load-balances uneven per-job cost); each tags its result with +/// the index, and the results are sorted on that index after the scope joins. +/// Only the cursor is shared; the results side is lock-free. `nthreads` is +/// clamped to `[1, n.max(1)]` (a 0-job call yields an empty vec; a 0 thread +/// count coerces to 1). +pub(crate) fn run_indexed(n: usize, nthreads: usize, run_one: F) -> Vec where - F: Fn(&[Scalar]) -> RunReport + Sync, + T: Send, + F: Fn(usize) -> T + Sync, { - let points = space.points(); - // `points.len() >= 1` always (GridSpace rejects empty axes), so the clamp - // lower bound never exceeds the upper; this also coerces a 0 to 1 worker. - let nthreads = nthreads.clamp(1, points.len()); + let nthreads = nthreads.clamp(1, n.max(1)); let cursor = AtomicUsize::new(0); - let mut results: Vec<(usize, RunReport)> = std::thread::scope(|scope| { + let mut results: Vec<(usize, T)> = std::thread::scope(|scope| { let handles: Vec<_> = (0..nthreads) .map(|_| { scope.spawn(|| { - let mut local: Vec<(usize, RunReport)> = Vec::new(); + let mut local: Vec<(usize, T)> = Vec::new(); loop { let i = cursor.fetch_add(1, Ordering::Relaxed); - if i >= points.len() { + if i >= n { break; } - local.push((i, run_one(&points[i]))); + local.push((i, run_one(i))); } local }) @@ -158,10 +158,25 @@ where }); results.sort_by_key(|&(i, _)| i); + results.into_iter().map(|(_, t)| t).collect() +} + +/// The thread-count-explicit core of [`sweep`]. Module-private: the public +/// `sweep` 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`]: it +/// enumerates the grid points, runs each disjointly, and zips the reports back +/// onto their points in enumeration (odometer) order. +fn sweep_with_threads(space: &GridSpace, nthreads: usize, run_one: F) -> SweepFamily +where + F: Fn(&[Scalar]) -> RunReport + Sync, +{ + let points = space.points(); + let reports = run_indexed(points.len(), nthreads, |i| run_one(&points[i])); SweepFamily { - points: results + points: points .into_iter() - .map(|(i, report)| SweepPoint { params: points[i].clone(), report }) + .zip(reports) + .map(|(params, report)| SweepPoint { params, report }) .collect(), } }