gitadmin/foxhunt
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

refactor(ml-alpha): gut Phase 1a sources for greenfield CfC rebuild

Browse Source 
Removes mlp/training/eval/backtest/metrics_detail/calibration and the
old example trainers. Preserves multi_horizon_labels, purged_split,
fxcache_reader (Phase A data path), mamba2_block (gate reference).
Subsequent commits populate cfc/heads/isv/pinned/trainer/data/gate.

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
This commit is contained in:
jgrusewski
2026-05-16 21:05:26 +02:00
parent 056b4abe52
commit 34739a53a9
12 changed files with 31 additions and 3409 deletions
Download Patch File Download Diff File Expand all files Collapse all files

129
crates/ml-alpha/examples/alpha_bar_baseline.rs
View File

@@ -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<String>,
/// 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<String>,
}
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(())
}

165
crates/ml-alpha/examples/alpha_bar_detailed.rs
View File

@@ -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<Vec<f32>> = (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(())
}

103
crates/ml-alpha/examples/alpha_calibrate.rs
View File

@@ -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<f32> {
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(())
}

729
crates/ml-alpha/examples/alpha_train_stacker.rs
View File

@@ -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<usize>,
/// 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<String>,
/// 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<usize>,
}
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<f32> = Vec::with_capacity(n_bars * alpha_dim);
let mut prices: Vec<f32> = 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::<f32>() / 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<usize> = (0..n_valid)
.filter(|&p| labels.valid_indices[p] < train_split_bar.saturating_sub(embargo))
.collect();
let val_label_range: Vec<usize> = (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<usize> {
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<f32>)> {
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<usize> = (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<usize> = 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<f32> = 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<f32> = Vec::with_capacity(val_pos.len());
let mut val_ys: Vec<f32> = 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<usize> = 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<u8> = 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::<f32>() / 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<u8> = 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<f32> {
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<f32> = 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<usize> = (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<f32> = 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<u8> =
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<f32> = 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<f32> = 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<f32> = 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::<f32>() / 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<f32> = stacker_logits.iter()
.map(|&z| 1.0_f32 / (1.0 + (-z.clamp(-50.0, 50.0)).exp()))
.collect();
let mut prices_t_host: Vec<f32> = Vec::with_capacity(stacker_test_n);
let mut prices_kt_host: Vec<f32> = 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<f32> = vec![0.00, 0.05, 0.10, 0.15, 0.20, 0.25];
let costs: Vec<f32> = 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<f32> = 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<f32> = (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<u8> =
val_ys.iter().map(|&y| if y > 0.5 { 1 } else { 0 }).collect();
let val_labels_f32: Vec<f32> = 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<f32> = 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(())
}

245
crates/ml-alpha/examples/gbm_baseline.rs
View File

@@ -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<String>,
/// 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<u8> = data
.val_labels
.iter()
.map(|&y| if y > 0.5 { 1 } else { 0 })
.collect();
let pred_logits: Vec<f32> = 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::<f32>() / 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(())
}

323
crates/ml-alpha/src/backtest.rs
View File

@@ -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<CudaStream>,
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<Self> {
// 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<f32>,
prices_t: &CudaSlice<f32>,
prices_kt: &CudaSlice<f32>,
thresholds: &[f32],
costs: &[f32],
test_time_span_s: f64,
) -> Result<Vec<BacktestStats>> {
let mut all_stats: Vec<BacktestStats> =
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<f32>,
prices_t: &CudaSlice<f32>,
prices_kt: &CudaSlice<f32>,
thresholds: &[f32],
cost_per_trade: f32,
test_time_span_s: f64,
) -> Result<Vec<BacktestStats>> {
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::<f32>(nt)
.map_err(|e| anyhow!("alloc pnl_grid: {e}"))?;
let mut trade_grid = self.stream
.alloc_zeros::<i32>(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::<f32>(t as usize)
.map_err(|e| anyhow!("alloc sum_ret: {e}"))?;
let mut sum_sq = self.stream.alloc_zeros::<f32>(t as usize)
.map_err(|e| anyhow!("alloc sum_sq: {e}"))?;
let mut sum_trades = self.stream.alloc_zeros::<i32>(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)
}
}

202
crates/ml-alpha/src/calibration.rs
View File

@@ -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<f32>;
}
/// 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<Self, &'static str> {
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<f32> {
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<f32>,
values: Vec<f32>,
}
impl IsotonicCalibrator {
pub fn fit(logits: &[f32], labels: &[f32]) -> Result<Self, &'static str> {
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<f32> = paired.iter().map(|p| p.0).collect();
let mut values: Vec<f32> = paired.iter().map(|p| p.1).collect();
let mut weights: Vec<f32> = 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<f32> {
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::<dyn Calibrator>();
}
#[test]
fn test_platt_perfect_separable() {
let logits: Vec<f32> = (0..200)
.map(|i| if i % 2 == 0 { 2.0 } else { -2.0 })
.collect();
let labels: Vec<f32> = (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<f32> = (0..100).map(|i| i as f32 * 0.1).collect();
let labels: Vec<f32> = (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);
}
}

210
crates/ml-alpha/src/eval.rs
View File

@@ -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<usize> = (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);
}
}

89
crates/ml-alpha/src/lib.rs
View File

@@ -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};

258
crates/ml-alpha/src/metrics_detail.rs
View File

@@ -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<ReliabilityBin> {
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<StratificationBin> {
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<f32> = feature_values.iter().copied().collect();
sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
let mut cuts: Vec<f32> = 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<f32> = (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<f32> = (-50..50).map(|i| i as f32 * 0.1).collect();
let labels: Vec<f32> = (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<f32> = (0..100).map(|i| (i % 2) as f32).collect();
let features: Vec<f32> = (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);
}
}
}

288
crates/ml-alpha/src/mlp.rs
View File

@@ -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<CudaStream>,
// 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<CudaStream>) -> Result<Self> {
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<f32>` 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<Vec<f32>> {
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<f32> {
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<String, GpuTensor> = 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);
}
}

699
crates/ml-alpha/src/training.rs
View File

@@ -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<MlpConfigSerde> 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<CudaStream>,
}
impl Phase1aTrainer {
/// Initialize trainer from config. Loads fxcache, computes purged split,
/// constructs MLP.
pub fn from_config(mut config: Phase1aConfig, stream: Arc<CudaStream>) -> Result<Self> {
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<EvalReport> {
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<Phase1aRunOutputs> {
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<f32> = vec![0.0; batch_buf_feats];
let mut batch_labels_buf: Vec<f32> = vec![0.0; batch_size];
let mut shuffled: Vec<usize> = (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<f32> = Vec::with_capacity(val_n);
let val_chunk = 8192usize.min(val_n);
let mut chunk_features: Vec<f32> = 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<u8> = 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<f32>,
/// Ground-truth binary labels for the same val samples (length = n_val).
pub val_labels: Vec<f32>,
/// 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<usize>,
}
/// 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<f32>,
/// 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<usize>,
/// Parallel labels in {0.0, 1.0} for `train_indices`.
pub train_labels: Vec<f32>,
/// Bar indices in the validation half whose features AND label are valid.
pub val_indices: Vec<usize>,
/// Parallel labels in {0.0, 1.0} for `val_indices`.
pub val_labels: Vec<f32>,
/// 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<Phase1aData> {
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<f32> = Vec::with_capacity(bar_count * in_dim);
let mut prices: Vec<f32> = 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<bool> = Vec::with_capacity(bar_count);
let mut dropped_for_features = 0usize;
let mut first_bad_bars: Vec<usize> = 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::<Vec<_>>()
.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<usize>,
horizon: usize,
valid_bar: &[bool],
) -> (Vec<usize>, Vec<f32>) {
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<ColumnStats> {
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<f64>,
std: Vec<f64>,
}
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<f64> = 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
}