feat(ml-alpha): GPU-native stacked regime head on top of Mamba2 (Phase 1d.3)

Builds a second-stage MLP that takes [mamba_logit, 6 Block-S features]
as input (7 dims) and learns the joint alpha-and-regime score in one
calibrated output. Trains on the cal half of val (same 50/50 split as
Platt/isotonic so all comparisons are on the same held-out test bars).

Architecture: 7 → hidden_dim → 1 sigmoid, GELU activation, BCE-with-logits
loss, AdamW. Uses the existing GPU-native `MlpModel` from
crates/ml-alpha/src/mlp.rs — same primitives used for the Phase 1c MLP
baseline. No CPU compute on the hot path (per feedback_cpu_is_read_only);
all weights, activations, gradients, optimizer state on GPU; host writes
the input matrix to a pinned buffer once per batch via GpuTensor::from_host.

The Block-S columns are z-score normalised using cal-half statistics
(then applied to the full val matrix) before training; mamba_logit is
left raw since it's already close to standard-normal scale via the
Mamba's natural calibration (see pearl_mamba_sss_state_yields_native_calibration).

After training, reports stacker held-out accuracy + AUC + Brier +
log-loss, plus stratified accuracy by Block-S feature so we can see
whether the stacker absorbed the regime conditioning (uniform accuracy
across quintiles) or just sharpened the Q4-gate (still elevated in Q4).

Why this matters for production deployment per pearl_mamba_inherits_regime_structure:
- Single calibrated score for conformal coverage gating downstream
- Retrainable when market regimes drift
- Captures interactions between regime features that a static threshold
  AND can't (e.g., spread-Q4 only when book is balanced)
- Replaces the planned Phase 1d.3 dual-head architecture with a smaller
  stacked-generalisation approach (no separate regime classifier)

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

crates/ml-alpha/examples/phase1d_long_horizon.rs

@@ -37,6 +37,7 @@ 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::multi_horizon_labels::generate_labels;
 
 use ml_core::cuda_autograd::gpu_tensor::GpuTensor;
@@ -64,6 +65,16 @@ struct Cli {
     /// 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,
+    /// Train a stacked regime head on the cal half? `0` skips.
+    #[arg(long, default_value_t = true)] stacker: bool,
+    /// 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,
 }
 
 fn main() -> Result<()> {
@@ -303,6 +314,160 @@ fn main() -> Result<()> {
                  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.
+        if cli.stacker {
+            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);
+
+            // 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).
+            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) ──