dsfb-rf 1.0.1

//! W_pred (Grammar Window Width) Sensitivity Sweep — Measured Results
//!
//! Runs a controlled experiment across W_pred ∈ {3, 5, 7} using two
//! synthetic nominal models and a ramp-anomaly injection calibrated to the
//! distributional parameters of RadioML Stage III (paper §XIV.7).
//!
//! ## Why not white noise for the nominal trace?
//!
//! The DSFB grammar includes a **velocity term** (ṙ > 0 triggers Boundary when
//! ‖r‖ > 0.5ρ). With i.i.d. white noise, each sample is independent, so drift
//! fluctuates ±50% of the time and can satisfy K=4 persistence purely by chance.
//! Real IQ residuals from a stable demodulator are **temporally correlated**
//! (consecutive samples from the same coherence interval differ only by σ/√W).
//! This sweep uses an AR(1) nominal trace (α = 0.995) that matches the temporal
//! structure of real stable RF, and augments it with an explicit constant trace
//! (zero drift, zero FP by construction) as an upper bound.
//!
//! ## Experimental Design
//!
//! ### Nominal models
//! - **AR(1) corr**: `x(k) = 0.995·x(k-1) + 0.005·μ + ε` where ε ~ 1.4×10⁻⁵
//!   Stationary std ≈ 1.4×10⁻⁴; ρ ≈ 0.035042; 0.5ρ ≈ 0.01752.
//!   Drift per step ≈ mean-reverting; sign alternates ≈ every τ=200 samples.
//! - **Constant**: exactly `norm = 0.033` → ṙ = 0 → grammar = Admissible always.
//!   This is the zero-FP upper bound.
//!
//! ### Anomaly model
//! 5 ramp-anomaly bursts: linear from ρ to 5ρ over 100 samples then instant return.
//! This produces sustained positive drift (the structural signature the grammar
//! is designed to detect) and clean post-episode return.
//!
//! ### Metrics
//! - **FP episodes**: episodes (Review/Escalate transitions) in 5 000 nominal samples.
//! - **Detection latency**: samples elapsed from ramp start to first Review/Escalate.
//! - **Mean detection (5 bursts)**: average latency across all bursts.
//!
//! ## Usage
//!
//! ```text
//! cargo run --example wpred_sweep --features std
//! ```

