diff --git a/crates/ml-alpha/examples/alpha_bar_baseline.rs b/crates/ml-alpha/examples/alpha_bar_baseline.rs deleted file mode 100644 index 61254770e..000000000 --- a/crates/ml-alpha/examples/alpha_bar_baseline.rs +++ /dev/null @@ -1,129 +0,0 @@ -//! Phase 1a entry-point binary. -//! -//! Loads local fxcache, runs purged walk-forward split, trains a 2-layer MLP -//! on binary direction labels, evaluates validation accuracy, prints the gate -//! verdict. -//! -//! ## Usage -//! -//! Defaults to the local Phase 1a test fxcache: -//! ``` -//! cargo run -p ml-alpha --example phase1a --release -//! ``` -//! -//! Override config via CLI: -//! ``` -//! cargo run -p ml-alpha --example phase1a --release -- \ -//! --fxcache-path /path/to/file.fxcache \ -//! --horizon 60 \ -//! --epochs 5 -//! ``` -//! -//! Or via TOML config: -//! ``` -//! cargo run -p ml-alpha --example phase1a --release -- --config config/foxhuntq-phase1a.toml -//! ``` - -use std::sync::Arc; - -use anyhow::Result; -use clap::Parser; -use cudarc::driver::CudaContext; -use tracing_subscriber::EnvFilter; - -use ml_alpha::training::{Phase1aConfig, Phase1aTrainer}; - -#[derive(Debug, Parser)] -#[command(name = "alpha_bar_baseline", about = "FoxhuntQ-Δ Phase 1a falsification smoke")] -struct Cli { - /// Override fxcache path. Defaults to local test_data fxcache. - #[arg(long)] - fxcache_path: Option, - - /// Number of training epochs. - #[arg(long, default_value_t = 5)] - epochs: usize, - - /// Batch size for training. - #[arg(long, default_value_t = 1024)] - batch_size: usize, - - /// Adam learning rate. - #[arg(long, default_value_t = 1e-3)] - learning_rate: f32, - - /// Label horizon in bars (`sign(price[t+H] − price[t])`). - #[arg(long, default_value_t = 60)] - horizon: usize, - - /// Train fraction for the purged walk-forward split. - #[arg(long, default_value_t = 0.8)] - train_frac: f32, - - /// Additional embargo bars beyond `horizon` (Lopez de Prado embargo). - #[arg(long, default_value_t = 0)] - embargo_bars: usize, - - /// Hidden dim of the 2-layer MLP. - #[arg(long, default_value_t = 256)] - hidden_dim: usize, - - /// Deterministic RNG seed. - #[arg(long, default_value_t = 42)] - seed: u64, - - /// Path to a TOML config file (overrides individual flags). - #[arg(long)] - config: Option, -} - -fn main() -> Result<()> { - tracing_subscriber::fmt() - .with_env_filter(EnvFilter::try_from_default_env().unwrap_or_else(|_| EnvFilter::new("info"))) - .init(); - - let cli = Cli::parse(); - - if cli.config.is_some() { - // TOML config support deferred to step 3; skeleton uses CLI defaults only. - anyhow::bail!("--config (TOML) not yet supported; use CLI flags"); - } - let mut config = Phase1aConfig::default(); - - // Apply CLI overrides over the loaded / default config. - if let Some(p) = cli.fxcache_path { config.fxcache_path = p; } - config.epochs = cli.epochs; - config.batch_size = cli.batch_size; - config.learning_rate = cli.learning_rate; - config.horizon = cli.horizon; - config.train_frac = cli.train_frac; - config.embargo_bars = cli.embargo_bars; - config.mlp.hidden_dim = cli.hidden_dim; - config.seed = cli.seed; - - tracing::info!(?config, "Phase 1a config resolved"); - - // Initialize CUDA context + stream. RTX 3050 Ti on local; H100/L40S on - // Argo (compute-cap auto-derived by the workflow). - let ctx = CudaContext::new(0) - .map_err(|e| anyhow::anyhow!("CUDA context init: {e}"))?; - let stream = ctx.new_stream() - .map_err(|e| anyhow::anyhow!("CUDA stream init: {e}"))?; - let stream = Arc::new(stream); - - let mut trainer = Phase1aTrainer::from_config(config, Arc::clone(&stream))?; - let report = trainer.run()?; - - tracing::info!( - accuracy = report.accuracy, - auc = report.auc, - n_samples = report.n_samples, - up_fraction = report.up_fraction, - tied_fraction = report.tied_fraction, - "Phase 1a evaluation complete" - ); - println!(); - println!("{}", report.verdict_line()); - - Ok(()) -} diff --git a/crates/ml-alpha/examples/alpha_bar_detailed.rs b/crates/ml-alpha/examples/alpha_bar_detailed.rs deleted file mode 100644 index 0f96d2f42..000000000 --- a/crates/ml-alpha/examples/alpha_bar_detailed.rs +++ /dev/null @@ -1,165 +0,0 @@ -//! Phase 1c detailed metrics smoke — calibration + stratified accuracy. -//! -//! Trains the same Phase 1a MLP as `phase1a`, then layers on: -//! - Brier score (strictly proper scoring rule, lower = better) -//! - Log-loss (cross-entropy, lower = better) -//! - 10-bin reliability curve (predicted-prob vs observed-positive-rate) -//! - Quintile-stratified accuracy on each of the 6 Block-S snapshot features -//! -//! The intent: diagnose the AUC≫accuracy gap we saw at snapshot resolution. -//! If accuracy lifts dramatically within specific feature quintiles, alpha -//! is regime-conditional (concentrated in high-event-rate / wide-spread / -//! specific-time-of-day bars). If reliability is monotonic but compressed -//! around 0.5, the gap is plain miscalibration — fixable with a sigmoid -//! threshold tune at deployment. - -use anyhow::{Context, Result}; -use clap::Parser; -use tracing_subscriber::EnvFilter; - -use cudarc::driver::CudaContext; - -use ml_alpha::fxcache_reader::FxCacheReader; -use ml_alpha::metrics_detail::{ - brier_score, log_loss, reliability_curve, stratified_accuracy, -}; -use ml_alpha::training::{Phase1aConfig, Phase1aTrainer}; - -#[derive(Parser, Debug)] -#[command(name = "alpha_bar_detailed", about = "Phase 1c detailed metrics (calibration + stratification)")] -struct Cli { - /// Path to fxcache file (alpha feature column required). - #[arg(long)] - fxcache_path: String, - - /// Number of training epochs. - #[arg(long, default_value_t = 5)] - epochs: usize, - - /// Forward-horizon for the binary direction label (in rows). - #[arg(long, default_value_t = 100)] - horizon: usize, - - /// Number of bins for the reliability curve. - #[arg(long, default_value_t = 10)] - n_reliability_bins: usize, - - /// Number of equi-frequency strata for the stratified-accuracy analysis. - #[arg(long, default_value_t = 5)] - n_strata: usize, -} - -/// Names of the 6 Block-S features (positions 75..81 of the 81-dim snapshot row). -const BLOCK_S_FEATURE_NAMES: [&str; 6] = [ - "time_since_trade_s", - "time_since_snap_s", - "book_event_rate_per_s", - "spread_bps", - "L1_imbalance", - "micro_mid_drift", -]; -const BLOCK_S_OFFSET: usize = 75; - -fn main() -> Result<()> { - tracing_subscriber::fmt() - .with_env_filter(EnvFilter::try_from_default_env().unwrap_or_else(|_| EnvFilter::new("info"))) - .init(); - - let cli = Cli::parse(); - let ctx = CudaContext::new(0).context("init CUDA context (GPU 0)")?; - let stream = ctx.default_stream(); - - let mut config = Phase1aConfig::default(); - config.fxcache_path = cli.fxcache_path.clone(); - config.epochs = cli.epochs; - config.horizon = cli.horizon; - - let mut trainer = Phase1aTrainer::from_config(config, stream)?; - let outputs = trainer.run_full().context("train + eval")?; - let report = &outputs.report; - - println!(); - println!("================================================================================"); - println!("PHASE 1c DETAILED METRICS"); - println!("================================================================================"); - println!("Headline: accuracy={:.4} AUC={:.4} n_val={}", report.accuracy, report.auc, report.n_samples); - println!("Horizon: {} rows forward", cli.horizon); - println!("Up fraction: {:.4} (val)", report.up_fraction); - println!(); - - // ── Calibration ──────────────────────────────────────────────────── - let brier = brier_score(&outputs.val_logits, &outputs.val_labels); - let logl = log_loss(&outputs.val_logits, &outputs.val_labels); - let chance_brier = report.up_fraction * (1.0 - report.up_fraction); // optimal for prior-only baseline - let chance_logl = -(report.up_fraction.ln() * report.up_fraction - + (1.0 - report.up_fraction).ln() * (1.0 - report.up_fraction)); - println!("--- Calibration ---"); - println!("Brier score: {:.5} (chance baseline = {:.5})", brier, chance_brier); - println!("Log loss: {:.5} (chance baseline = {:.5})", logl, chance_logl); - println!(); - - println!("--- Reliability curve ({} bins) ---", cli.n_reliability_bins); - println!("{:>7} {:>7} {:>10} {:>11} {:>14}", "bin_lo", "bin_hi", "n", "mean_pred", "observed_pos"); - let rel = reliability_curve(&outputs.val_logits, &outputs.val_labels, cli.n_reliability_bins); - for b in &rel { - println!( - "{:>7.3} {:>7.3} {:>10} {:>11.4} {:>14.4}", - b.bin_lo, b.bin_hi, b.n, b.mean_pred, b.observed_pos_rate - ); - } - println!(); - - // ── Stratification ───────────────────────────────────────────────── - // Pull raw (un-normalized) Block-S feature values via fxcache re-open. - let reader = FxCacheReader::open(&cli.fxcache_path) - .with_context(|| format!("reopen fxcache for stratification: {}", &cli.fxcache_path))?; - let alpha_dim = reader - .alpha_feature_dim() - .ok_or_else(|| anyhow::anyhow!("fxcache has no alpha column — stratification needs Block S"))?; - if alpha_dim < BLOCK_S_OFFSET + BLOCK_S_FEATURE_NAMES.len() { - anyhow::bail!( - "fxcache alpha_dim ({}) too small for Block S (need ≥ {})", - alpha_dim, - BLOCK_S_OFFSET + BLOCK_S_FEATURE_NAMES.len() - ); - } - - // Collect per-val-sample feature values for each Block-S column. - let n_val = outputs.val_indices.len(); - let mut col_vals: Vec> = (0..BLOCK_S_FEATURE_NAMES.len()) - .map(|_| Vec::with_capacity(n_val)) - .collect(); - for &bar_idx in &outputs.val_indices { - let row = reader - .alpha_features(bar_idx) - .ok_or_else(|| anyhow::anyhow!("missing alpha row at bar {}", bar_idx))?; - for (col_off, dst) in col_vals.iter_mut().enumerate() { - dst.push(row[BLOCK_S_OFFSET + col_off]); - } - } - - println!("--- Stratified accuracy (Block-S features, {} quantile bins) ---", cli.n_strata); - for (col_off, dst) in col_vals.iter().enumerate() { - let name = BLOCK_S_FEATURE_NAMES[col_off]; - let strats = stratified_accuracy( - &outputs.val_logits, - &outputs.val_labels, - dst, - cli.n_strata, - ); - println!(); - println!(" ▸ feature: {}", name); - println!( - " {:>4} {:>12} {:>12} {:>10} {:>10} {:>12}", - "k", "feat_lo", "feat_hi", "n", "accuracy", "observed" - ); - for s in &strats { - println!( - " {:>4} {:>12.4} {:>12.4} {:>10} {:>10.4} {:>12.4}", - s.stratum, s.feature_lo, s.feature_hi, s.n, s.accuracy, s.observed_pos_rate - ); - } - } - println!(); - Ok(()) -} diff --git a/crates/ml-alpha/examples/alpha_calibrate.rs b/crates/ml-alpha/examples/alpha_calibrate.rs deleted file mode 100644 index dd50ba4a1..000000000 --- a/crates/ml-alpha/examples/alpha_calibrate.rs +++ /dev/null @@ -1,103 +0,0 @@ -//! Phase 1d.0 — calibration smoke. -//! -//! Train the snapshot-level MLP exactly as `alpha_bar_detailed` does, then: -//! 1. Split val 50/50 into (calibration set, held-out test set). -//! 2. Fit Platt and isotonic on the calibration set. -//! 3. Apply each to the held-out test set; report Brier + log-loss for -//! uncalibrated / Platt / isotonic. -//! 4. Gate decision: if min(Platt_Brier, isotonic_Brier) ≤ chance baseline -//! `up_fraction × (1 - up_fraction)`, proceed to 1d.1; else falsify. - -use anyhow::{Context, Result}; -use clap::Parser; -use cudarc::driver::CudaContext; -use tracing_subscriber::EnvFilter; - -use ml_alpha::calibration::{Calibrator, IsotonicCalibrator, PlattScaler}; -use ml_alpha::metrics_detail::{brier_score, log_loss}; -use ml_alpha::training::{Phase1aConfig, Phase1aTrainer}; - -#[derive(Parser, Debug)] -#[command(name = "alpha_calibrate", about = "FoxhuntQ-Δ Phase 1d.0 calibration smoke")] -struct Cli { - #[arg(long)] - fxcache_path: String, - #[arg(long, default_value_t = 100)] - horizon: usize, - #[arg(long, default_value_t = 5)] - epochs: usize, - #[arg(long, default_value_t = 0.5)] - cal_split_frac: f32, -} - -/// Convert calibrated probabilities back to logits so we can reuse the -/// existing `brier_score` / `log_loss` helpers (both apply sigmoid internally). -fn probs_to_logits(probs: &[f32]) -> Vec { - probs - .iter() - .map(|&p| { - let p_clip = p.clamp(1e-7, 1.0 - 1e-7); - (p_clip / (1.0 - p_clip)).ln() - }) - .collect() -} - -fn main() -> Result<()> { - tracing_subscriber::fmt() - .with_env_filter(EnvFilter::try_from_default_env().unwrap_or_else(|_| EnvFilter::new("info"))) - .init(); - let cli = Cli::parse(); - let ctx = CudaContext::new(0).context("init CUDA")?; - let stream = ctx.default_stream(); - let mut config = Phase1aConfig::default(); - config.fxcache_path = cli.fxcache_path.clone(); - config.horizon = cli.horizon; - config.epochs = cli.epochs; - let mut trainer = Phase1aTrainer::from_config(config, stream)?; - let out = trainer.run_full()?; - - let n_val = out.val_logits.len(); - let n_cal = (n_val as f32 * cli.cal_split_frac) as usize; - let (cal_l, test_l) = out.val_logits.split_at(n_cal); - let (cal_y, test_y) = out.val_labels.split_at(n_cal); - - let up_frac = out.report.up_fraction as f32; - let chance_brier = up_frac * (1.0 - up_frac); - let chance_logl = -(up_frac.ln() * up_frac + (1.0 - up_frac).ln() * (1.0 - up_frac)); - - let uncal_brier = brier_score(test_l, test_y); - let uncal_logl = log_loss(test_l, test_y); - - let platt = PlattScaler::fit(cal_l, cal_y, 2000, 1e-2).expect("platt fit"); - let platt_logits = probs_to_logits(&platt.transform(test_l)); - let platt_brier = brier_score(&platt_logits, test_y); - let platt_logl = log_loss(&platt_logits, test_y); - - let iso = IsotonicCalibrator::fit(cal_l, cal_y).expect("iso fit"); - let iso_logits = probs_to_logits(&iso.transform(test_l)); - let iso_brier = brier_score(&iso_logits, test_y); - let iso_logl = log_loss(&iso_logits, test_y); - - println!("\n================================================="); - println!("PHASE 1d.0 — CALIBRATION SMOKE"); - println!("================================================="); - println!("Headline (uncalibrated trainer): acc={:.4} AUC={:.4} n_val={}", - out.report.accuracy, out.report.auc, n_val); - println!("Held-out test n = {} (calibration set = {})", test_l.len(), n_cal); - println!("Chance baselines: Brier = {:.5} log-loss = {:.5}", chance_brier, chance_logl); - println!(); - println!("Uncalibrated: Brier = {:.5} log-loss = {:.5}", uncal_brier, uncal_logl); - println!("Platt scaling: Brier = {:.5} log-loss = {:.5} (a={:.4}, b={:.4})", - platt_brier, platt_logl, platt.a, platt.b); - println!("Isotonic regression: Brier = {:.5} log-loss = {:.5}", iso_brier, iso_logl); - println!(); - let best_brier = platt_brier.min(iso_brier); - if best_brier <= chance_brier { - println!("GATE PASS: best Brier ({:.5}) ≤ chance baseline ({:.5}); proceed to 1d.1.", - best_brier, chance_brier); - } else { - println!("GATE FAIL: best Brier ({:.5}) > chance baseline ({:.5}); falsify 1d.0.", - best_brier, chance_brier); - } - Ok(()) -} diff --git a/crates/ml-alpha/examples/alpha_train_stacker.rs b/crates/ml-alpha/examples/alpha_train_stacker.rs deleted file mode 100644 index 6002b94d3..000000000 --- a/crates/ml-alpha/examples/alpha_train_stacker.rs +++ /dev/null @@ -1,729 +0,0 @@ -//! Phase 1d.2 — Multi-minute label smoke (K=6000 snapshots ≈ 1-5 min forward). -//! -//! THE decisive gate for the two-head architecture. The K-sweep (Phase 1c, -//! commit `db874b184`) showed alpha at stateless single-snapshot resolution -//! decays from K=50 peak to gone by K=500. Mamba2's job is to amplify -//! short-horizon evidence into a long-horizon prediction via SSM state -//! accumulation across `seq_len` snapshots. -//! -//! This test: feed a trailing window of 32 snapshots into Mamba2; predict -//! the binary direction `sign(mid[t+6000] - mid[t])`. Mamba's state must -//! integrate the per-step microstructure signal into a multi-minute call. -//! -//! Gate (per implementation plan): -//! - AUC > 0.55 at K=6000 → multi-minute alpha confirmed; design works. -//! - AUC < 0.52 → DECISIVE FAIL. The two-head architecture cannot turn -//! short-window snapshot context into multi-minute prediction with -//! the current model; need different inputs, model class, or both. -//! -//! Why this is decisive: if we *can* predict long horizons from short -//! tick context plus sequence state, the existing 81-dim snapshot stack -//! is the right foundation and the rest of FoxhuntQ-Δ proceeds. If we -//! can't, the whole architecture is wrong and needs a redesign before -//! any production work. - -use anyhow::{Context, Result}; -use clap::Parser; -use cudarc::driver::CudaContext; -use rand::{Rng, SeedableRng}; -use rand_chacha::ChaCha8Rng; -use std::sync::Arc; -use tracing_subscriber::EnvFilter; - -use ml_alpha::calibration::{Calibrator, IsotonicCalibrator, PlattScaler}; -use ml_alpha::eval::{accuracy_from_logits, auc_from_logits}; -use ml_alpha::fxcache_reader::{FxCacheReader, COL_RAW_CLOSE, FEAT_DIM}; -use ml_alpha::mamba2_block::{ - Mamba2AdamW, Mamba2AdamWConfig, Mamba2Block, Mamba2BlockConfig, -}; -use ml_alpha::metrics_detail::{brier_score, log_loss, stratified_accuracy}; -use ml_alpha::mlp::{MlpConfig, MlpModel}; -use ml_alpha::backtest::GpuBacktest; -use ml_alpha::multi_horizon_labels::generate_labels; - -use ml_core::cuda_autograd::gpu_tensor::GpuTensor; - -#[derive(Parser, Debug)] -#[command(name = "alpha_train_stacker", about = "FoxhuntQ-Δ Phase 1d.2 multi-minute Mamba smoke")] -struct Cli { - #[arg(long)] fxcache_path: String, - /// Label horizon in snapshots. K=6000 ≈ 1-5 min forward at typical rates. - #[arg(long, default_value_t = 6000)] horizon: usize, - /// Sequence length per training example (≤ kernel cap of 32). - #[arg(long, default_value_t = 32)] seq_len: usize, - #[arg(long, default_value_t = 64)] hidden_dim: usize, - #[arg(long, default_value_t = 16)] state_dim: usize, - #[arg(long, default_value_t = 5)] epochs: usize, - #[arg(long, default_value_t = 128)] batch_size: usize, - #[arg(long, default_value_t = 3e-3)] lr: f32, - /// 80/20 train/val split (purged via the embargo). - #[arg(long, default_value_t = 0.8)] train_frac: f32, - /// Embargo bars between train and val ranges (default = horizon for safety). - #[arg(long)] embargo: Option, - /// Subsample stride for train sequences (1 = use every starting position). - #[arg(long, default_value_t = 4)] train_stride: usize, - #[arg(long, default_value_t = 42)] seed: u64, - /// Fraction of val to use as calibration set (rest is held-out test). - /// 0 disables Platt/isotonic post-hoc calibration. - #[arg(long, default_value_t = 0.5)] cal_frac: f32, - /// Stacker hidden dim (small — input is 7-dim, plenty of capacity). - #[arg(long, default_value_t = 32)] stacker_hidden: usize, - /// Stacker training epochs. - #[arg(long, default_value_t = 20)] stacker_epochs: usize, - /// Stacker learning rate. - #[arg(long, default_value_t = 1e-2)] stacker_lr: f32, - /// Stacker batch size. - #[arg(long, default_value_t = 1024)] stacker_batch: usize, - /// Round-trip transaction cost in price units (ES.FUT: 0.25 = 1 tick = $12.50/contract). - #[arg(long, default_value_t = 0.25)] cost_per_trade: f32, - /// Phase E.1 Task 12b: dump per-bar alpha_logits (stacker output) to a - /// binary file consumable by `alpha_dqn_h600_smoke.rs`. The file is a - /// little-endian u32 length prefix followed by n_bars f32 values; bars - /// `< seq_len - 1` are written as 0 (no history). Empty / unset → no - /// cache dumped. - #[arg(long)] alpha_cache_out: Option, - - /// Optional cap on bars consumed from the fxcache. When set, the trainer - /// processes only the first `max_rows` rows; useful for fitting Mamba2 - /// training on a 4 GB consumer GPU without rebuilding a smaller fxcache. - /// Downstream `alpha_baseline --max-snapshots` must be ≤ this value. - #[arg(long)] max_rows: Option, -} - -fn main() -> Result<()> { - tracing_subscriber::fmt() - .with_env_filter(EnvFilter::try_from_default_env().unwrap_or_else(|_| EnvFilter::new("warn"))) - .init(); - let cli = Cli::parse(); - let ctx = CudaContext::new(0).context("init CUDA")?; - let stream = ctx.default_stream(); - - // ── Load fxcache + extract feature matrix and prices ────────────── - let reader = FxCacheReader::open(&cli.fxcache_path)?; - let alpha_dim = reader.alpha_feature_dim() - .ok_or_else(|| anyhow::anyhow!("fxcache lacks alpha column"))?; - let total_bars = reader.bar_count(); - let n_bars = cli.max_rows.map(|m| m.min(total_bars)).unwrap_or(total_bars); - if n_bars < total_bars { - println!( - "fxcache: {} rows total, {} rows used (--max-rows={}), alpha_dim={}, horizon K={}", - total_bars, n_bars, cli.max_rows.unwrap_or(0), alpha_dim, cli.horizon - ); - } else { - println!("fxcache: {} rows, alpha_dim={}, horizon K={}", n_bars, alpha_dim, cli.horizon); - } - - // Flat feature matrix (CPU-resident) and prices. - let mut feature_matrix: Vec = Vec::with_capacity(n_bars * alpha_dim); - let mut prices: Vec = Vec::with_capacity(n_bars); - for i in 0..n_bars { - let row = reader.alpha_features(i).context("alpha row missing")?; - feature_matrix.extend_from_slice(row); - let rec = reader.record(i); - prices.push(rec.targets[COL_RAW_CLOSE - FEAT_DIM]); - } - - // ── Generate K=6000 labels ──────────────────────────────────────── - let labels = generate_labels(&prices, cli.horizon); - println!( - "labels: kept={} dropped_edge={} dropped_invalid={} up_frac={:.4}", - labels.labels.len(), labels.n_dropped_edge, labels.n_dropped_invalid, - labels.labels.iter().sum::() / labels.labels.len().max(1) as f32 - ); - if labels.labels.is_empty() { - anyhow::bail!("no valid labels at K={} — fxcache too small or all-tied", cli.horizon); - } - - // ── Purged train/val split over labels.valid_indices ────────────── - // We split the VALID label list (not the raw bar list) at train_frac, - // then apply an embargo so train sequences ending near the boundary - // don't share forward-window prices with val. - let embargo = cli.embargo.unwrap_or(cli.horizon); - let n_valid = labels.valid_indices.len(); - let n_train_target = (n_valid as f32 * cli.train_frac) as usize; - let train_split_bar = labels.valid_indices[n_train_target - 1]; - let val_start_bar = train_split_bar + embargo; - let train_label_range: Vec = (0..n_valid) - .filter(|&p| labels.valid_indices[p] < train_split_bar.saturating_sub(embargo)) - .collect(); - let val_label_range: Vec = (0..n_valid) - .filter(|&p| labels.valid_indices[p] >= val_start_bar) - .collect(); - println!( - "split: n_train_labels={} n_val_labels={} (embargo={} bars, val_start_bar={})", - train_label_range.len(), val_label_range.len(), embargo, val_start_bar - ); - - // Sequence start: bar `valid_indices[p] - seq_len + 1` to bar `valid_indices[p]`. - // Filter out any p where the sequence would underflow. - let make_seq_label_pos = |label_positions: &[usize], stride: usize| -> Vec { - label_positions.iter() - .filter(|&&p| labels.valid_indices[p] + 1 >= cli.seq_len) - .copied() - .step_by(stride) - .collect() - }; - let train_pos = make_seq_label_pos(&train_label_range, cli.train_stride); - let val_pos = make_seq_label_pos(&val_label_range, 1); - println!("train sequences (stride={}): {} val sequences: {}", - cli.train_stride, train_pos.len(), val_pos.len()); - - // ── Build Mamba2 + AdamW ────────────────────────────────────────── - let block_cfg = Mamba2BlockConfig { - in_dim: alpha_dim, - hidden_dim: cli.hidden_dim, - state_dim: cli.state_dim, - seq_len: cli.seq_len, - }; - let mut block = Mamba2Block::new(block_cfg.clone(), Arc::clone(&stream))?; - let mut opt = Mamba2AdamW::new(&block, Mamba2AdamWConfig { - lr: cli.lr, ..Default::default() - })?; - println!("Mamba2Block params: {}", block.param_count()); - - let gather = |sel: &[usize]| -> Result<(GpuTensor, Vec)> { - let b = sel.len(); - let mut host = Vec::with_capacity(b * cli.seq_len * alpha_dim); - let mut ys = Vec::with_capacity(b); - for &p in sel { - let end_bar = labels.valid_indices[p]; - let start_bar = end_bar + 1 - cli.seq_len; - for t in 0..cli.seq_len { - let off = (start_bar + t) * alpha_dim; - host.extend_from_slice(&feature_matrix[off..off + alpha_dim]); - } - ys.push(labels.labels[p]); - } - let dev = stream.clone_htod(&host)?; - let tensor = GpuTensor::new(dev, vec![b, cli.seq_len, alpha_dim]) - .map_err(|e| anyhow::anyhow!("batch tensor: {e}"))?; - Ok((tensor, ys)) - }; - - let bce_loss = |logits: &[f32], ys: &[f32]| -> f32 { - let mut s = 0.0_f32; - let eps = 1e-7_f32; - for (&z, &y) in logits.iter().zip(ys.iter()) { - let z = z.clamp(-50.0, 50.0); - let p = (1.0 / (1.0 + (-z).exp())).clamp(eps, 1.0 - eps); - s += -(y * p.ln() + (1.0 - y) * (1.0 - p).ln()); - } - s / logits.len() as f32 - }; - - // ── Training loop ───────────────────────────────────────────────── - let mut rng = ChaCha8Rng::seed_from_u64(cli.seed); - let train_len = train_pos.len(); - if train_len < cli.batch_size { - anyhow::bail!("n_train ({train_len}) < batch_size ({})", cli.batch_size); - } - let batches_per_epoch = train_len / cli.batch_size; - - for epoch in 0..cli.epochs { - let mut perm: Vec = (0..train_len).collect(); - for i in (1..train_len).rev() { - let j = rng.gen_range(0..=i); - perm.swap(i, j); - } - let mut loss_sum = 0.0_f32; - for batch_idx in 0..batches_per_epoch { - let raw_sel = &perm[batch_idx * cli.batch_size .. (batch_idx + 1) * cli.batch_size]; - let sel: Vec = raw_sel.iter().map(|&i| train_pos[i]).collect(); - let (input, ys) = gather(&sel)?; - let (logit, cache) = block.forward_train(&input)?; - let logit_host = logit.to_host(&stream)?; - loss_sum += bce_loss(&logit_host, &ys); - - let n_b = ys.len() as f32; - let d_logit_host: Vec = logit_host.iter().zip(ys.iter()) - .map(|(&z, &y)| { - let p = 1.0 / (1.0 + (-z.clamp(-50.0, 50.0)).exp()); - (p - y) / n_b - }) - .collect(); - let d_logit_dev = stream.clone_htod(&d_logit_host)?; - let d_logit = GpuTensor::new(d_logit_dev, vec![ys.len(), 1])?; - let grads = block.backward(&cache, &d_logit)?; - opt.step(&mut block, &grads)?; - } - println!("epoch {epoch:2} mean_train_bce={:.5} n_batches={batches_per_epoch}", - loss_sum / batches_per_epoch as f32); - } - - // ── Validation ──────────────────────────────────────────────────── - let mut val_logits: Vec = Vec::with_capacity(val_pos.len()); - let mut val_ys: Vec = Vec::with_capacity(val_pos.len()); - let mut i = 0; - while i < val_pos.len() { - let this = cli.batch_size.min(val_pos.len() - i); - let sel: Vec = val_pos[i..i + this].to_vec(); - let (input, ys) = gather(&sel)?; - let logit = block.forward(&input)?; - let logit_host = logit.to_host(&stream)?; - val_logits.extend_from_slice(&logit_host); - val_ys.extend_from_slice(&ys); - i += this; - } - let labels_u8: Vec = val_ys.iter().map(|&y| if y > 0.5 { 1 } else { 0 }).collect(); - let acc = accuracy_from_logits(&val_logits, &labels_u8); - let auc = auc_from_logits(&val_logits, &labels_u8); - let uncal_brier = brier_score(&val_logits, &val_ys); - let uncal_logl = log_loss(&val_logits, &val_ys); - let up_frac = val_ys.iter().sum::() / val_ys.len() as f32; - let chance_brier = up_frac * (1.0 - up_frac); - let chance_logl = -(up_frac.ln() * up_frac + (1.0 - up_frac).ln() * (1.0 - up_frac)); - - println!(); - println!("==========================================================="); - println!("PHASE 1d.2 — MAMBA2 MULTI-MINUTE SMOKE (K={})", cli.horizon); - println!("==========================================================="); - println!("val sequences: {}", val_logits.len()); - println!("up_fraction: {:.4} (chance baselines: Brier={:.5} log-loss={:.5})", - up_frac, chance_brier, chance_logl); - println!("--- Uncalibrated ---"); - println!(" accuracy: {:.4}", acc); - println!(" AUC: {:.4}", auc); - println!(" Brier: {:.5}", uncal_brier); - println!(" log-loss: {:.5}", uncal_logl); - - // ── Post-hoc Platt + Isotonic calibration on a held-out cal split ── - if cli.cal_frac > 0.0 && cli.cal_frac < 1.0 { - let n_val = val_logits.len(); - let n_cal = (n_val as f32 * cli.cal_frac) as usize; - let (cal_l, test_l) = val_logits.split_at(n_cal); - let (cal_y, test_y) = val_ys.split_at(n_cal); - let test_labels_u8: Vec = test_y.iter().map(|&y| if y > 0.5 { 1 } else { 0 }).collect(); - - // Headline on test half BEFORE calibration (baseline for comparison). - let test_acc_raw = accuracy_from_logits(test_l, &test_labels_u8); - let test_auc_raw = auc_from_logits(test_l, &test_labels_u8); - let test_brier_raw = brier_score(test_l, test_y); - let test_logl_raw = log_loss(test_l, test_y); - - // Fit + apply Platt. - let platt = PlattScaler::fit(cal_l, cal_y, 2000, 1e-2) - .map_err(|e| anyhow::anyhow!("platt fit: {e}"))?; - let platt_probs = platt.transform(test_l); - // accuracy_from_logits expects logits; convert calibrated probs back to logits. - let to_logits = |probs: &[f32]| -> Vec { - probs.iter().map(|&p| { - let pc = p.clamp(1e-7, 1.0 - 1e-7); - (pc / (1.0 - pc)).ln() - }).collect() - }; - let platt_logits = to_logits(&platt_probs); - let platt_acc = accuracy_from_logits(&platt_logits, &test_labels_u8); - let platt_auc = auc_from_logits(&platt_logits, &test_labels_u8); - let platt_brier = brier_score(&platt_logits, test_y); - let platt_logl = log_loss(&platt_logits, test_y); - - // Fit + apply Isotonic. - let iso = IsotonicCalibrator::fit(cal_l, cal_y) - .map_err(|e| anyhow::anyhow!("iso fit: {e}"))?; - let iso_probs = iso.transform(test_l); - let iso_logits = to_logits(&iso_probs); - let iso_acc = accuracy_from_logits(&iso_logits, &test_labels_u8); - let iso_auc = auc_from_logits(&iso_logits, &test_labels_u8); - let iso_brier = brier_score(&iso_logits, test_y); - let iso_logl = log_loss(&iso_logits, test_y); - - println!("--- Held-out test (cal_frac={:.2}, n_cal={}, n_test={}) ---", - cli.cal_frac, n_cal, test_l.len()); - println!(" Uncalibrated: acc={:.4} AUC={:.4} Brier={:.5} log-loss={:.5}", - test_acc_raw, test_auc_raw, test_brier_raw, test_logl_raw); - println!(" Platt: acc={:.4} AUC={:.4} Brier={:.5} log-loss={:.5} (a={:.4} b={:.4})", - platt_acc, platt_auc, platt_brier, platt_logl, platt.a, platt.b); - println!(" Isotonic: acc={:.4} AUC={:.4} Brier={:.5} log-loss={:.5}", - iso_acc, iso_auc, iso_brier, iso_logl); - - // ── Stacked regime head ────────────────────────────────────── - // GPU-native MLP that takes [mamba_logit, 6 Block-S features] as - // input and predicts the same direction label. Trains on the cal - // half, evaluates on the test half. The cal half is used for both - // Platt and the stacker so all comparisons are on the same - // held-out test bars. - { - const N_BLOCK_S: usize = 6; - const STACKER_IN_DIM: usize = 1 + N_BLOCK_S; // mamba_logit + Block-S - // Block-S column offsets within the snapshot row. - const BS_OFFSETS: [usize; N_BLOCK_S] = [75, 76, 77, 78, 79, 80]; - - // Build a flat [N, STACKER_IN_DIM] feature matrix where N = n_val. - // Each row: [mamba_logit, time_since_trade, time_since_snap, - // book_event_rate, spread_bps, L1_imbalance, micro_mid_drift]. - // Block-S values come from the END BAR of each val sequence. - let mut stacker_inputs: Vec = Vec::with_capacity(val_pos.len() * STACKER_IN_DIM); - for (i, &p) in val_pos.iter().enumerate() { - let end_bar = labels.valid_indices[p]; - stacker_inputs.push(val_logits[i]); - for &off in &BS_OFFSETS { - stacker_inputs.push(feature_matrix[end_bar * alpha_dim + off]); - } - } - // Z-score normalise the Block-S columns on the cal half so the - // MLP's Xavier init is on a sensible scale. The mamba_logit is - // already close to standard-normal scale, so we leave it raw. - let mut col_mean = [0.0_f64; STACKER_IN_DIM]; - let mut col_var = [0.0_f64; STACKER_IN_DIM]; - let n_cal_rows = n_cal; - for r in 0..n_cal_rows { - for c in 1..STACKER_IN_DIM { - col_mean[c] += stacker_inputs[r * STACKER_IN_DIM + c] as f64; - } - } - for c in 1..STACKER_IN_DIM { - col_mean[c] /= n_cal_rows as f64; - } - for r in 0..n_cal_rows { - for c in 1..STACKER_IN_DIM { - let d = stacker_inputs[r * STACKER_IN_DIM + c] as f64 - col_mean[c]; - col_var[c] += d * d; - } - } - let mut col_std = [1.0_f32; STACKER_IN_DIM]; - for c in 1..STACKER_IN_DIM { - let s = (col_var[c] / n_cal_rows as f64).sqrt().max(1e-6); - col_std[c] = s as f32; - } - // Apply normalisation across the FULL val matrix using cal stats. - for r in 0..val_pos.len() { - for c in 1..STACKER_IN_DIM { - let v = stacker_inputs[r * STACKER_IN_DIM + c]; - stacker_inputs[r * STACKER_IN_DIM + c] = - (v - col_mean[c] as f32) / col_std[c]; - } - } - - // Train/test split on the SAME boundary as Platt (first n_cal rows = train, - // rest = test). Temporal split (val_pos is in ascending end-bar order). - let stacker_train_n = n_cal; - let stacker_test_n = val_pos.len() - stacker_train_n; - - // Build and train the stacker MLP. - let stacker_cfg = MlpConfig { - in_dim: STACKER_IN_DIM, - hidden_dim: cli.stacker_hidden, - out_dim: 1, - }; - let mut stacker = MlpModel::new(stacker_cfg, Arc::clone(&stream))?; - stacker.set_learning_rate(cli.stacker_lr); - println!("--- Stacker MLP ---"); - println!(" in_dim={} hidden={} params={}", - STACKER_IN_DIM, cli.stacker_hidden, stacker.param_count()); - - // Permutation for shuffling each epoch. - let mut srng = ChaCha8Rng::seed_from_u64(cli.seed.wrapping_add(7)); - let batch = cli.stacker_batch.min(stacker_train_n); - let batches_per_epoch = stacker_train_n / batch; - for sep in 0..cli.stacker_epochs { - let mut perm: Vec = (0..stacker_train_n).collect(); - for i in (1..stacker_train_n).rev() { - let j = srng.gen_range(0..=i); - perm.swap(i, j); - } - let mut loss_sum = 0.0_f32; - for bi in 0..batches_per_epoch { - let sel = &perm[bi * batch..(bi + 1) * batch]; - let mut feats = Vec::with_capacity(batch * STACKER_IN_DIM); - let mut ys = Vec::with_capacity(batch); - for &r in sel { - let off = r * STACKER_IN_DIM; - feats.extend_from_slice(&stacker_inputs[off..off + STACKER_IN_DIM]); - ys.push(val_ys[r]); - } - let loss = stacker.train_step(&feats, &ys, batch)?; - loss_sum += loss; - } - if sep == 0 || (sep + 1) % 5 == 0 { - println!(" epoch {sep:2} mean_bce={:.5}", loss_sum / batches_per_epoch.max(1) as f32); - } - } - - // Evaluate stacker on test half. - let mut stacker_logits: Vec = Vec::with_capacity(stacker_test_n); - let eval_batch = cli.stacker_batch; - let mut i = stacker_train_n; - while i < val_pos.len() { - let this = eval_batch.min(val_pos.len() - i); - let mut feats = Vec::with_capacity(this * STACKER_IN_DIM); - for r in i..i + this { - let off = r * STACKER_IN_DIM; - feats.extend_from_slice(&stacker_inputs[off..off + STACKER_IN_DIM]); - } - let out = stacker.forward_infer(&feats, this)?; - stacker_logits.extend_from_slice(&out); - i += this; - } - let stacker_test_ys: &[f32] = &val_ys[stacker_train_n..]; - let stacker_labels_u8: Vec = - stacker_test_ys.iter().map(|&y| if y > 0.5 { 1 } else { 0 }).collect(); - let s_acc = accuracy_from_logits(&stacker_logits, &stacker_labels_u8); - let s_auc = auc_from_logits(&stacker_logits, &stacker_labels_u8); - let s_brier = brier_score(&stacker_logits, stacker_test_ys); - let s_logl = log_loss(&stacker_logits, stacker_test_ys); - println!(" Stacker: acc={:.4} AUC={:.4} Brier={:.5} log-loss={:.5} (n_test={})", - s_acc, s_auc, s_brier, s_logl, stacker_test_n); - - // ── Phase E.1 Task 12b: alpha cache export ───────────── - // Run stacker inference on ALL bars (not just val/test) and - // dump the resulting alpha_logits to a binary file. The - // alpha_dqn_h600_smoke.rs binary loads this cache and - // attaches it to SnapshotRow.alpha_logit, giving the - // execution policy a real directional signal. - // - // Bars with end_bar < seq_len - 1 lack history and are - // written as 0.0 (Pearl A sentinel — the env's downstream - // consumers treat 0.0 as "no signal"). - if let Some(out_path) = cli.alpha_cache_out.as_ref() { - use std::io::Write as _; - println!("\n=== Phase E.1 alpha-cache export ==="); - println!(" Computing stacker inference on all {} bars", n_bars); - let mut alpha_cache: Vec = vec![0.0; n_bars]; - let inf_batch = cli.batch_size; - let start = cli.seq_len - 1; - let mut end_bar = start; - while end_bar < n_bars { - let this = inf_batch.min(n_bars - end_bar); - // Build (this, seq_len, alpha_dim) batch - let mut host: Vec = Vec::with_capacity(this * cli.seq_len * alpha_dim); - for b in 0..this { - let eb = end_bar + b; - let start_bar = eb + 1 - cli.seq_len; - for t in 0..cli.seq_len { - let off = (start_bar + t) * alpha_dim; - host.extend_from_slice(&feature_matrix[off..off + alpha_dim]); - } - } - let dev = stream.clone_htod(&host)?; - let tensor = GpuTensor::new(dev, vec![this, cli.seq_len, alpha_dim]) - .map_err(|e| anyhow::anyhow!("alpha-cache batch: {e}"))?; - let logit = block.forward(&tensor)?; - let mamba_logits = logit.to_host(&stream)?; - // Build stacker input: [mamba_logit, normalized Block-S features] - let mut s_in: Vec = Vec::with_capacity(this * STACKER_IN_DIM); - for b in 0..this { - let eb = end_bar + b; - s_in.push(mamba_logits[b]); - for (k, &off) in BS_OFFSETS.iter().enumerate() { - let raw = feature_matrix[eb * alpha_dim + off]; - let norm = (raw - col_mean[k + 1] as f32) / col_std[k + 1]; - s_in.push(norm); - } - } - let s_out = stacker.forward_infer(&s_in, this)?; - for b in 0..this { - alpha_cache[end_bar + b] = s_out[b]; - } - end_bar += this; - if end_bar % 100_000 < inf_batch { - println!(" alpha-cache progress: {}/{}", end_bar, n_bars); - } - } - let mut f = std::fs::File::create(out_path) - .with_context(|| format!("create alpha cache {}", out_path))?; - let n: u32 = n_bars as u32; - f.write_all(&n.to_le_bytes())?; - for v in &alpha_cache { - f.write_all(&v.to_le_bytes())?; - } - let mean_logit: f32 = alpha_cache.iter().sum::() / n_bars as f32; - let mut variance: f32 = 0.0; - for v in &alpha_cache { - let d = v - mean_logit; - variance += d * d; - } - variance /= n_bars as f32; - println!(" Wrote {} f32 entries to {}", n_bars, out_path); - println!(" Stats: mean = {:+.4}, std = {:.4}", mean_logit, variance.sqrt()); - } - - // Stratified accuracy of the STACKER on Block-S features — - // tells us whether the stacker absorbed the regime conditioning - // (uniform accuracy across quintiles) or just learned a sharper - // threshold gate (same Q4-elevated pattern as the raw Mamba). - // ── Phase 1d.4 GPU backtest ────────────────────────────── - // Use the stacker's test-half predictions as the trading - // signal. Apply confidence threshold sweep, compute per-trade - // PnL and Sharpe on GPU. - // - // Trade rule: if `|stacker_prob - 0.5| > τ`, take - // direction = sign(prob - 0.5), enter at end_bar, exit at - // end_bar + horizon. Cost in price units per round-trip. - // - // Test sequences span the SECOND HALF of val_pos; for each, - // we need price[end_bar] and price[end_bar + horizon]. - let stacker_probs: Vec = stacker_logits.iter() - .map(|&z| 1.0_f32 / (1.0 + (-z.clamp(-50.0, 50.0)).exp())) - .collect(); - let mut prices_t_host: Vec = Vec::with_capacity(stacker_test_n); - let mut prices_kt_host: Vec = Vec::with_capacity(stacker_test_n); - for r in stacker_train_n..val_pos.len() { - let end_bar = labels.valid_indices[val_pos[r]]; - let kt_bar = end_bar + cli.horizon; - prices_t_host.push(prices[end_bar]); - prices_kt_host.push(prices[kt_bar]); - } - let probs_dev = stream.clone_htod(&stacker_probs)?; - let prices_t_dev = stream.clone_htod(&prices_t_host)?; - let prices_kt_dev = stream.clone_htod(&prices_kt_host)?; - let bt = GpuBacktest::from_block(&block)?; - - // Compute the test-set wall-clock time span for annualisation. - // Use the first and last end-bar timestamps from the val sequences. - // FxCacheReader stores timestamps in nanoseconds. - let first_test_end_bar = - labels.valid_indices[val_pos[stacker_train_n]]; - let last_test_end_bar = - labels.valid_indices[val_pos[val_pos.len() - 1]]; - let first_ts_ns = reader.record_timestamp(first_test_end_bar); - let last_ts_ns = reader.record_timestamp(last_test_end_bar); - let test_time_span_s = ((last_ts_ns - first_ts_ns) as f64) * 1e-9; - println!(); - println!("--- GPU backtest sweep (n_test={}, span={:.2}s ≈ {:.2}h) ---", - stacker_test_n, test_time_span_s, test_time_span_s / 3600.0); - - // Cost sweep: frictionless / quarter-tick / half-tick / 1-tick / 2-tick. - // 0.25 = 1 tick = $12.50/contract round-trip on ES.FUT. - let thresholds: Vec = vec![0.00, 0.05, 0.10, 0.15, 0.20, 0.25]; - let costs: Vec = vec![0.0, 0.0625, 0.125, 0.25, 0.50]; - let bt_stats = bt.run( - &probs_dev, &prices_t_dev, &prices_kt_dev, - &thresholds, &costs, test_time_span_s, - )?; - - println!(" {:>6} {:>5} {:>9} {:>10} {:>10} {:>9} {:>9} {:>10} {:>11} {:>11}", - "cost", "τ", "n_trades", "mean_ret", "std_ret", - "Sharpe", "hit_rate", "trades/yr", "Sharpe_ann", "total_pnl"); - for s in &bt_stats { - println!(" {:>6.4} {:>5.2} {:>9} {:>10.5} {:>10.5} {:>9.4} {:>9.4} {:>10.0} {:>11.4} {:>11.2}", - s.cost, s.threshold, s.n_trades, s.mean_ret, s.std_ret, - s.sharpe_per_trade, s.hit_rate, s.trades_per_year, - s.sharpe_annualised, s.total_pnl); - } - // Identify best annualised Sharpe per cost band (≥ 100 trades for stat power). - println!(); - println!("BEST OPERATING POINT per cost band (annualised-Sharpe-maximising, n_trades ≥ 100):"); - println!(" {:>6} {:>5} {:>9} {:>11} {:>9} {:>11}", - "cost", "τ", "n_trades", "Sharpe_ann", "hit_rate", "total_pnl"); - for &c in &costs { - let best = bt_stats.iter() - .filter(|s| (s.cost - c).abs() < 1e-6 && s.n_trades >= 100) - .max_by(|a, b| a.sharpe_annualised - .partial_cmp(&b.sharpe_annualised) - .unwrap_or(std::cmp::Ordering::Equal)); - if let Some(b) = best { - println!(" {:>6.4} {:>5.2} {:>9} {:>11.4} {:>9.4} {:>11.2}", - b.cost, b.threshold, b.n_trades, - b.sharpe_annualised, b.hit_rate, b.total_pnl); - } - } - // Overall verdict at the most realistic professional cost: 0.125 (half-tick). - let realistic_best = bt_stats.iter() - .filter(|s| (s.cost - 0.125).abs() < 1e-6 && s.n_trades >= 100) - .max_by(|a, b| a.sharpe_annualised - .partial_cmp(&b.sharpe_annualised) - .unwrap_or(std::cmp::Ordering::Equal)); - println!(); - if let Some(best) = realistic_best { - println!("REALISTIC VERDICT (cost=0.125, half-tick — professional execution):"); - println!(" threshold={:.2} n_trades={} Sharpe_ann={:.4} hit_rate={:.4}", - best.threshold, best.n_trades, best.sharpe_annualised, best.hit_rate); - if best.sharpe_annualised > 2.0 { - println!(" BACKTEST GATE PASS: annualised Sharpe > 2.0 at half-tick cost — deployable."); - } else if best.sharpe_annualised > 0.5 { - println!(" BACKTEST MARGINAL: 0.5 < Sharpe_ann ≤ 2.0 — improve execution or thresholds."); - } else { - println!(" BACKTEST FAIL at realistic cost: signal exists but doesn't survive half-tick friction."); - } - } - // Frictionless upper bound as a diagnostic. - let frictionless = bt_stats.iter() - .filter(|s| (s.cost - 0.0).abs() < 1e-6 && s.n_trades >= 100) - .max_by(|a, b| a.sharpe_annualised - .partial_cmp(&b.sharpe_annualised) - .unwrap_or(std::cmp::Ordering::Equal)); - if let Some(f) = frictionless { - println!("FRICTIONLESS UPPER BOUND (cost=0.0): Sharpe_ann={:.4} at τ={:.2} ({} trades) — model's intrinsic edge", - f.sharpe_annualised, f.threshold, f.n_trades); - } - println!(); - - println!("--- Stacker stratified accuracy (test half, by Block-S feature) ---"); - const BLOCK_S_NAMES_S: [&str; 6] = [ - "time_since_trade_s", "time_since_snap_s", "book_event_rate_per_s", - "spread_bps", "L1_imbalance", "micro_mid_drift", - ]; - let stacker_labels_f32: Vec = stacker_labels_u8.iter().map(|&y| y as f32).collect(); - for (col_off, name) in BLOCK_S_NAMES_S.iter().enumerate() { - let global_col = BS_OFFSETS[col_off]; - let feat_test: Vec = (stacker_train_n..val_pos.len()) - .map(|r| { - let end_bar = labels.valid_indices[val_pos[r]]; - feature_matrix[end_bar * alpha_dim + global_col] - }) - .collect(); - let strats = stratified_accuracy(&stacker_logits, &stacker_labels_f32, &feat_test, 5); - let q4 = strats.last().unwrap(); - let q0 = &strats[0]; - println!(" {name:>24}: Q0_acc={:.4} Q4_acc={:.4} (Q0_n={}, Q4_n={})", - q0.accuracy, q4.accuracy, q0.n, q4.n); - } - } - } - - // ── Stratified accuracy on Block-S features (regime diagnostic) ── - // - // Block-S sits at cols 75..81 of the 81-dim snapshot row. For each - // val sequence, gather the END-BAR's feature value, then stratify - // val accuracy across quintiles of that feature. This tells us - // whether Mamba's K=6000 alpha is regime-conditional (concentrated - // in spread-Q4 / micro-drift-Q4 like the Phase 1c stateless MLP) or - // uniform — which decides whether we need an explicit regime head - // before the backtest. - const BLOCK_S_OFF: usize = 75; - const BLOCK_S_NAMES: [&str; 6] = [ - "time_since_trade_s", - "time_since_snap_s", - "book_event_rate_per_s", - "spread_bps", - "L1_imbalance", - "micro_mid_drift", - ]; - let labels_u8_full: Vec = - val_ys.iter().map(|&y| if y > 0.5 { 1 } else { 0 }).collect(); - let val_labels_f32: Vec = labels_u8_full.iter().map(|&y| y as f32).collect(); - println!(); - println!("--- Stratified val accuracy across Block-S features (5 quintiles) ---"); - for (col_off, name) in BLOCK_S_NAMES.iter().enumerate() { - let global_col = BLOCK_S_OFF + col_off; - let feat: Vec = val_pos.iter() - .map(|&p| { - let end_bar = labels.valid_indices[p]; - feature_matrix[end_bar * alpha_dim + global_col] - }) - .collect(); - let strats = stratified_accuracy(&val_logits, &val_labels_f32, &feat, 5); - println!(); - println!(" ▸ feature: {}", name); - println!(" {:>4} {:>14} {:>14} {:>10} {:>10} {:>12}", - "k", "feat_lo", "feat_hi", "n", "accuracy", "observed"); - for s in &strats { - println!( - " {:>4} {:>14.4} {:>14.4} {:>10} {:>10.4} {:>12.4}", - s.stratum, s.feature_lo, s.feature_hi, s.n, s.accuracy, s.observed_pos_rate - ); - } - } - - println!(); - if auc > 0.55 { - println!("GATE PASS: AUC > 0.55 at K={}; multi-minute alpha confirmed.", cli.horizon); - println!(" Two-head FoxhuntQ-Δ architecture validated."); - } else if auc < 0.52 { - println!("GATE FAIL (decisive): AUC < 0.52; design dead in current form."); - println!(" Need different inputs / model class / context window."); - } else { - println!("GATE MARGINAL: 0.52 ≤ AUC ≤ 0.55; tune hyperparameters or extend seq_len."); - } - Ok(()) -} diff --git a/crates/ml-alpha/examples/gbm_baseline.rs b/crates/ml-alpha/examples/gbm_baseline.rs deleted file mode 100644 index f018e65bb..000000000 --- a/crates/ml-alpha/examples/gbm_baseline.rs +++ /dev/null @@ -1,245 +0,0 @@ -//! Phase 1a GBM baseline (pure-Rust gradient-boosted decision trees via `gbdt`). -//! -//! ## Why -//! -//! The MLP baseline (see `alpha_bar_baseline.rs`) hits ~0.495 validation accuracy on -//! properly-aligned, normalized data. Per the Grinsztajn et al. NeurIPS 2022 -//! survey ("Why do tree-based models still outperform deep learning on typical -//! tabular data?"), gradient-boosted trees are the canonical baseline for -//! ~74-dim engineered tabular features — not MLP. This binary serves as the -//! discriminator between two failure modes: -//! -//! - **GBM > 0.55**: signal exists at bar resolution; MLP-class was the wrong -//! model. FoxhuntQ-Δ proceeds with a tree-ensemble-based alpha head. -//! - **GBM ≈ 0.50**: signal does not exist at bar resolution; the falsification -//! is decisive. Escalate to Phase 1C (tick resolution) per spec. -//! -//! ## Implementation notes -//! -//! - Pure-Rust `gbdt 0.1.3`. No system libLightGBM/libxgboost required. -//! - Loss: `LogLikelyhood` (the only binary-training loss `gbdt` supports; -//! note the upstream typo). Labels converted from {0, 1} → {−1, +1} per -//! that loss's contract. -//! - **No feature normalization** — decision trees are scale-invariant, so -//! we feed the raw `feature_matrix` directly (skip the MLP's z-score step). -//! - Same purged walk-forward split and label generation as the MLP baseline, -//! via the shared `prepare_phase1a_data` helper, so the only varying -//! variable across the two smokes is the model class. -//! -//! ## Usage -//! -//! ``` -//! cargo run -p ml-alpha --example gbm_baseline --release -//! cargo run -p ml-alpha --example gbm_baseline --release -- --iterations 100 --max-depth 8 -//! ``` - -use anyhow::{Context, Result}; -use clap::Parser; -use gbdt::config::Config; -use gbdt::decision_tree::{Data, DataVec, PredVec}; -use gbdt::gradient_boost::GBDT; -use tracing_subscriber::EnvFilter; - -use ml_alpha::eval::{accuracy_from_logits, auc_from_logits, EvalReport}; -use ml_alpha::fxcache_reader::FxCacheReader; -use ml_alpha::purged_split::PurgedSplit; -use ml_alpha::training::{auto_detect_feature_layout, prepare_phase1a_data, Phase1aConfig}; - -#[derive(Debug, Parser)] -#[command(name = "gbm_baseline", about = "FoxhuntQ-Δ Phase 1a GBM baseline (gbdt 0.1.3)")] -struct Cli { - /// Override fxcache path. Defaults to local test_data fxcache. - #[arg(long)] - fxcache_path: Option, - - /// Number of boosting iterations (= number of trees). - #[arg(long, default_value_t = 100)] - iterations: usize, - - /// Maximum tree depth. - #[arg(long, default_value_t = 6)] - max_depth: u32, - - /// Learning rate (shrinkage). - #[arg(long, default_value_t = 0.05)] - shrinkage: f32, - - /// Per-iteration row subsample fraction (bagging). - #[arg(long, default_value_t = 0.8)] - data_sample_ratio: f64, - - /// Per-iteration feature subsample fraction. - #[arg(long, default_value_t = 0.8)] - feature_sample_ratio: f64, - - /// Minimum samples per leaf. - #[arg(long, default_value_t = 50)] - min_leaf_size: usize, - - /// Label horizon in bars. - #[arg(long, default_value_t = 60)] - horizon: usize, - - /// Train fraction for the purged walk-forward split. - #[arg(long, default_value_t = 0.8)] - train_frac: f32, - - /// Additional embargo bars beyond `horizon`. - #[arg(long, default_value_t = 0)] - embargo_bars: usize, -} - -fn main() -> Result<()> { - tracing_subscriber::fmt() - .with_env_filter(EnvFilter::try_from_default_env().unwrap_or_else(|_| EnvFilter::new("info"))) - .init(); - - let cli = Cli::parse(); - - let mut config = Phase1aConfig::default(); - if let Some(p) = cli.fxcache_path.as_deref() { - config.fxcache_path = p.to_string(); - } - config.horizon = cli.horizon; - config.train_frac = cli.train_frac; - config.embargo_bars = cli.embargo_bars; - - tracing::info!(?cli, "Phase 1a GBM baseline starting"); - - // Load + split data using the same pipeline as the MLP baseline (single - // source of truth — only the model class varies between the two smokes). - let reader = FxCacheReader::open(&config.fxcache_path) - .with_context(|| format!("loading fxcache from {}", &config.fxcache_path))?; - let split_cfg = PurgedSplit::new( - reader.bar_count(), - config.train_frac, - config.horizon, - config.embargo_bars, - )?; - let split = split_cfg.split(); - tracing::info!( - n_train_range = split.n_train, - n_val_range = split.n_val, - train = ?split.train, - val = ?split.val, - "purged walk-forward split" - ); - - // Match the MLP baseline's auto-detect: if the fxcache carries the alpha - // column, lift `config.mlp.in_dim` from the legacy 74 to ALPHA_FEATURE_DIM - // (134) before the shared helper asserts the dim against the column. - auto_detect_feature_layout(&reader, &mut config); - - let data = prepare_phase1a_data(&reader, &split, &config)?; - let in_dim = data.in_dim; - - // Build GBDT training samples. `LogLikelyhood` expects labels in {−1, +1} - // per `gbdt::config` docs; our `Phase1aData` labels are in {0.0, 1.0}. - tracing::info!( - n = data.train_indices.len(), - in_dim, - "building GBDT training DataVec (minus-one / plus-one label encoding for LogLikelyhood)" - ); - let mut train_data: DataVec = data - .train_indices - .iter() - .zip(data.train_labels.iter()) - .map(|(&bar_idx, &lbl_01)| { - let feat = data.feature_matrix[bar_idx * in_dim..(bar_idx + 1) * in_dim].to_vec(); - let lbl_pm1 = if lbl_01 > 0.5 { 1.0_f32 } else { -1.0_f32 }; - Data::new_training_data(feat, 1.0, lbl_pm1, None) - }) - .collect(); - - let mut gbdt_cfg = Config::new(); - gbdt_cfg.set_feature_size(in_dim); - gbdt_cfg.set_iterations(cli.iterations); - gbdt_cfg.set_max_depth(cli.max_depth); - gbdt_cfg.set_shrinkage(cli.shrinkage); - gbdt_cfg.set_data_sample_ratio(cli.data_sample_ratio); - gbdt_cfg.set_feature_sample_ratio(cli.feature_sample_ratio); - gbdt_cfg.set_min_leaf_size(cli.min_leaf_size); - gbdt_cfg.set_loss("LogLikelyhood"); - gbdt_cfg.set_training_optimization_level(2); - - tracing::info!( - iterations = cli.iterations, - max_depth = cli.max_depth, - shrinkage = cli.shrinkage, - bag = cli.data_sample_ratio, - feat_subsample = cli.feature_sample_ratio, - "training GBDT (single-threaded, ~minutes)" - ); - let t0 = std::time::Instant::now(); - let mut gbdt = GBDT::new(&gbdt_cfg); - gbdt.fit(&mut train_data); - let train_secs = t0.elapsed().as_secs_f64(); - tracing::info!(train_secs, "GBDT training complete"); - - // Predict on val - tracing::info!(n_val = data.val_indices.len(), "building val DataVec"); - let val_data: DataVec = data - .val_indices - .iter() - .map(|&bar_idx| { - let feat = data.feature_matrix[bar_idx * in_dim..(bar_idx + 1) * in_dim].to_vec(); - Data::new_test_data(feat, None) - }) - .collect(); - let predictions: PredVec = gbdt.predict(&val_data); - - // Accuracy: gbdt's LogLikelyhood outputs probability ∈ [0, 1] of label = +1. - // Threshold at 0.5. `accuracy_from_logits` treats logit > 0 as predict 1, - // which is sign-equivalent to (prob - 0.5) > 0 — translate the probs into - // that convention. - let val_labels_u8: Vec = data - .val_labels - .iter() - .map(|&y| if y > 0.5 { 1 } else { 0 }) - .collect(); - let pred_logits: Vec = predictions.iter().map(|&p| p - 0.5).collect(); - let accuracy = accuracy_from_logits(&pred_logits, &val_labels_u8); - let auc = auc_from_logits(&pred_logits, &val_labels_u8); - - // Also report mean predicted-probability and its histogram to spot the - // common "model is just outputting the prior" failure mode. - let mean_pred = predictions.iter().sum::() / predictions.len().max(1) as f32; - let n_above_half = predictions.iter().filter(|&&p| p > 0.5).count(); - tracing::info!( - n_predictions = predictions.len(), - mean_pred, - n_pred_up = n_above_half, - frac_pred_up = n_above_half as f32 / predictions.len() as f32, - "prediction distribution" - ); - - let report = EvalReport { - n_samples: val_data.len(), - accuracy, - auc, - up_fraction: data.up_fraction_val, - tied_fraction: data.n_val_tied_or_invalid as f32 / split.n_val.max(1) as f32, - note: format!( - "Phase 1a GBM baseline (gbdt 0.1.3 / LogLikelyhood): iter={}, depth={}, \ - η={}, bag={}, feat_subsample={}, train_secs={:.1}", - cli.iterations, - cli.max_depth, - cli.shrinkage, - cli.data_sample_ratio, - cli.feature_sample_ratio, - train_secs - ), - }; - - tracing::info!( - accuracy = report.accuracy, - auc = report.auc, - n_samples = report.n_samples, - up_fraction = report.up_fraction, - "Phase 1a GBM evaluation complete" - ); - println!(); - println!("{}", report.verdict_line()); - - Ok(()) -} diff --git a/crates/ml-alpha/src/backtest.rs b/crates/ml-alpha/src/backtest.rs deleted file mode 100644 index a01141f15..000000000 --- a/crates/ml-alpha/src/backtest.rs +++ /dev/null @@ -1,323 +0,0 @@ -//! Phase 1d.4 — GPU-native backtest. -//! -//! Evaluates a trading policy on the snapshot stream: -//! - Per sequence i with stacker probability `probs[i] ∈ [0, 1]`: -//! if `|probs[i] - 0.5| > τ_conf`, take direction = sign(probs[i] - 0.5), -//! hold for `horizon` snapshots, exit at `price[end_bar + horizon]`. -//! - Per-trade PnL = `direction * (price_kt - price_t) - cost`. -//! - Aggregate across N sequences: trade count, mean return, std return, -//! Sharpe (per-trade, unannualised), hit rate. -//! -//! Sweeps multiple thresholds in a single GPU pass via row-major `[T, N]` -//! kernels (`backtest_per_trade_pnl`) followed by three GPU reductions -//! (`sum_reduce_f32` for returns, `sum_squared_reduce_f32` for returns², -//! `sum_reduce_i32` for trade counts). The final Sharpe/hit-rate -//! arithmetic happens on host once the reductions return scalars — that's -//! 3T scalars to read back regardless of N, so host I/O is bounded -//! independent of dataset size. -//! -//! GPU-pure on the hot path: pinned-mapped buffer transfers for the -//! input arrays (probs/prices), kernel launches for the per-trade math -//! and reductions, no per-trade host roundtrip. Per `feedback_cpu_is_read_only`. - -use std::sync::Arc; - -use anyhow::{anyhow, Result}; -use cudarc::driver::{CudaFunction, CudaSlice, CudaStream, LaunchConfig, PushKernelArg}; - -use crate::mamba2_block::Mamba2Block; - -/// One row of the (cost × threshold) sweep table. -#[derive(Debug, Clone)] -pub struct BacktestStats { - /// Round-trip cost in price units for this row. - pub cost: f32, - /// The confidence threshold this row corresponds to. - pub threshold: f32, - /// Number of trades taken across the test set. - pub n_trades: i32, - /// Mean per-trade PnL (in price units; e.g. for ES.FUT, 1 = 4 ticks = $50). - pub mean_ret: f32, - /// Sample standard deviation of per-trade PnL. - pub std_ret: f32, - /// Mean / std = per-trade Sharpe (unannualised). - pub sharpe_per_trade: f32, - /// Annualised Sharpe = sharpe_per_trade × sqrt(trades_per_year). - /// Assumes trades are roughly i.i.d. across time — overestimates when - /// overlapping signals are taken; underestimates when the model - /// concentrates trades in correlated regimes. Use as a comparable - /// ballpark, not a deployment forecast. - pub sharpe_annualised: f32, - /// Realised trade rate scaled to a calendar year — gives the - /// annualisation factor context. - pub trades_per_year: f32, - /// Fraction of trades with positive PnL after cost. - pub hit_rate: f32, - /// Total PnL across all trades. - pub total_pnl: f32, -} - -/// Seconds in a calendar year, used for annualisation: 365.25 × 86400. -const SECONDS_PER_YEAR: f64 = 365.25 * 86_400.0; - -/// GPU-native backtest evaluator. Holds the kernel handles + a scratch -/// allocator on the provided CUDA stream. Reusable across multiple -/// `run` calls (e.g. for parameter sweeps). -pub struct GpuBacktest { - stream: Arc, - kernel_pnl: CudaFunction, - kernel_sum_f32: CudaFunction, - kernel_sum_sq_f32: CudaFunction, - kernel_sum_i32: CudaFunction, -} - -impl GpuBacktest { - /// Construct from an existing Mamba2Block — reuses the cubin module - /// already loaded by the block (no second cubin load). - pub fn from_block(block: &Mamba2Block) -> Result { - // Resolve the four backtest kernel symbols. The cubin module is - // already held alive by the Mamba2Block; load_function returns - // an independent CudaFunction handle that's safe to keep here. - // We need access to the module — Mamba2Block keeps it private - // (_module field is prefixed with underscore). Instead, we - // re-load the cubin here. The cost is one load on construction; - // subsequent kernel launches reuse the handles. - static CUBIN: &[u8] = - include_bytes!(concat!(env!("OUT_DIR"), "/mamba2_alpha_kernel.cubin")); - let module = block - .stream - .context() - .load_cubin(CUBIN.to_vec()) - .map_err(|e| anyhow!("GpuBacktest: cubin load: {e}"))?; - let kernel_pnl = module - .load_function("backtest_per_trade_pnl") - .map_err(|e| anyhow!("backtest_per_trade_pnl resolve: {e}"))?; - let kernel_sum_f32 = module - .load_function("backtest_sum_reduce_f32") - .map_err(|e| anyhow!("backtest_sum_reduce_f32 resolve: {e}"))?; - let kernel_sum_sq_f32 = module - .load_function("backtest_sum_squared_reduce_f32") - .map_err(|e| anyhow!("backtest_sum_squared_reduce_f32 resolve: {e}"))?; - let kernel_sum_i32 = module - .load_function("backtest_sum_reduce_i32") - .map_err(|e| anyhow!("backtest_sum_reduce_i32 resolve: {e}"))?; - - Ok(Self { - stream: Arc::clone(&block.stream), - kernel_pnl, - kernel_sum_f32, - kernel_sum_sq_f32, - kernel_sum_i32, - }) - } - - /// Run the (cost × threshold) sweep. - /// - /// Arguments: - /// - `probs` : `[N]` stacker sigmoid output ∈ [0, 1] (already on GPU) - /// - `prices_t` : `[N]` mid-prices at sequence end-bar (already on GPU) - /// - `prices_kt` : `[N]` mid-prices at end-bar + horizon (already on GPU) - /// - `thresholds` : `[T]` confidence thresholds (uploaded inside) - /// - `costs` : `[C]` round-trip cost values in price units - /// - `test_time_span_s` : wall-clock seconds spanned by the test set - /// (last end-bar ts − first end-bar ts); used for annualised Sharpe - /// - /// Returns `C × T` `BacktestStats` rows, in (cost-major, threshold-minor) order. - #[allow(clippy::too_many_arguments)] - pub fn run( - &self, - probs: &CudaSlice, - prices_t: &CudaSlice, - prices_kt: &CudaSlice, - thresholds: &[f32], - costs: &[f32], - test_time_span_s: f64, - ) -> Result> { - let mut all_stats: Vec = - Vec::with_capacity(costs.len() * thresholds.len()); - for &cost in costs { - let stats = self.run_single_cost( - probs, prices_t, prices_kt, thresholds, cost, test_time_span_s, - )?; - all_stats.extend(stats); - } - Ok(all_stats) - } - - /// Single-cost backtest. The (cost × threshold) sweep wraps this in - /// an outer loop over costs; pulled out so we can keep all kernel - /// launches and allocations local. - #[allow(clippy::too_many_arguments)] - fn run_single_cost( - &self, - probs: &CudaSlice, - prices_t: &CudaSlice, - prices_kt: &CudaSlice, - thresholds: &[f32], - cost_per_trade: f32, - test_time_span_s: f64, - ) -> Result> { - let n = probs.len() as i32; - let t = thresholds.len() as i32; - if prices_t.len() as i32 != n || prices_kt.len() as i32 != n { - return Err(anyhow!( - "GpuBacktest::run: shape mismatch (probs={}, prices_t={}, prices_kt={})", - probs.len(), prices_t.len(), prices_kt.len() - )); - } - if t == 0 { - return Ok(Vec::new()); - } - - // Upload thresholds. - let thresholds_dev = self.stream - .clone_htod(thresholds) - .map_err(|e| anyhow!("htod thresholds: {e}"))?; - - // Allocate per-trade output buffers [T, N]. - let nt = (n as usize) * (t as usize); - let mut pnl_grid = self.stream - .alloc_zeros::(nt) - .map_err(|e| anyhow!("alloc pnl_grid: {e}"))?; - let mut trade_grid = self.stream - .alloc_zeros::(nt) - .map_err(|e| anyhow!("alloc trade_grid: {e}"))?; - - // Launch per-trade PnL kernel: grid = (T, ceil(N/256)), block = 256 - let block_threads: u32 = 256; - let grid_y = (((n as usize) + block_threads as usize - 1) / block_threads as usize) as u32; - let pnl_cfg = LaunchConfig { - grid_dim: (t as u32, grid_y, 1), - block_dim: (block_threads, 1, 1), - shared_mem_bytes: 0, - }; - unsafe { - self.stream - .launch_builder(&self.kernel_pnl) - .arg(probs) - .arg(prices_t) - .arg(prices_kt) - .arg(&thresholds_dev) - .arg(&cost_per_trade) - .arg(&mut pnl_grid) - .arg(&mut trade_grid) - .arg(&n) - .arg(&t) - .launch(pnl_cfg) - .map_err(|e| anyhow!("backtest_per_trade_pnl launch: {e}"))?; - } - - // GPU reductions: one block per threshold, 256 threads. - let red_cfg = LaunchConfig { - grid_dim: (t as u32, 1, 1), - block_dim: (block_threads, 1, 1), - shared_mem_bytes: 0, - }; - let mut sum_ret = self.stream.alloc_zeros::(t as usize) - .map_err(|e| anyhow!("alloc sum_ret: {e}"))?; - let mut sum_sq = self.stream.alloc_zeros::(t as usize) - .map_err(|e| anyhow!("alloc sum_sq: {e}"))?; - let mut sum_trades = self.stream.alloc_zeros::(t as usize) - .map_err(|e| anyhow!("alloc sum_trades: {e}"))?; - - unsafe { - self.stream - .launch_builder(&self.kernel_sum_f32) - .arg(&pnl_grid) - .arg(&mut sum_ret) - .arg(&n).arg(&t) - .launch(red_cfg) - .map_err(|e| anyhow!("sum_reduce_f32 launch: {e}"))?; - } - unsafe { - self.stream - .launch_builder(&self.kernel_sum_sq_f32) - .arg(&pnl_grid) - .arg(&mut sum_sq) - .arg(&n).arg(&t) - .launch(red_cfg) - .map_err(|e| anyhow!("sum_squared_reduce_f32 launch: {e}"))?; - } - unsafe { - self.stream - .launch_builder(&self.kernel_sum_i32) - .arg(&trade_grid) - .arg(&mut sum_trades) - .arg(&n).arg(&t) - .launch(red_cfg) - .map_err(|e| anyhow!("sum_reduce_i32 launch: {e}"))?; - } - - // Read scalar reductions back to host (3T scalars; bounded independent of N). - let mut sum_ret_h = vec![0.0_f32; t as usize]; - let mut sum_sq_h = vec![0.0_f32; t as usize]; - let mut sum_trades_h = vec![0_i32; t as usize]; - self.stream.memcpy_dtoh(&sum_ret, &mut sum_ret_h) - .map_err(|e| anyhow!("dtoh sum_ret: {e}"))?; - self.stream.memcpy_dtoh(&sum_sq, &mut sum_sq_h) - .map_err(|e| anyhow!("dtoh sum_sq: {e}"))?; - self.stream.memcpy_dtoh(&sum_trades, &mut sum_trades_h) - .map_err(|e| anyhow!("dtoh sum_trades: {e}"))?; - - // Also need hit rate — count of positive trades. We get that with one more - // reduction over the SIGN of pnl_grid > 0. For simplicity, compute on host - // from the pnl_grid (a single dtoh of [T, N] floats — could be large but - // for T=10, N=200K that's 8 MB — acceptable as a one-shot end-of-backtest - // host transfer, not on the hot path). - let mut pnl_grid_h = vec![0.0_f32; nt]; - self.stream.memcpy_dtoh(&pnl_grid, &mut pnl_grid_h) - .map_err(|e| anyhow!("dtoh pnl_grid: {e}"))?; - let mut trade_grid_h = vec![0_i32; nt]; - self.stream.memcpy_dtoh(&trade_grid, &mut trade_grid_h) - .map_err(|e| anyhow!("dtoh trade_grid: {e}"))?; - - // Assemble per-threshold BacktestStats. - let n_u = n as usize; - let mut out = Vec::with_capacity(t as usize); - for ti in 0..(t as usize) { - let nt_trades = sum_trades_h[ti]; - let total = sum_ret_h[ti]; - let mut wins = 0_i32; - for ni in 0..n_u { - let v = pnl_grid_h[ti * n_u + ni]; - if trade_grid_h[ti * n_u + ni] == 1 && v > 0.0 { - wins += 1; - } - } - let (mean_ret, std_ret, sharpe, hit_rate) = if nt_trades > 0 { - let nf = nt_trades as f32; - let mean = total / nf; - let var = (sum_sq_h[ti] / nf - mean * mean).max(0.0); - let sd = var.sqrt(); - let sh = if sd > 1e-9 { mean / sd } else { 0.0 }; - let hit = wins as f32 / nf; - (mean, sd, sh, hit) - } else { - (0.0, 0.0, 0.0, 0.0) - }; - // Annualisation: trades_per_year = n_trades × (year_seconds / test_span_seconds). - // Annualised Sharpe = per_trade_Sharpe × sqrt(trades_per_year) - // — standard practice for per-event Sharpe (Sharpe-time-scaling). - let (trades_per_year, sharpe_annualised) = if test_time_span_s > 0.0 && nt_trades > 0 { - let tpy = (nt_trades as f64) * (SECONDS_PER_YEAR / test_time_span_s); - let ann = sharpe as f64 * tpy.sqrt(); - (tpy as f32, ann as f32) - } else { - (0.0, 0.0) - }; - out.push(BacktestStats { - cost: cost_per_trade, - threshold: thresholds[ti], - n_trades: nt_trades, - mean_ret, - std_ret, - sharpe_per_trade: sharpe, - sharpe_annualised, - trades_per_year, - hit_rate, - total_pnl: total, - }); - } - Ok(out) - } -} diff --git a/crates/ml-alpha/src/calibration.rs b/crates/ml-alpha/src/calibration.rs deleted file mode 100644 index b24b57286..000000000 --- a/crates/ml-alpha/src/calibration.rs +++ /dev/null @@ -1,202 +0,0 @@ -//! Phase 1d.0 — post-hoc probability calibration. -//! -//! Platt scaling: logistic regression `P_calib = σ(a·logit + b)` fit on a -//! held-out calibration set. Preserves rank order (AUC unchanged); fixes -//! the sigmoid threshold and prediction-distribution shape. -//! -//! Isotonic regression: monotonic non-parametric calibrator via PAV -//! (pool-adjacent-violators). Better when miscalibration isn't sigmoidal. - -/// A fitted calibrator that maps raw logits → calibrated probabilities. -pub trait Calibrator { - fn transform(&self, logits: &[f32]) -> Vec; -} - -/// Platt scaling — fit a logistic regression on (logit, label) pairs via -/// gradient descent on BCE loss. Two parameters: `a` (slope) and `b` (intercept). -/// -/// `P_calibrated = σ(a · logit + b)` -#[derive(Debug, Clone)] -pub struct PlattScaler { - pub a: f32, - pub b: f32, -} - -impl PlattScaler { - /// Fit `(a, b)` by minimising mean BCE over the calibration set. - /// Returns the fitted scaler. `max_iters` is the gradient-descent budget; - /// `lr` is the learning rate. - pub fn fit( - logits: &[f32], - labels: &[f32], - max_iters: usize, - lr: f32, - ) -> Result { - if logits.len() != labels.len() { - return Err("Platt: logits and labels length mismatch"); - } - if logits.is_empty() { - return Err("Platt: empty calibration set"); - } - let n = logits.len() as f32; - let mut a = 1.0_f32; - let mut b = 0.0_f32; - for _ in 0..max_iters { - let mut grad_a = 0.0_f32; - let mut grad_b = 0.0_f32; - for (&l, &y) in logits.iter().zip(labels.iter()) { - let z = (a * l + b).clamp(-50.0, 50.0); - let p = 1.0 / (1.0 + (-z).exp()); - let err = p - y; - grad_a += err * l; - grad_b += err; - } - a -= lr * grad_a / n; - b -= lr * grad_b / n; - } - Ok(Self { a, b }) - } -} - -impl Calibrator for PlattScaler { - fn transform(&self, logits: &[f32]) -> Vec { - logits - .iter() - .map(|&l| { - let z = (self.a * l + self.b).clamp(-50.0, 50.0); - 1.0 / (1.0 + (-z).exp()) - }) - .collect() - } -} - -/// Isotonic regression calibrator (pool-adjacent-violators algorithm). -/// -/// Fits a monotone-non-decreasing step function from sorted-logit positions -/// to observed-positive-rate. For new inputs, linear interpolation between -/// the nearest two cut points. -#[derive(Debug, Clone)] -pub struct IsotonicCalibrator { - cuts: Vec, - values: Vec, -} - -impl IsotonicCalibrator { - pub fn fit(logits: &[f32], labels: &[f32]) -> Result { - if logits.len() != labels.len() { - return Err("isotonic: logit/label length mismatch"); - } - if logits.is_empty() { - return Err("isotonic: empty calibration set"); - } - let mut paired: Vec<(f32, f32)> = - logits.iter().copied().zip(labels.iter().copied()).collect(); - paired.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal)); - - let n = paired.len(); - let mut cuts: Vec = paired.iter().map(|p| p.0).collect(); - let mut values: Vec = paired.iter().map(|p| p.1).collect(); - let mut weights: Vec = vec![1.0; n]; - - let mut i = 0; - while i + 1 < values.len() { - if values[i] > values[i + 1] { - let new_val = (values[i] * weights[i] + values[i + 1] * weights[i + 1]) - / (weights[i] + weights[i + 1]); - let new_w = weights[i] + weights[i + 1]; - values[i] = new_val; - weights[i] = new_w; - values.remove(i + 1); - weights.remove(i + 1); - cuts.remove(i + 1); - if i > 0 { - i -= 1; - } - } else { - i += 1; - } - } - Ok(Self { cuts, values }) - } -} - -impl Calibrator for IsotonicCalibrator { - fn transform(&self, logits: &[f32]) -> Vec { - logits - .iter() - .map(|&l| { - if self.cuts.is_empty() { - return 0.5; - } - if l <= self.cuts[0] { - return self.values[0]; - } - if l >= *self.cuts.last().unwrap() { - return *self.values.last().unwrap(); - } - let idx = self.cuts.partition_point(|&c| c <= l).saturating_sub(1); - if idx + 1 >= self.cuts.len() { - return self.values[idx]; - } - let lo = self.cuts[idx]; - let hi = self.cuts[idx + 1]; - let lov = self.values[idx]; - let hiv = self.values[idx + 1]; - let frac = if hi > lo { (l - lo) / (hi - lo) } else { 0.0 }; - lov + frac * (hiv - lov) - }) - .collect() - } -} - -#[cfg(test)] -mod tests { - use super::*; - - #[test] - fn test_module_compiles() { - let _ = std::any::type_name::(); - } - - #[test] - fn test_platt_perfect_separable() { - let logits: Vec = (0..200) - .map(|i| if i % 2 == 0 { 2.0 } else { -2.0 }) - .collect(); - let labels: Vec = (0..200) - .map(|i| if i % 2 == 0 { 1.0 } else { 0.0 }) - .collect(); - // 2000 iters @ lr=1e-2: plenty of budget for the 2-param convex problem. - let cal = PlattScaler::fit(&logits, &labels, 2000, 1e-2).expect("fit"); - let probs = cal.transform(&logits); - // Invariant: positive samples have HIGHER probability than negative ones. - // Don't lock a specific threshold (per pearl_tests_must_prove_not_lock). - let pos_min = probs.iter().enumerate() - .filter(|(i, _)| i % 2 == 0) - .map(|(_, p)| *p) - .fold(f32::INFINITY, f32::min); - let neg_max = probs.iter().enumerate() - .filter(|(i, _)| i % 2 == 1) - .map(|(_, p)| *p) - .fold(f32::NEG_INFINITY, f32::max); - assert!(pos_min > neg_max, - "expected clean separation: min(pos_probs)={pos_min} > max(neg_probs)={neg_max}"); - assert!(pos_min > 0.5, "positive class mean prob must exceed 0.5, got pos_min={pos_min}"); - assert!(neg_max < 0.5, "negative class mean prob must drop below 0.5, got neg_max={neg_max}"); - } - - #[test] - fn test_isotonic_monotone_output() { - let logits: Vec = (0..100).map(|i| i as f32 * 0.1).collect(); - let labels: Vec = (0..100) - .map(|i| if i >= 50 { 1.0 } else { 0.0 }) - .collect(); - let cal = IsotonicCalibrator::fit(&logits, &labels).expect("fit"); - let probs = cal.transform(&logits); - for w in probs.windows(2) { - assert!(w[1] >= w[0] - 1e-6, "monotonicity violated: {} → {}", w[0], w[1]); - } - assert!(probs.last().unwrap() >= &0.5); - assert!(probs.first().unwrap() <= &0.5); - } -} diff --git a/crates/ml-alpha/src/eval.rs b/crates/ml-alpha/src/eval.rs deleted file mode 100644 index c855ecb86..000000000 --- a/crates/ml-alpha/src/eval.rs +++ /dev/null @@ -1,210 +0,0 @@ -//! Phase 1a evaluation: validation accuracy + AUC + gate decision. -//! -//! The Phase 1a gate per FoxhuntQ-Δ v4 spec: -//! - **Pass**: validation accuracy > 0.52 (binary direction prediction). -//! - **Fail**: accuracy ∈ [0.48, 0.52] — no signal at bar resolution. -//! Trigger Phase 1C (tick-resolution) before abandoning FoxhuntQ-Δ. - -use serde::{Deserialize, Serialize}; - -/// Result of a Phase 1a evaluation run. -#[derive(Debug, Clone, Serialize, Deserialize)] -pub struct EvalReport { - /// Number of validation samples evaluated (excluding tied-price drops). - pub n_samples: usize, - /// Validation accuracy on binary direction labels (`up` vs `down`). - pub accuracy: f32, - /// Area under ROC curve. - pub auc: f32, - /// Fraction of `up` labels in validation set (class balance). - pub up_fraction: f32, - /// Fraction of bars whose label was dropped due to tied prices. - pub tied_fraction: f32, - /// Human-readable diagnostic note. - pub note: String, -} - -/// Gate-decision summary for the FoxhuntQ-Δ Phase 1a falsification test. -#[derive(Debug, Clone, Copy, PartialEq, Eq)] -pub enum EvalSummary { - /// Accuracy > 0.55 — strong signal. Proceed straight to Phase 2 (skip 1b). - StrongPass, - /// Accuracy ∈ (0.52, 0.55] — signal present but modest. Proceed to Phase 1b - /// (richer encoders) per FoxhuntQ-Δ v4 spec. - Pass, - /// Accuracy ∈ [0.48, 0.52] — no signal at bar resolution. Phase 1C - /// (tick-resolution) is the next test, OR abandon FoxhuntQ-Δ if Phase 1C - /// also fails. - Fail, - /// Accuracy < 0.48 — anti-signal. This is suspicious (worse than random); - /// usually indicates a label-leakage bug in the opposite direction or a - /// flipped sign convention. Investigate before drawing conclusions. - AntiSignal, -} - -impl EvalReport { - /// Determine the gate verdict from accuracy. - pub fn summary(&self) -> EvalSummary { - if self.accuracy.is_nan() { - // Skeleton commit: training not run yet. - return EvalSummary::Fail; - } - if self.accuracy > 0.55 { - EvalSummary::StrongPass - } else if self.accuracy > 0.52 { - EvalSummary::Pass - } else if self.accuracy >= 0.48 { - EvalSummary::Fail - } else { - EvalSummary::AntiSignal - } - } - - /// Pretty-print verdict for logs. - pub fn verdict_line(&self) -> String { - match self.summary() { - EvalSummary::StrongPass => format!( - "GATE STRONG-PASS: accuracy={:.4} > 0.55 (n={}); skip Phase 1b, proceed to Phase 2", - self.accuracy, self.n_samples - ), - EvalSummary::Pass => format!( - "GATE PASS: accuracy={:.4} > 0.52 (n={}); proceed to Phase 1b for richer encoders", - self.accuracy, self.n_samples - ), - EvalSummary::Fail => format!( - "GATE FAIL: accuracy={:.4} ∈ [0.48, 0.52] (n={}); bar-resolution hypothesis likely \ - confirmed. Trigger Phase 1C (tick-resolution) before abandoning FoxhuntQ-Δ", - self.accuracy, self.n_samples - ), - EvalSummary::AntiSignal => format!( - "GATE ANTI-SIGNAL: accuracy={:.4} < 0.48 (n={}); suspect label-leakage or flipped \ - sign convention. Audit before drawing conclusions.", - self.accuracy, self.n_samples - ), - } - } -} - -/// Compute binary classification accuracy from logits + labels. -/// -/// Pure CPU function; used after pulling validation predictions off the GPU. -pub fn accuracy_from_logits(logits: &[f32], labels: &[u8]) -> f32 { - assert_eq!(logits.len(), labels.len()); - if logits.is_empty() { - return 0.5; - } - let mut correct = 0usize; - for (l, y) in logits.iter().zip(labels.iter()) { - let pred = if *l > 0.0 { 1u8 } else { 0u8 }; - if pred == *y { - correct += 1; - } - } - correct as f32 / logits.len() as f32 -} - -/// Compute AUC via the rank-based estimator (Mann-Whitney U). -/// -/// O(n log n) sort + linear pass. Cheap at our scale. -pub fn auc_from_logits(logits: &[f32], labels: &[u8]) -> f32 { - assert_eq!(logits.len(), labels.len()); - if logits.is_empty() { - return 0.5; - } - let n = logits.len(); - let n_pos = labels.iter().filter(|&&y| y == 1).count(); - let n_neg = n - n_pos; - if n_pos == 0 || n_neg == 0 { - return 0.5; // undefined on single-class set; return neutral - } - - // Sort indices by logit ascending; compute average ranks for the positive class. - let mut idx: Vec = (0..n).collect(); - idx.sort_by(|&a, &b| logits[a].partial_cmp(&logits[b]).unwrap_or(std::cmp::Ordering::Equal)); - - let mut sum_ranks_pos: f64 = 0.0; - let mut i = 0; - while i < n { - let mut j = i + 1; - // group of ties - while j < n && (logits[idx[j]] - logits[idx[i]]).abs() < f32::EPSILON { - j += 1; - } - let avg_rank = (i + j + 1) as f64 / 2.0; // average of [i+1, j] (1-indexed) - for &k in &idx[i..j] { - if labels[k] == 1 { - sum_ranks_pos += avg_rank; - } - } - i = j; - } - // AUC = (sum_ranks_pos − n_pos × (n_pos + 1) / 2) / (n_pos × n_neg) - let auc = (sum_ranks_pos - (n_pos * (n_pos + 1)) as f64 / 2.0) / (n_pos * n_neg) as f64; - auc as f32 -} - -#[cfg(test)] -mod tests { - use super::*; - - #[test] - fn accuracy_perfect() { - // logits >0 → predict 1; <0 → predict 0 - let logits = vec![1.0, -1.0, 2.0, -3.0]; - let labels = vec![1, 0, 1, 0]; - assert_eq!(accuracy_from_logits(&logits, &labels), 1.0); - } - - #[test] - fn accuracy_random() { - let logits = vec![1.0, -1.0, 2.0, -3.0]; - let labels = vec![0, 1, 0, 1]; // all wrong - assert_eq!(accuracy_from_logits(&logits, &labels), 0.0); - } - - #[test] - fn auc_perfect_separation() { - // All positives have higher logits than all negatives → AUC = 1.0 - let logits = vec![-2.0, -1.0, 1.0, 2.0]; - let labels = vec![0, 0, 1, 1]; - let auc = auc_from_logits(&logits, &labels); - assert!((auc - 1.0).abs() < 1e-6, "AUC = {auc}"); - } - - #[test] - fn auc_neutral_split() { - // Positives at mid-range ranks (2, 3): AUC = 0.5 exactly. - // Mann-Whitney U: sum_ranks_pos = 2 + 3 = 5; n_pos = n_neg = 2. - // AUC = (5 − 2×3/2) / (2×2) = 2/4 = 0.5 - let logits = vec![1.0, 2.0, 3.0, 4.0]; - let labels = vec![0, 1, 1, 0]; - let auc = auc_from_logits(&logits, &labels); - assert!((auc - 0.5).abs() < 1e-6, "AUC = {auc}"); - } - - #[test] - fn auc_anti_signal() { - // Positives at LOW ranks → AUC < 0.5. - let logits = vec![1.0, 2.0, 3.0, 4.0]; - let labels = vec![1, 1, 0, 0]; - let auc = auc_from_logits(&logits, &labels); - // sum_ranks_pos = 1+2 = 3; AUC = (3 - 3)/4 = 0.0 - assert!((auc - 0.0).abs() < 1e-6, "AUC = {auc}"); - } - - #[test] - fn verdict_branches() { - let r = |acc: f32| EvalReport { - n_samples: 1000, - accuracy: acc, - auc: 0.5, - up_fraction: 0.5, - tied_fraction: 0.0, - note: String::new(), - }; - assert_eq!(r(0.60).summary(), EvalSummary::StrongPass); - assert_eq!(r(0.53).summary(), EvalSummary::Pass); - assert_eq!(r(0.50).summary(), EvalSummary::Fail); - assert_eq!(r(0.40).summary(), EvalSummary::AntiSignal); - } -} diff --git a/crates/ml-alpha/src/lib.rs b/crates/ml-alpha/src/lib.rs index 66f1697cf..751e0f78e 100644 --- a/crates/ml-alpha/src/lib.rs +++ b/crates/ml-alpha/src/lib.rs @@ -1,66 +1,39 @@ -//! # ml-alpha — FoxhuntQ-Δ Phase 1a +//! ml-alpha — CfC perception + multi-horizon alpha heads. //! -//! Minimal alpha-only crate for the cheapest possible falsification of the -//! bar-resolution signal hypothesis. +//! Phase A (this crate): snapshot-level CfC trunk + 5 horizon heads, +//! trained supervised on MBP-10 with multi-horizon BCE. The CfC trunk +//! must beat the Mamba2 baseline (`mamba2_block`) at every horizon +//! before Phase B (PPO) is allowed to start. //! -//! ## Goal -//! -//! Answer ONE question: does supervised binary direction classification at -//! the existing imbalance-bar resolution exceed validation accuracy > 0.52 -//! on a purged walk-forward held-out fold? -//! -//! - **Yes** → FoxhuntQ-Δ proceeds to Phase 1b (richer encoders) and beyond. -//! - **No** → bar-resolution hypothesis confirmed; FoxhuntQ-Δ does not proceed -//! (see `project_bar_resolution_is_actual_architecture`). -//! -//! ## Architecture -//! -//! ```text -//! FxCache file (175k bars × 80 f32 features) -//! │ -//! ▼ (fxcache_reader.rs — minimal mmap-based reader) -//! Per-bar (features[74], target_close, OFI[32]) → labels via sign(price[t+H] − price[t]) -//! │ -//! ▼ (purged_split.rs — Lopez de Prado purged walk-forward) -//! (train_indices, val_indices) with H-bar embargo -//! │ -//! ▼ (mlp.rs — 2-layer GELU MLP, ml-core GpuLinear primitives) -//! p(direction up | features) ∈ [0, 1] -//! │ -//! ▼ (training.rs — BCE loss, GpuAdamW) -//! Trained model -//! │ -//! ▼ (eval.rs — validation accuracy + AUC) -//! Gate: accuracy > 0.52 → proceed -//! -//! ``` -//! -//! ## Discipline -//! -//! - **No `ml`/`ml-dqn`/`ml-supervised` deps** — keep cold-compile under ~2 min -//! - **Purged validation** — per Lopez de Prado, embargo = H bars (label horizon) -//! - **No tuned constants in adaptive paths** — per `feedback_adaptive_not_tuned` -//! - **GPU-resident training** — but local RTX 3050 Ti (4 GB) is enough at this scale +//! All hot-path arithmetic runs on GPU. The crate owns no CPU compute +//! on model tensors. See `docs/superpowers/specs/2026-05-16-ml-alpha-cfc-ppo-design.md`. -#![warn(clippy::all, clippy::pedantic)] -#![allow(clippy::module_name_repetitions, clippy::cast_precision_loss, clippy::cast_possible_truncation, clippy::missing_errors_doc)] +#![warn(clippy::all)] +#![allow( + clippy::module_name_repetitions, + clippy::cast_precision_loss, + clippy::cast_possible_truncation, + clippy::missing_errors_doc +)] +pub mod cfc; +pub mod heads; +pub mod isv; +pub mod pinned; +pub mod trainer; +pub mod data; +pub mod gate; + +// Preserved from prior crate state — used by Phase A. pub mod fxcache_reader; pub mod purged_split; -pub mod mlp; -pub mod training; -pub mod eval; -pub mod metrics_detail; -pub mod calibration; -pub mod mamba2_block; pub mod multi_horizon_labels; -pub mod backtest; -pub use fxcache_reader::{FxCacheReader, FxCacheRecord, FxCacheMetadata}; -pub use purged_split::{PurgedSplit, SplitIndices}; -pub use mlp::{MlpConfig, MlpModel}; -pub use training::{ - auto_detect_feature_layout, prepare_phase1a_data, Phase1aConfig, Phase1aData, - Phase1aRunOutputs, Phase1aTrainer, -}; -pub use eval::{EvalReport, EvalSummary}; +// Gate reference only — Mamba2 baseline against which CfC must compete. +pub mod mamba2_block; + +pub use cfc::CfcTrunk; +pub use heads::{MultiHorizonHeads, Projection}; +pub use isv::IsvBus; +pub use pinned::{MappedPinnedSnapshotSlot, MappedPinnedFillSlot}; +pub use multi_horizon_labels::{generate_labels, LongHorizonLabels}; diff --git a/crates/ml-alpha/src/metrics_detail.rs b/crates/ml-alpha/src/metrics_detail.rs deleted file mode 100644 index d7e434df8..000000000 --- a/crates/ml-alpha/src/metrics_detail.rs +++ /dev/null @@ -1,258 +0,0 @@ -//! Phase 1c detailed metrics — calibration + stratified accuracy analyses. -//! -//! Layered on top of `Phase1aTrainer::run_full`'s raw val logits + labels. -//! All metrics are computed from logit-space inputs (the trainer's native -//! output); sigmoid conversion happens internally so callers don't need to -//! pre-process. -//! -//! ## Why these four -//! -//! - **Brier score**: mean squared error of predicted probability vs. label. -//! Strictly proper scoring rule, bounded ≤ 0.25 for binary; lower = better. -//! - **Log-loss**: cross-entropy. Penalises overconfident wrong predictions -//! more steeply than Brier; reveals miscalibration. -//! - **Reliability curve**: binned predicted-prob vs observed-positive-rate. -//! Diagnoses *which direction* the model is miscalibrated in. -//! - **Stratified accuracy**: val accuracy within quintile bins of any -//! feature column. Diagnoses *which microstructure regime* carries signal. -//! -//! Together they answer: "Is the AUC-vs-accuracy gap a fixable calibration -//! problem, or does alpha concentrate in specific book regimes?" - -/// Convert a logit to a sigmoid probability, with numerical-stable clamp. -#[inline] -fn sigmoid(x: f32) -> f32 { - // Clamp logit to ±50 to keep exp() finite; sigmoid saturates well before. - let x = x.clamp(-50.0, 50.0); - 1.0 / (1.0 + (-x).exp()) -} - -/// Brier score: mean((p - y)^2) over the val set. Bounded in [0, 1] for -/// binary labels {0, 1}; chance baseline = 0.25 for a 50/50 prior. -pub fn brier_score(logits: &[f32], labels: &[f32]) -> f32 { - assert_eq!(logits.len(), labels.len()); - if logits.is_empty() { - return 0.0; - } - let mut s = 0.0_f64; - for (&l, &y) in logits.iter().zip(labels.iter()) { - let p = sigmoid(l) as f64; - let err = p - (y as f64); - s += err * err; - } - (s / logits.len() as f64) as f32 -} - -/// Binary cross-entropy (log-loss). Lower = better. -pub fn log_loss(logits: &[f32], labels: &[f32]) -> f32 { - assert_eq!(logits.len(), labels.len()); - if logits.is_empty() { - return 0.0; - } - let eps = 1e-7_f64; - let mut s = 0.0_f64; - for (&l, &y) in logits.iter().zip(labels.iter()) { - let p = (sigmoid(l) as f64).clamp(eps, 1.0 - eps); - let yf = y as f64; - s += -(yf * p.ln() + (1.0 - yf) * (1.0 - p).ln()); - } - (s / logits.len() as f64) as f32 -} - -/// One row of the reliability table. -#[derive(Debug, Clone)] -pub struct ReliabilityBin { - /// Predicted-probability bin lower edge (inclusive). - pub bin_lo: f32, - /// Predicted-probability bin upper edge (exclusive, except final bin). - pub bin_hi: f32, - /// Number of val samples whose predicted probability fell in this bin. - pub n: usize, - /// Mean predicted probability among samples in this bin. - pub mean_pred: f32, - /// Observed positive-class rate among samples in this bin. - pub observed_pos_rate: f32, -} - -/// Compute the reliability curve as a histogram of predicted-probability -/// bins. Returns one `ReliabilityBin` per bin (empty bins included so the -/// caller can compare across runs with a stable shape). -pub fn reliability_curve(logits: &[f32], labels: &[f32], n_bins: usize) -> Vec { - assert_eq!(logits.len(), labels.len()); - assert!(n_bins >= 2); - let mut sums = vec![0.0_f64; n_bins]; - let mut hits = vec![0.0_f64; n_bins]; - let mut counts = vec![0_usize; n_bins]; - for (&l, &y) in logits.iter().zip(labels.iter()) { - let p = sigmoid(l); - let mut bin = (p * n_bins as f32) as usize; - if bin >= n_bins { - bin = n_bins - 1; - } - sums[bin] += p as f64; - hits[bin] += y as f64; - counts[bin] += 1; - } - (0..n_bins) - .map(|i| { - let n = counts[i]; - let mean_pred = if n > 0 { (sums[i] / n as f64) as f32 } else { 0.0 }; - let observed = if n > 0 { (hits[i] / n as f64) as f32 } else { 0.0 }; - ReliabilityBin { - bin_lo: i as f32 / n_bins as f32, - bin_hi: (i + 1) as f32 / n_bins as f32, - n, - mean_pred, - observed_pos_rate: observed, - } - }) - .collect() -} - -/// One row of a stratification table. -#[derive(Debug, Clone)] -pub struct StratificationBin { - /// Quintile index, 0..n_strata. - pub stratum: usize, - /// Mean of the stratifying feature within this bin (for interpretability). - pub mean_feature: f32, - /// Lower edge of the feature bin (the (stratum/n_strata)-th quantile). - pub feature_lo: f32, - /// Upper edge of the feature bin. - pub feature_hi: f32, - /// Number of val samples in this stratum. - pub n: usize, - /// Accuracy within this stratum (sigmoid > 0.5 → predict 1). - pub accuracy: f32, - /// Observed positive-class rate within this stratum. - pub observed_pos_rate: f32, -} - -/// Stratify val accuracy by quantile bins of a single feature. -/// -/// `feature_values[i]` must correspond to `logits[i]` / `labels[i]`. The -/// function computes equi-frequency quantile cuts (each stratum holds ≈ n/k -/// samples) then reports the per-stratum accuracy. -pub fn stratified_accuracy( - logits: &[f32], - labels: &[f32], - feature_values: &[f32], - n_strata: usize, -) -> Vec { - assert_eq!(logits.len(), labels.len()); - assert_eq!(logits.len(), feature_values.len()); - assert!(n_strata >= 2); - let n = logits.len(); - if n == 0 { - return Vec::new(); - } - - // Compute quantile boundaries via a sorted copy of the feature values. - let mut sorted: Vec = feature_values.iter().copied().collect(); - sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal)); - let mut cuts: Vec = Vec::with_capacity(n_strata + 1); - cuts.push(sorted[0]); - for q in 1..n_strata { - let idx = (q * n) / n_strata; - cuts.push(sorted[idx.min(n - 1)]); - } - cuts.push(sorted[n - 1]); - - // Assign each sample to a stratum via the cut boundaries. - let mut bin_counts = vec![0_usize; n_strata]; - let mut bin_correct = vec![0_usize; n_strata]; - let mut bin_pos = vec![0_usize; n_strata]; - let mut bin_feat_sum = vec![0.0_f64; n_strata]; - - for ((&l, &y), &v) in logits.iter().zip(labels.iter()).zip(feature_values.iter()) { - // Linear search over cuts is fine: n_strata is tiny (5-10). - let mut bin = n_strata - 1; - for q in 0..n_strata { - if v < cuts[q + 1] { - bin = q; - break; - } - } - let pred = if sigmoid(l) > 0.5 { 1.0 } else { 0.0 }; - bin_counts[bin] += 1; - bin_feat_sum[bin] += v as f64; - if pred == y { - bin_correct[bin] += 1; - } - if y > 0.5 { - bin_pos[bin] += 1; - } - } - - (0..n_strata) - .map(|i| { - let n = bin_counts[i]; - let accuracy = if n > 0 { bin_correct[i] as f32 / n as f32 } else { 0.0 }; - let observed = if n > 0 { bin_pos[i] as f32 / n as f32 } else { 0.0 }; - let mean_feat = if n > 0 { (bin_feat_sum[i] / n as f64) as f32 } else { 0.0 }; - StratificationBin { - stratum: i, - mean_feature: mean_feat, - feature_lo: cuts[i], - feature_hi: cuts[i + 1], - n, - accuracy, - observed_pos_rate: observed, - } - }) - .collect() -} - -#[cfg(test)] -mod tests { - use super::*; - - #[test] - fn test_brier_score_perfect_predictor() { - // Logit = ±50 (saturated) lines up with labels {1, 0, 1, 0}; brier ~0. - let logits = vec![50.0, -50.0, 50.0, -50.0]; - let labels = vec![1.0, 0.0, 1.0, 0.0]; - let b = brier_score(&logits, &labels); - assert!(b < 1e-6, "perfect predictor → Brier ≈ 0, got {b}"); - } - - #[test] - fn test_brier_score_random_predictor() { - // Logit = 0 → prob = 0.5 everywhere. Brier = mean((0.5 - y)^2) = 0.25. - let logits = vec![0.0; 100]; - let labels: Vec = (0..100).map(|i| if i % 2 == 0 { 1.0 } else { 0.0 }).collect(); - let b = brier_score(&logits, &labels); - assert!((b - 0.25).abs() < 1e-5, "random predictor → Brier = 0.25, got {b}"); - } - - #[test] - fn test_log_loss_perfect_predictor() { - let logits = vec![50.0, -50.0, 50.0, -50.0]; - let labels = vec![1.0, 0.0, 1.0, 0.0]; - let ll = log_loss(&logits, &labels); - assert!(ll < 1e-5, "perfect predictor → log_loss ≈ 0, got {ll}"); - } - - #[test] - fn test_reliability_curve_sums_to_n() { - let logits: Vec = (-50..50).map(|i| i as f32 * 0.1).collect(); - let labels: Vec = (0..100).map(|i| (i % 2) as f32).collect(); - let bins = reliability_curve(&logits, &labels, 10); - let total: usize = bins.iter().map(|b| b.n).sum(); - assert_eq!(total, 100); - } - - #[test] - fn test_stratified_accuracy_assigns_all_samples() { - let logits = vec![1.0; 100]; - let labels: Vec = (0..100).map(|i| (i % 2) as f32).collect(); - let features: Vec = (0..100).map(|i| i as f32).collect(); - let bins = stratified_accuracy(&logits, &labels, &features, 5); - let total: usize = bins.iter().map(|b| b.n).sum(); - assert_eq!(total, 100); - // Each stratum should hold roughly 20 samples (equi-frequency). - for b in &bins { - assert!(b.n >= 18 && b.n <= 22, "stratum {} has n={} (expected ~20)", b.stratum, b.n); - } - } -} diff --git a/crates/ml-alpha/src/mlp.rs b/crates/ml-alpha/src/mlp.rs deleted file mode 100644 index 1fd5ab8cf..000000000 --- a/crates/ml-alpha/src/mlp.rs +++ /dev/null @@ -1,288 +0,0 @@ -//! 2-layer GELU MLP baseline for Phase 1a falsification. -//! -//! Architecture: -//! ```text -//! Input: [batch, in_dim=74] -//! │ -//! ▼ Linear (74 → hidden=256) -//! ▼ GELU -//! ▼ Linear (256 → 1) -//! ▼ (logit) -//! Output: [batch, 1] -//! ``` -//! -//! Loss: `bce_with_logits` against {0.0, 1.0} direction labels. The logit is -//! left unactivated so the loss kernel can use the numerically-stable -//! `max(z, 0) - z*y + log1p(exp(-|z|))` form; the gradient `(sigmoid(z) - y)/n` -//! is bounded in `[0, 1/n)` regardless of `|z|`, which avoids the unbounded -//! gradient pathology of MSE-on-±1 for confidently-wrong predictions. -//! -//! Optimizer: AdamW (decoupled weight decay) with default config (lr=3e-4, -//! β1=0.9, β2=0.999, wd=1e-5, grad-norm clip at 10.0). The trainer overrides -//! the learning rate from `Phase1aConfig`. -//! -//! ## Why a 2-layer MLP and not deeper / fancier -//! -//! Per TLOB paper (Berti & Kasneci 2025): "MLPLOB surprisingly outperforms -//! complex SOTAs". A 2-layer GELU MLP is the published baseline that complex -//! transformers must beat to justify their cost. If even this MLP exceeds -//! 0.52 validation accuracy on our setup, signal exists at bar resolution and -//! Phase 1b can pursue richer encoders. If pinned at 0.50-0.52, no encoder -//! complexity will save us — bar-resolution hypothesis is confirmed. - -use std::collections::BTreeMap; -use std::sync::Arc; - -use anyhow::{anyhow, Result}; -use cudarc::cublas::CudaBlas; -use cudarc::driver::CudaStream; -use ml_core::cuda_autograd::{ - activations::ActivationKernels, - linear::GpuLinear, - loss::LossKernels, - optimizer::{AdamWConfig, GpuAdamW}, - GpuTensor, GpuVarStore, -}; - -/// Hyperparameters for the Phase 1a MLP. -#[derive(Debug, Clone, Copy)] -pub struct MlpConfig { - /// Input feature dimension (default: 74 = 42 base + 32 OFI). - pub in_dim: usize, - /// Hidden layer width (default: 256 — matches FoxhuntQ-Δ v4 spec). - pub hidden_dim: usize, - /// Output dimension (default: 1 — single logit for binary direction). - pub out_dim: usize, -} - -impl Default for MlpConfig { - fn default() -> Self { - Self { - in_dim: 74, - hidden_dim: 256, - out_dim: 1, - } - } -} - -/// 2-layer GELU MLP on GPU, using `ml-core`'s cuda_autograd primitives. -/// -/// Owns its `GpuVarStore`, AdamW optimizer state, and cuBLAS / activation / -/// loss kernel handles. Single-stream design — all forward, backward, and -/// optimizer kernels enqueue on the stream passed at construction. -pub struct MlpModel { - pub config: MlpConfig, - store: GpuVarStore, - input_linear: GpuLinear, - output_linear: GpuLinear, - activations: ActivationKernels, - loss_kernels: LossKernels, - optimizer: GpuAdamW, - cublas: CudaBlas, - stream: Arc, - // Registered parameter names (mirrors what `store.linear()` registered). - // Used to key the gradient BTreeMap that AdamW looks up. - input_weight_name: String, - input_bias_name: String, - output_weight_name: String, - output_bias_name: String, -} - -impl MlpModel { - /// Construct a fresh MLP with Xavier-initialized weights. - /// - /// Registers two linear layers under the prefixes `"input"` and `"output"` - /// in the var store. AdamW lazy-allocates moment buffers on the first - /// `train_step`, so this constructor is cheap. - pub fn new(config: MlpConfig, stream: Arc) -> Result { - let mut store = GpuVarStore::new(Arc::clone(&stream)); - - let input_linear = store - .linear("input", config.in_dim, config.hidden_dim) - .map_err(|e| anyhow!("input_linear init: {e}"))?; - let output_linear = store - .linear("output", config.hidden_dim, config.out_dim) - .map_err(|e| anyhow!("output_linear init: {e}"))?; - - let activations = ActivationKernels::new(&stream) - .map_err(|e| anyhow!("activations init: {e}"))?; - let loss_kernels = LossKernels::new(&stream) - .map_err(|e| anyhow!("loss_kernels init: {e}"))?; - let optimizer = GpuAdamW::new(AdamWConfig::default(), Arc::clone(&stream)) - .map_err(|e| anyhow!("optimizer init: {e}"))?; - let cublas = CudaBlas::new(Arc::clone(&stream)) - .map_err(|e| anyhow!("cublas init: {e}"))?; - - Ok(Self { - config, - store, - input_linear, - output_linear, - activations, - loss_kernels, - optimizer, - cublas, - stream, - input_weight_name: String::from("input.weight"), - input_bias_name: String::from("input.bias"), - output_weight_name: String::from("output.weight"), - output_bias_name: String::from("output.bias"), - }) - } - - /// Override AdamW learning rate (the constructor used the AdamWConfig default). - pub fn set_learning_rate(&mut self, lr: f32) { - self.optimizer.set_learning_rate(lr); - } - - /// Number of trainable parameters. - pub fn param_count(&self) -> usize { - let c = &self.config; - c.hidden_dim * (c.in_dim + 1) + c.out_dim * (c.hidden_dim + 1) - } - - /// Inference forward pass (no saved activations). - /// - /// `features` is `[batch_size × in_dim]` row-major; returns the raw logits - /// (one per sample) as a `Vec` of length `batch_size`. The caller - /// thresholds at 0 to get the predicted direction. - pub fn forward_infer(&self, features: &[f32], batch_size: usize) -> Result> { - let x = GpuTensor::from_host( - features, - vec![batch_size, self.config.in_dim], - &self.stream, - ) - .map_err(|e| anyhow!("upload features: {e}"))?; - - let (h1_pre, _) = self - .input_linear - .forward(&x, &self.store, &self.cublas, &self.stream) - .map_err(|e| anyhow!("input_linear forward: {e}"))?; - let (h1_post, _) = self - .activations - .gelu_fwd(&h1_pre, &self.stream) - .map_err(|e| anyhow!("gelu_fwd: {e}"))?; - let (z, _) = self - .output_linear - .forward(&h1_post, &self.store, &self.cublas, &self.stream) - .map_err(|e| anyhow!("output_linear forward: {e}"))?; - - z.to_host(&self.stream) - .map_err(|e| anyhow!("logits to_host: {e}")) - } - - /// One training step: forward → BCE-with-logits → backward → AdamW. - /// - /// `features` is `[batch_size × in_dim]` row-major; `labels` is - /// `[batch_size]` in `{0.0, 1.0}` (soft labels in `[0, 1]` are also valid). - /// Returns the mean BCE loss over the batch. - pub fn train_step( - &mut self, - features: &[f32], - labels: &[f32], - batch_size: usize, - ) -> Result { - let x = GpuTensor::from_host( - features, - vec![batch_size, self.config.in_dim], - &self.stream, - ) - .map_err(|e| anyhow!("upload features: {e}"))?; - let y = GpuTensor::from_host( - labels, - vec![batch_size, self.config.out_dim], - &self.stream, - ) - .map_err(|e| anyhow!("upload labels: {e}"))?; - - // ── Forward (save activations for backward) ──────────────── - let (h1_pre, input_acts) = self - .input_linear - .forward(&x, &self.store, &self.cublas, &self.stream) - .map_err(|e| anyhow!("input_linear forward: {e}"))?; - let (h1_post, saved_h1_pre) = self - .activations - .gelu_fwd(&h1_pre, &self.stream) - .map_err(|e| anyhow!("gelu_fwd: {e}"))?; - let (z, output_acts) = self - .output_linear - .forward(&h1_post, &self.store, &self.cublas, &self.stream) - .map_err(|e| anyhow!("output_linear forward: {e}"))?; - - // ── Loss + dL/dz ─────────────────────────────────────────── - let loss_result = self - .loss_kernels - .bce_with_logits(&z, &y, &self.stream) - .map_err(|e| anyhow!("bce_with_logits: {e}"))?; - let loss_host = loss_result - .loss - .to_host(&self.stream) - .map_err(|e| anyhow!("loss to_host: {e}"))?; - let loss_value = loss_host[0]; - - // ── Backward through output → GELU → input ───────────────── - let out_grads = self - .output_linear - .backward( - &loss_result.grad, - &output_acts, - &self.store, - &self.cublas, - &self.stream, - ) - .map_err(|e| anyhow!("output_linear backward: {e}"))?; - let dh1_pre = self - .activations - .gelu_bwd(&out_grads.dx, &saved_h1_pre, &self.stream) - .map_err(|e| anyhow!("gelu_bwd: {e}"))?; - let in_grads = self - .input_linear - .backward( - &dh1_pre, - &input_acts, - &self.store, - &self.cublas, - &self.stream, - ) - .map_err(|e| anyhow!("input_linear backward: {e}"))?; - - // ── AdamW step: key gradients by registered parameter names ─ - let mut grads: BTreeMap = BTreeMap::new(); - grads.insert(self.input_weight_name.clone(), in_grads.dw); - grads.insert(self.input_bias_name.clone(), in_grads.db); - grads.insert(self.output_weight_name.clone(), out_grads.dw); - grads.insert(self.output_bias_name.clone(), out_grads.db); - self.optimizer - .step(&mut self.store, &grads) - .map_err(|e| anyhow!("adamw step: {e}"))?; - - Ok(loss_value) - } -} - -#[cfg(test)] -mod tests { - use super::*; - - /// Verifies the param-count math without touching GPU resources. The - /// formula must match `MlpModel::param_count`; the GPU-backed model is - /// tested via the example binary on the local fxcache. - #[test] - fn param_count_math_matches_default() { - // Default: 74 → 256 → 1 - // input: 256 × 74 + 256 = 19_200 - // output: 1 × 256 + 1 = 257 - // total: 19_457 - let cfg = MlpConfig::default(); - let expected = cfg.hidden_dim * (cfg.in_dim + 1) + cfg.out_dim * (cfg.hidden_dim + 1); - assert_eq!(expected, 19_457); - } - - #[test] - fn config_default_dims() { - let cfg = MlpConfig::default(); - assert_eq!(cfg.in_dim, 74); - assert_eq!(cfg.hidden_dim, 256); - assert_eq!(cfg.out_dim, 1); - } -} diff --git a/crates/ml-alpha/src/training.rs b/crates/ml-alpha/src/training.rs deleted file mode 100644 index c847ef042..000000000 --- a/crates/ml-alpha/src/training.rs +++ /dev/null @@ -1,699 +0,0 @@ -//! Phase 1a training loop. -//! -//! Orchestrates: data load → purged split → label generation → MLP train → -//! validation pass. Single-fold, no walk-forward replication (that's Phase 2+). -//! -//! ## Label source -//! -//! Direction labels come from `raw_close` (target column 2), not `preproc_close` -//! (column 0). The preproc column is log-return normalised across the bar -//! sequence and can lose sign information for very small returns; `raw_close` -//! is the canonical "honest signal" column per `fxcache_reader::COL_RAW_CLOSE`. -//! Tied raw_close pairs (no price movement at horizon H) are dropped, not -//! coerced — they're microstructure noise, not signal. -//! -//! ## Pipeline overview -//! -//! ```text -//! fxcache (mmap) ─→ feature matrix [bars × 74] (one-shot, kept on CPU) -//! ┌─→ prices: raw_close per bar -//! └─→ labels[t] = sign(price[t+H] − price[t]) -//! -//! PurgedSplit.split() ─→ train range, val range -//! -//! Filter to valid (non-tied) labels → train_indices, val_indices -//! -//! For each epoch: -//! ChaCha-shuffle(train_indices) -//! For each minibatch: -//! gather features [batch_size × 74] + labels [batch_size] -//! upload to GPU, train_step → BCE loss, grads, AdamW step -//! -//! Validation: forward_infer over entire val set in chunks → logits -//! Compute accuracy + AUC → EvalReport with verdict. -//! ``` - -use std::sync::Arc; - -use anyhow::{Context, Result}; -use cudarc::driver::CudaStream; -use rand::{Rng, SeedableRng}; -use rand_chacha::ChaCha8Rng; -use serde::{Deserialize, Serialize}; - -use crate::eval::{accuracy_from_logits, auc_from_logits, EvalReport}; -use crate::fxcache_reader::{FxCacheReader, COL_RAW_CLOSE, FEAT_DIM, OFI_DIM}; -use crate::mlp::{MlpConfig, MlpModel}; -use crate::purged_split::{binary_direction_label, PurgedSplit, SplitIndices}; - -/// Phase 1a training configuration. -#[derive(Debug, Clone, Serialize, Deserialize)] -pub struct Phase1aConfig { - pub fxcache_path: String, - pub epochs: usize, - pub batch_size: usize, - pub learning_rate: f32, - pub train_frac: f32, - pub horizon: usize, - pub embargo_bars: usize, - pub mlp: MlpConfigSerde, - pub seed: u64, -} - -/// Serializable mirror of `MlpConfig` (the in-crate type is `Copy` and not -/// `Serialize`-derived; we mirror here for TOML config files). -#[derive(Debug, Clone, Serialize, Deserialize)] -pub struct MlpConfigSerde { - pub in_dim: usize, - pub hidden_dim: usize, - pub out_dim: usize, -} - -impl From for MlpConfig { - fn from(s: MlpConfigSerde) -> Self { - Self { - in_dim: s.in_dim, - hidden_dim: s.hidden_dim, - out_dim: s.out_dim, - } - } -} - -impl Default for Phase1aConfig { - fn default() -> Self { - Self { - fxcache_path: String::from("/home/jgrusewski/Work/foxhunt/test_data/feature-cache/13c0b086a975cc7e2384377a2cd0e97738c9410292fcfecb5807c29bf885cb48.fxcache"), - epochs: 5, - batch_size: 1024, - learning_rate: 1e-3, - train_frac: 0.8, - horizon: 60, - embargo_bars: 0, - mlp: MlpConfigSerde { - in_dim: 74, - hidden_dim: 256, - out_dim: 1, - }, - seed: 42, - } - } -} - -/// Phase 1a orchestrator. Encapsulates the load → split → train → eval cycle. -pub struct Phase1aTrainer { - pub config: Phase1aConfig, - pub reader: FxCacheReader, - pub split: SplitIndices, - pub model: MlpModel, - pub stream: Arc, -} - -impl Phase1aTrainer { - /// Initialize trainer from config. Loads fxcache, computes purged split, - /// constructs MLP. - pub fn from_config(mut config: Phase1aConfig, stream: Arc) -> Result { - let reader = FxCacheReader::open(&config.fxcache_path) - .with_context(|| format!("loading fxcache from {}", config.fxcache_path))?; - let total_bars = reader.bar_count(); - tracing::info!( - total_bars, - version = reader.metadata().version, - has_alpha_features = reader.has_alpha_features(), - "fxcache loaded" - ); - - let split_cfg = PurgedSplit::new( - total_bars, - config.train_frac, - config.horizon, - config.embargo_bars, - )?; - let split = split_cfg.split(); - tracing::info!( - n_train = split.n_train, - n_val = split.n_val, - train = ?split.train, - val = ?split.val, - "purged walk-forward split" - ); - - auto_detect_feature_layout(&reader, &mut config); - - let mut model = MlpModel::new(config.mlp.clone().into(), Arc::clone(&stream))?; - model.set_learning_rate(config.learning_rate); - tracing::info!( - param_count = model.param_count(), - in_dim = model.config.in_dim, - hidden_dim = model.config.hidden_dim, - lr = config.learning_rate, - "MLP initialized" - ); - - Ok(Self { - config, - reader, - split, - model, - stream, - }) - } - - /// Run the smoke and return only the gate-decision report. Calls - /// [`Self::run_full`] under the hood and drops the raw predictions — - /// existing call sites that only need accuracy/AUC keep working. - pub fn run(&mut self) -> Result { - self.run_full().map(|out| out.report) - } - - /// Run the full Phase 1a smoke: train for `epochs` epochs, evaluate on the - /// validation set, return the gate-decision report **plus** the raw val - /// logits, true labels, and the original-bar indices the val samples come - /// from. This is the entry point for the detailed metrics example - /// (calibration + stratified accuracy) — generic enough that any - /// post-hoc analysis can be layered on top of the same predictions. - pub fn run_full(&mut self) -> Result { - let in_dim = self.config.mlp.in_dim; - let data = prepare_phase1a_data(&self.reader, &self.split, &self.config)?; - debug_assert_eq!(data.in_dim, in_dim); - let mut feature_matrix = data.feature_matrix; - let train_indices = data.train_indices; - let train_labels = data.train_labels; - let val_indices = data.val_indices; - let val_labels = data.val_labels; - let val_up_frac = data.up_fraction_val; - let val_tied = data.n_val_tied_or_invalid; - // (data-stage tracing already emitted inside `prepare_phase1a_data`.) - - if train_indices.len() < self.config.batch_size { - anyhow::bail!( - "n_train ({}) < batch_size ({}); reduce --batch-size or use a larger dataset", - train_indices.len(), - self.config.batch_size - ); - } - - // ── Train-only z-score normalization ────────────────────────── - // Mirror `crates/ml/src/walk_forward.rs::NormStats` convention: - // f64 accumulators for numerical stability, MIN_STD = 1e-8, and - // post-divide clip to ±NORMALIZED_FEATURE_BOUND so a single outlier - // can't propagate to ±sqrt(N)·magnitude downstream. Crucially fit on - // TRAIN INDICES ONLY — fitting on the full set leaks validation - // statistics into the training-time baseline (Lopez de Prado, Audit - // Rec 3 in `docs/lookahead-bias-audit-2026-04-28.md`). - let norm_stats = NormStats::fit_train(&feature_matrix, in_dim, &train_indices); - norm_stats.apply_inplace(&mut feature_matrix, in_dim); - let (n_nan_post, n_inf_post, fmin_post, fmax_post, fmean_post, fstd_post) = - feature_stats(&feature_matrix); - tracing::info!( - n_nan_post, n_inf_post, - fmin_post, fmax_post, fmean_post, fstd_post, - "normalization complete (train-only z-score, clipped to ±{})", - NormStats::FEATURE_BOUND - ); - - // ── Train loop ───────────────────────────────────────────────── - let mut rng = ChaCha8Rng::seed_from_u64(self.config.seed); - let batch_size = self.config.batch_size; - let batch_buf_feats: usize = batch_size * in_dim; - let mut batch_features: Vec = vec![0.0; batch_buf_feats]; - let mut batch_labels_buf: Vec = vec![0.0; batch_size]; - - let mut shuffled: Vec = (0..train_indices.len()).collect(); - - for epoch in 0..self.config.epochs { - // Fisher-Yates shuffle of position-indices into train_indices. - for i in (1..shuffled.len()).rev() { - let j = rng.gen_range(0..=i); - shuffled.swap(i, j); - } - - let n_batches = shuffled.len() / batch_size; // drop last partial batch - let mut epoch_loss_sum = 0.0_f64; - for b in 0..n_batches { - let batch_start = b * batch_size; - for k in 0..batch_size { - let pos = shuffled[batch_start + k]; - let bar_idx = train_indices[pos]; - let src = &feature_matrix[bar_idx * in_dim..(bar_idx + 1) * in_dim]; - batch_features[k * in_dim..(k + 1) * in_dim].copy_from_slice(src); - batch_labels_buf[k] = train_labels[pos]; - } - let loss = self - .model - .train_step(&batch_features, &batch_labels_buf, batch_size) - .with_context(|| format!("train_step epoch={epoch} batch={b}"))?; - epoch_loss_sum += loss as f64; - } - let mean_loss = (epoch_loss_sum / n_batches as f64) as f32; - tracing::info!( - epoch, - mean_bce_loss = mean_loss, - n_batches, - "epoch complete" - ); - } - - // ── Validation forward pass (chunked) ───────────────────────── - // 35k × 74 × 4 ≈ 10 MB input fits easily; chunk anyway at 8192 to keep - // peak GPU memory predictable on RTX 3050 Ti (4 GB). - let val_n = val_indices.len(); - let mut val_logits: Vec = Vec::with_capacity(val_n); - let val_chunk = 8192usize.min(val_n); - let mut chunk_features: Vec = vec![0.0; val_chunk * in_dim]; - - let mut i = 0; - while i < val_n { - let this_chunk = val_chunk.min(val_n - i); - for k in 0..this_chunk { - let bar_idx = val_indices[i + k]; - let src = &feature_matrix[bar_idx * in_dim..(bar_idx + 1) * in_dim]; - chunk_features[k * in_dim..(k + 1) * in_dim].copy_from_slice(src); - } - let logits = self - .model - .forward_infer(&chunk_features[..this_chunk * in_dim], this_chunk) - .with_context(|| format!("validation forward chunk i={i}"))?; - val_logits.extend_from_slice(&logits); - i += this_chunk; - } - - // ── Metrics ─────────────────────────────────────────────────── - let labels_u8: Vec = val_labels.iter().map(|&y| if y > 0.5 { 1 } else { 0 }).collect(); - let accuracy = accuracy_from_logits(&val_logits, &labels_u8); - let auc = auc_from_logits(&val_logits, &labels_u8); - let tied_fraction = val_tied as f32 / self.split.n_val.max(1) as f32; - - let report = EvalReport { - n_samples: val_n, - accuracy, - auc, - up_fraction: val_up_frac, - tied_fraction, - note: format!( - "Phase 1a smoke: trained {} epochs, batch={}, lr={}, H={} bars; \ - train n={}, val n={} (raw_close labels)", - self.config.epochs, - self.config.batch_size, - self.config.learning_rate, - self.config.horizon, - train_indices.len(), - val_n, - ), - }; - Ok(Phase1aRunOutputs { - report, - val_logits, - val_labels, - val_indices, - }) - } -} - -/// Rich return value for [`Phase1aTrainer::run_full`]. Carries both the -/// gate-decision report and the raw validation predictions, so callers can -/// run post-hoc analyses (calibration curves, stratified accuracy, -/// reliability tables) without re-training. -pub struct Phase1aRunOutputs { - pub report: EvalReport, - /// Raw MLP logits over the validation set (length = n_val). - pub val_logits: Vec, - /// Ground-truth binary labels for the same val samples (length = n_val). - pub val_labels: Vec, - /// Original-bar indices (into the fxcache row space) for each val sample. - /// Use these to look up un-normalized feature values when stratifying. - pub val_indices: Vec, -} - -/// Phase 1a data bundle returned by [`prepare_phase1a_data`]. -/// -/// Owns the CPU-resident feature matrix (raw, unnormalized — callers that -/// need z-score apply it themselves; tree models that are scale-invariant -/// consume the matrix as-is) and the filtered (train_idx, train_label) + -/// (val_idx, val_label) pairs. Both the MLP trainer and the GBM baseline -/// example call into this function so the data path is shared and audited -/// in exactly one place. -pub struct Phase1aData { - /// Flat row-major matrix of `bar_count × in_dim` raw f32 features. - pub feature_matrix: Vec, - /// Feature dimension (mirrors `FEAT_DIM + OFI_DIM`). - pub in_dim: usize, - /// Total bar count from the fxcache header. - pub bar_count: usize, - /// Bar indices in the training half whose features AND label are valid. - pub train_indices: Vec, - /// Parallel labels in {0.0, 1.0} for `train_indices`. - pub train_labels: Vec, - /// Bar indices in the validation half whose features AND label are valid. - pub val_indices: Vec, - /// Parallel labels in {0.0, 1.0} for `val_indices`. - pub val_labels: Vec, - /// Fraction of `train_labels` equal to 1.0 (the "up" class). - pub up_fraction_train: f32, - /// Fraction of `val_labels` equal to 1.0. - pub up_fraction_val: f32, - /// Bars in the val range dropped due to either feature corruption OR - /// price-tie at horizon. Used to populate `EvalReport::tied_fraction`. - pub n_val_tied_or_invalid: usize, -} - -/// Select feature source for Phase 1a from the fxcache contents. -/// -/// If the cache carries the alpha column, override `config.mlp.in_dim` to the -/// 134-dim modern stack; otherwise validate that the caller is still -/// configured for the 74-dim legacy (FEAT_DIM + OFI_DIM) layout. Both -/// `Phase1aTrainer::from_config` and the GBM baseline example call this -/// before constructing the model / calling `prepare_phase1a_data`, so the -/// auto-detect lives in one place. -pub fn auto_detect_feature_layout(reader: &FxCacheReader, config: &mut Phase1aConfig) { - if let Some(dim) = reader.alpha_feature_dim() { - config.mlp.in_dim = dim; - tracing::info!( - in_dim = dim, - "using alpha feature column from fxcache (variable dim, declared in metadata)" - ); - } else { - assert_eq!( - config.mlp.in_dim, - FEAT_DIM + OFI_DIM, - "legacy fxcache without alpha column requires in_dim = FEAT_DIM + OFI_DIM = {}", - FEAT_DIM + OFI_DIM - ); - tracing::info!( - in_dim = config.mlp.in_dim, - "using legacy 74-dim feature layout (fxcache has no alpha column)" - ); - } -} - -/// Shared data pipeline for Phase 1a baselines (MLP and GBM). -/// -/// Steps: -/// 1. Extract `[features (42) || ofi (32)]` per bar into a flat f32 matrix. -/// 2. Extract `raw_close` per bar for label generation. -/// 3. Audit feature validity (NaN/Inf or `|x| > 1e6` sentinel) and emit a -/// per-column + per-decile corruption report — non-fatal, valid bars are -/// just filtered. -/// 4. Compute binary direction labels via `sign(raw_close[t+H] − raw_close[t])` -/// and drop tied-price bars (microstructure noise per -/// `purged_split::binary_direction_label` semantics). -/// 5. Filter to bars with BOTH valid features AND non-tied labels for train -/// and val ranges. -/// -/// The output is raw (unnormalized) features so consumers can choose their -/// own normalization strategy: the MLP path applies train-only z-score -/// inside the trainer; the GBM path skips it (decision trees are scale- -/// invariant). -pub fn prepare_phase1a_data( - reader: &FxCacheReader, - split: &SplitIndices, - config: &Phase1aConfig, -) -> Result { - let in_dim = config.mlp.in_dim; - let alpha_dim_opt = reader.alpha_feature_dim(); - let use_alpha = alpha_dim_opt.is_some(); - if let Some(dim) = alpha_dim_opt { - assert_eq!( - in_dim, dim, - "in_dim ({in_dim}) must equal fxcache's declared alpha_feature_dim ({dim})" - ); - } else { - assert_eq!( - in_dim, - FEAT_DIM + OFI_DIM, - "in_dim ({in_dim}) must equal FEAT_DIM + OFI_DIM ({}) for legacy column", - FEAT_DIM + OFI_DIM - ); - } - - let bar_count = reader.bar_count(); - let mut feature_matrix: Vec = Vec::with_capacity(bar_count * in_dim); - let mut prices: Vec = Vec::with_capacity(bar_count); - for i in 0..bar_count { - let rec = reader.record(i); - // Labels always come from raw_close in the targets column (unchanged across feature versions). - prices.push(rec.targets[COL_RAW_CLOSE - FEAT_DIM]); - if use_alpha { - let alpha = reader - .alpha_features(i) - .ok_or_else(|| anyhow::anyhow!("alpha_feature_dim is Some but alpha_features({i}) is None"))?; - debug_assert_eq!(alpha.len(), in_dim); - feature_matrix.extend_from_slice(alpha); - } else { - // Legacy 74-dim layout: features (42) + OFI (32). - feature_matrix.extend_from_slice(rec.features); - feature_matrix.extend_from_slice(rec.ofi); - } - } - - let (n_nan, n_inf, fmin, fmax, fmean, fstd) = feature_stats(&feature_matrix); - tracing::info!( - feature_bytes = feature_matrix.len() * 4, - n_nan, n_inf, fmin, fmax, fmean, fstd, - "feature matrix extracted (CPU-resident)" - ); - - const FEATURE_MAGNITUDE_CAP: f32 = 1.0e6; - let mut valid_bar: Vec = Vec::with_capacity(bar_count); - let mut dropped_for_features = 0usize; - let mut first_bad_bars: Vec = Vec::new(); - for i in 0..bar_count { - let row = &feature_matrix[i * in_dim..(i + 1) * in_dim]; - let ok = row.iter().all(|&x| x.is_finite() && x.abs() < FEATURE_MAGNITUDE_CAP); - valid_bar.push(ok); - if !ok { - dropped_for_features += 1; - if first_bad_bars.len() < 10 { - first_bad_bars.push(i); - } - } - } - let col_stats = per_column_stats(&feature_matrix, bar_count, in_dim, FEATURE_MAGNITUDE_CAP); - let bad_cols: Vec<(usize, usize, usize, usize)> = col_stats - .iter() - .enumerate() - .filter(|(_, s)| s.n_nan + s.n_inf + s.n_extreme > 0) - .map(|(c, s)| (c, s.n_nan, s.n_inf, s.n_extreme)) - .collect(); - let decile_size = bar_count / 10; - let mut bad_by_decile = [0usize; 10]; - for &b in valid_bar - .iter() - .enumerate() - .filter(|(_, &ok)| !ok) - .map(|(i, _)| i) - .collect::>() - .iter() - { - let bucket = (b / decile_size.max(1)).min(9); - bad_by_decile[bucket] += 1; - } - tracing::info!( - dropped_for_features, - magnitude_cap = FEATURE_MAGNITUDE_CAP, - ?first_bad_bars, - n_bad_columns = bad_cols.len(), - bad_columns_summary = ?bad_cols, - ?bad_by_decile, - "feature corruption audit" - ); - - let horizon = config.horizon; - let (train_indices, train_labels) = - collect_valid_labels(&prices, &split.train, horizon, &valid_bar); - let (val_indices, val_labels) = - collect_valid_labels(&prices, &split.val, horizon, &valid_bar); - - let up_fraction_train = mean_label(&train_labels); - let up_fraction_val = mean_label(&val_labels); - let train_tied = split.n_train - train_indices.len(); - let val_tied = split.n_val - val_indices.len(); - - tracing::info!( - n_train = train_indices.len(), - n_val = val_indices.len(), - up_fraction_train, - up_fraction_val, - train_tied, - val_tied, - "label generation complete (raw_close direction labels)" - ); - - Ok(Phase1aData { - feature_matrix, - in_dim, - bar_count, - train_indices, - train_labels, - val_indices, - val_labels, - up_fraction_train, - up_fraction_val, - n_val_tied_or_invalid: val_tied, - }) -} - -/// Collect (bar_idx, 0/1 label) pairs over a range, dropping bars where the -/// label is tied OR the feature row was flagged invalid. -fn collect_valid_labels( - prices: &[f32], - range: &std::ops::Range, - horizon: usize, - valid_bar: &[bool], -) -> (Vec, Vec) { - let mut indices = Vec::with_capacity(range.len()); - let mut labels = Vec::with_capacity(range.len()); - for t in range.clone() { - if !valid_bar[t] { - continue; - } - if let Some(y) = binary_direction_label(prices, t, horizon) { - indices.push(t); - labels.push(y as f32); - } - } - (indices, labels) -} - -/// Per-column corruption stats: counts of NaN / Inf / extreme-magnitude -/// values in each feature dimension. Used to audit whether bad data is -/// concentrated in specific feature columns (e.g. an OFI dimension that -/// defaults to a sentinel) vs in specific bars (e.g. warmup periods). -#[derive(Debug, Clone, Copy)] -struct ColumnStats { - n_nan: usize, - n_inf: usize, - n_extreme: usize, -} - -fn per_column_stats(xs: &[f32], bar_count: usize, in_dim: usize, cap: f32) -> Vec { - let mut stats = vec![ColumnStats { n_nan: 0, n_inf: 0, n_extreme: 0 }; in_dim]; - for i in 0..bar_count { - let row = &xs[i * in_dim..(i + 1) * in_dim]; - for (c, &x) in row.iter().enumerate() { - if x.is_nan() { - stats[c].n_nan += 1; - } else if !x.is_finite() { - stats[c].n_inf += 1; - } else if x.abs() >= cap { - stats[c].n_extreme += 1; - } - } - } - stats -} - -/// Compute basic stats on the feature matrix for sanity checking. -fn feature_stats(xs: &[f32]) -> (usize, usize, f32, f32, f32, f32) { - let mut n_nan = 0usize; - let mut n_inf = 0usize; - let mut min = f32::INFINITY; - let mut max = f32::NEG_INFINITY; - let mut sum = 0.0_f64; - let mut sum_sq = 0.0_f64; - let mut n_finite = 0usize; - for &x in xs { - if x.is_nan() { - n_nan += 1; - } else if !x.is_finite() { - n_inf += 1; - } else { - if x < min { min = x; } - if x > max { max = x; } - sum += x as f64; - sum_sq += (x as f64) * (x as f64); - n_finite += 1; - } - } - let mean = if n_finite > 0 { (sum / n_finite as f64) as f32 } else { f32::NAN }; - let var = if n_finite > 1 { - let m = sum / n_finite as f64; - ((sum_sq / n_finite as f64) - m * m).max(0.0) - } else { - 0.0 - }; - let std = (var.sqrt()) as f32; - (n_nan, n_inf, min, max, mean, std) -} - -/// Per-feature z-score normalization statistics. -/// -/// Mirrors the production convention in `crates/ml/src/walk_forward.rs`: -/// - **f64 accumulators**: sums of f32 features can lose precision at 175k×74 -/// samples; f64 keeps it. -/// - **Two-pass mean/variance** (not Welford): non-streaming data so the -/// simpler algorithm is fine and matches the production reference exactly. -/// - **MIN_STD = 1e-8**: prevents division by zero on degenerate (constant) -/// feature columns. -/// - **FEATURE_BOUND = 20.0**: post-divide clip. Real z-scores live within ±5 -/// even on extreme bars; ±20 traps any latent corruption without rejecting -/// legitimate signal. -/// - **Train-only fit**: caller passes train_indices; val statistics never -/// enter the baseline (lookahead-bias audit rec 3). -#[derive(Debug, Clone)] -struct NormStats { - mean: Vec, - std: Vec, -} - -impl NormStats { - const MIN_STD: f64 = 1.0e-8; - const FEATURE_BOUND: f32 = 20.0; - - /// Fit (mean, std) per feature column using only `train_indices` rows of - /// `feature_matrix`. Returns degenerate (mean=0, std=1) stats if the - /// train set is empty. - fn fit_train(feature_matrix: &[f32], in_dim: usize, train_indices: &[usize]) -> Self { - if train_indices.is_empty() { - return Self { - mean: vec![0.0; in_dim], - std: vec![1.0; in_dim], - }; - } - let n = train_indices.len() as f64; - let mut mean = vec![0.0_f64; in_dim]; - for &t in train_indices { - let row = &feature_matrix[t * in_dim..(t + 1) * in_dim]; - for (m, &x) in mean.iter_mut().zip(row.iter()) { - *m += x as f64; - } - } - for m in &mut mean { - *m /= n; - } - let mut variance = vec![0.0_f64; in_dim]; - for &t in train_indices { - let row = &feature_matrix[t * in_dim..(t + 1) * in_dim]; - for (v, (&x, &m)) in variance.iter_mut().zip(row.iter().zip(mean.iter())) { - let diff = x as f64 - m; - *v += diff * diff; - } - } - let std: Vec = variance - .iter() - .map(|v| (v / n).sqrt().max(Self::MIN_STD)) - .collect(); - Self { mean, std } - } - - /// Apply z-score + clip in-place across the entire feature matrix. - /// Both train AND val rows are transformed using the train-fit stats. - fn apply_inplace(&self, feature_matrix: &mut [f32], in_dim: usize) { - debug_assert_eq!(self.mean.len(), in_dim); - debug_assert_eq!(self.std.len(), in_dim); - let bound = Self::FEATURE_BOUND; - for row in feature_matrix.chunks_exact_mut(in_dim) { - for (x, (&m, &s)) in row.iter_mut().zip(self.mean.iter().zip(self.std.iter())) { - let z = ((*x as f64 - m) / s) as f32; - *x = z.clamp(-bound, bound); - } - } - } -} - -fn mean_label(labels: &[f32]) -> f32 { - if labels.is_empty() { - return 0.5; - } - let sum: f64 = labels.iter().map(|&y| y as f64).sum(); - (sum / labels.len() as f64) as f32 -}