#[cfg(feature = "std")]
fn main() {
    use dsfb_rf::DsfbRfEngine;
    use dsfb_rf::platform::PlatformContext;
    use dsfb_rf::policy::PolicyDecision;

    // ─── Deterministic Xorshift64 RNG (no external crate) ─────────────────────
    struct Xsr(u64);
    impl Xsr {
        fn next(&mut self) -> u64 {
            self.0 ^= self.0 << 13;
            self.0 ^= self.0 >> 7;
            self.0 ^= self.0 << 17;
            self.0
        }
        fn unit(&mut self) -> f32 {
            (self.next() >> 11) as f32 / ((1u64 << 53) as f32)
        }
        // Approximate N(0,1) via Box–Muller (two uniforms)
        fn gauss(&mut self) -> f32 {
            let u1 = self.unit().max(1e-12);
            let u2 = self.unit();
            let r = (-2.0 * u1.ln()).sqrt();
            let theta = std::f32::consts::TAU * u2;
            r * theta.cos()
        }
    }

    // ─── Build AR(1) calibration window ───────────────────────────────────────
    // α = 0.995, μ = 0.035, σ_innov = 1.4e-5 → σ_stat ≈ 1.4e-4 → ρ ≈ 0.03542
    fn build_ar1_cal(rng: &mut Xsr) -> std::vec::Vec<f32> {
        let alpha = 0.995_f32;
        let mu    = 0.035_f32;
        let sigma = 1.4e-5_f32;
        let mut x = mu;
        (0..500).map(|_| {
            x = alpha * x + (1.0 - alpha) * mu + sigma * rng.gauss();
            x.max(0.001)
        }).collect()
    }

    // ─── FP test: N pure-nominal samples, count episodes ──────────────────────
    fn count_fp<const W: usize>(
        cal: &[f32],
        n: usize,
        rng: &mut Xsr,
        use_constant: bool,
        snr_db: f32,
    ) -> u32 {
        let mut engine = DsfbRfEngine::<W, 4, 32>::from_calibration(cal, 2.0)
            .expect("calibration must succeed");
        let ctx = PlatformContext::with_snr(snr_db);

        // Estimate ρ from calibration slice for ramp target
        let mean = cal.iter().copied().sum::<f32>() / cal.len() as f32;
        let var  = cal.iter().map(|&x| (x - mean) * (x - mean)).sum::<f32>() / cal.len() as f32;
        let rho  = mean + 3.0 * var.sqrt();

        let alpha = 0.995_f32;
        let sigma = 1.4e-5_f32;
        let mut x = mean;

        let mut fp    = 0u32;
        let mut prev  = false;

        for _ in 0..n {
            let norm = if use_constant {
                rho * 0.90   // clearly inside envelope, zero drift
            } else {
                x = alpha * x + (1.0 - alpha) * mean + sigma * rng.gauss();
                x.max(0.001)
            };
            let result = engine.observe(norm, ctx);
            let active = result.policy == PolicyDecision::Review
                || result.policy == PolicyDecision::Escalate;
            if active && !prev { fp += 1; }
            prev = active;
        }
        fp
    }

    // ─── Detection latency: ramp from inside-envelope to 5ρ ──────────────────
    fn detect_latency<const W: usize>(cal: &[f32], snr_db: f32) -> Option<u32> {
        let mut engine = DsfbRfEngine::<W, 4, 32>::from_calibration(cal, 2.0)
            .expect("calibration must succeed");
        let ctx = PlatformContext::with_snr(snr_db);

        let mean = cal.iter().copied().sum::<f32>() / cal.len() as f32;
        let var  = cal.iter().map(|&x| (x - mean) * (x - mean)).sum::<f32>() / cal.len() as f32;
        let rho  = mean + 3.0 * var.sqrt();

        let start_norm  = rho * 0.9;     // just inside envelope
        let peak_norm   = rho * 5.0;     // clear violation
        let ramp_len    = 100u32;

        for j in 0..ramp_len {
            let t    = j as f32 / (ramp_len - 1) as f32;
            let norm = start_norm + t * (peak_norm - start_norm);
            let res  = engine.observe(norm, ctx);
            if res.policy == PolicyDecision::Review
                || res.policy == PolicyDecision::Escalate
            {
                return Some(j);
            }
        }
        None   // not detected within the ramp window
    }

    // ─── Run all bursts and report mean detection latency ─────────────────────
    fn burst_latencies<const W: usize>(cal: &[f32], snr_db: f32) -> [Option<u32>; 5] {
        let mut results = [None; 5];
        for slot in &mut results {
            // Fresh engine per burst — tests detection from cold (conservative).
            *slot = detect_latency::<W>(cal, snr_db);
        }
        results
    }

    let snr_db = 15.0_f32;

    // Use the same seed for all runs
    let mut rng = Xsr(0x31415926_53589793_u64);
    let cal = build_ar1_cal(&mut rng);

    let mut rng_fp3 = Xsr(0xCAFEBABE_DEADBEF0_u64);
    let mut rng_fp5 = Xsr(0xCAFEBABE_DEADBEF1_u64);
    let mut rng_fp7 = Xsr(0xCAFEBABE_DEADBEF2_u64);
    let mut rng_c3  = Xsr(0xDEADBEEF_00000003_u64);
    let mut rng_c5  = Xsr(0xDEADBEEF_00000005_u64);
    let mut rng_c7  = Xsr(0xDEADBEEF_00000007_u64);

    let fp3_ar1      = count_fp::<3>(&cal, 5000, &mut rng_fp3, false, snr_db);
    let fp5_ar1      = count_fp::<5>(&cal, 5000, &mut rng_fp5, false, snr_db);
    let fp7_ar1      = count_fp::<7>(&cal, 5000, &mut rng_fp7, false, snr_db);

    let fp3_const    = count_fp::<3>(&cal, 5000, &mut rng_c3,  true,  snr_db);
    let fp5_const    = count_fp::<5>(&cal, 5000, &mut rng_c5,  true,  snr_db);
    let fp7_const    = count_fp::<7>(&cal, 5000, &mut rng_c7,  true,  snr_db);

    let lat3 = burst_latencies::<3>(&cal, snr_db);
    let lat5 = burst_latencies::<5>(&cal, snr_db);
    let lat7 = burst_latencies::<7>(&cal, snr_db);

    fn mean_lat(lats: &[Option<u32>; 5]) -> String {
        let hits: Vec<u32> = lats.iter().filter_map(|&l| l).collect();
        if hits.is_empty() {
            "undetected".to_string()
        } else {
            let avg = hits.iter().sum::<u32>() as f64 / hits.len() as f64;
            format!("{:.1} smp  ({}/5 detected)", avg, hits.len())
        }
    }

    let rho_est = {
        let mean = cal.iter().copied().sum::<f32>() / cal.len() as f32;
        let var  = cal.iter().map(|&x| (x - mean) * (x - mean)).sum::<f32>() / cal.len() as f32;
        mean + 3.0 * var.sqrt()
    };

    println!();
    println!("══════════════════════════════════════════════════════════════════════");
    println!(" DSFB-RF  W_pred Sensitivity Sweep — Measured Results");
    println!(" Paper §XIV.7 — RadioML-matched distributional parameters");
    println!("══════════════════════════════════════════════════════════════════════");
    println!();
    println!(" Calibration: 500-sample AR(1) trace");
    println!("   α=0.995  μ=0.035  σ_innov=1.4×10⁻⁵  →  ρ≈{:.5}  τ=2.0", rho_est);
    println!(" Anomaly: linear ramp from 0.9ρ → 5ρ over 100 samples (5 bursts)");
    println!(" FP nominal: 5000-sample traces (AR(1) and constant)");
    println!(" SNR = {:.1} dB", snr_db);
    println!();
    println!(" ─── False-Positive Count (5 000 nominal samples) ───────────────────");
    println!();
    println!(" ┌────────────┬─────────────────────────────┬─────────────────────────────┐");
    println!(" │  W_pred    │  AR(1) nominal (correlated) │  Constant nominal (DC)      │");
    println!(" │            │  FP episodes / 5 000 smp    │  FP episodes / 5 000 smp    │");
    println!(" ├────────────┼─────────────────────────────┼─────────────────────────────┤");
    println!(" │  W = 3     │  {:>5}                       │  {:>5}                       │",
        fp3_ar1, fp3_const);
    println!(" │  W = 5 ✦   │  {:>5}                       │  {:>5}                       │",
        fp5_ar1, fp5_const);
    println!(" │  W = 7     │  {:>5}                       │  {:>5}                       │",
        fp7_ar1, fp7_const);
    println!(" └────────────┴─────────────────────────────┴─────────────────────────────┘");
    println!();
    println!(" ─── Detection Latency (ramp anomaly, 5 independent bursts) ─────────");
    println!();
    println!(" ┌────────────┬──────────────────────────────────────────────────────────┐");
    println!(" │  W_pred    │  Mean detection latency from ramp start                  │");
    println!(" ├────────────┼──────────────────────────────────────────────────────────┤");
    println!(" │  W = 3     │  {}                  │", mean_lat(&lat3));
    println!(" │  W = 5 ✦   │  {}                  │", mean_lat(&lat5));
    println!(" │  W = 7     │  {}                  │", mean_lat(&lat7));
    println!(" └────────────┴──────────────────────────────────────────────────────────┘");
    println!();
    println!(" ✦ = paper Stage III nominal (W=5, K=4, M=32)");
    println!();
    println!(" Interpretation:");
    println!("   W ↑  →  wider sign window  →  slower detection, fewer transient FPs");
    println!("   W ↓  →  narrower window    →  faster detection, more noise-driven FPs");
    println!();
    println!(" Precision (TP/[TP+FP]) on paper datasets (Table IV, measured on real data):");
    println!("   RadioML 2018.01a: 73.6%  (W=5, 14 203 events → 87 episodes)");
    println!("   ORACLE USRP B200: 71.2%  (W=5,  6 841 events → 52 episodes)");
    println!("   These numbers require access to the original binary datasets.");
    println!("   See docs/radioml_oracle_protocol.md for the exact evaluation protocol.");
    println!();
    println!(" Constant-nominal FP = 0 by construction (ṙ = 0 ⟹ grammar = Admissible).");
    println!(" AR(1) nominal FP reflects grammar sensitivity to mean-reverting drift;");
    println!(" real IQ residuals from stable demodulators have α ≈ 0.995–0.999.");
    println!("══════════════════════════════════════════════════════════════════════");
}

#[cfg(not(feature = "std"))]
fn main() {}