foxhunt/ml/tests/cusum_test.rs

//! CUSUM Structural Break Detector Tests
//!
//! Comprehensive test suite for two-sided CUSUM algorithm following TDD methodology.
//! Tests cover:
//! - Algorithm correctness (mean shifts, threshold sensitivity)
//! - Real market data integration (ES.FUT, NQ.FUT, ZN.FUT, 6E.FUT)
//! - Performance benchmarks (<50μs per update)
//! - Statistical validation (false positive rate, detection delay)
//! - Property-based testing (invariants, edge cases)

use approx::assert_relative_eq;
use ml::regime::cusum::{CUSUMDetector, StructuralBreak};
use proptest::prelude::*;
use statrs::distribution::{ContinuousCDF, Normal};
use std::time::Instant;

// ===== Basic Functionality Tests =====

#[test]
fn test_cusum_no_change_stable() {
    // Stable data with no mean shift should not trigger detection
    let mut detector = CUSUMDetector::new(0.0, 1.0, 0.5, 5.0);

    // Generate 1000 samples from N(0, 1)
    let mut rng = fastrand::Rng::with_seed(42);
    for _ in 0..1000 {
        let value = rng.f64() * 2.0 - 1.0; // Approx uniform [-1, 1]
        let result = detector.update(value);

        // Should not detect breaks in stable data
        assert!(result.is_none(), "False positive on stable data");
    }

    // CUSUM sums should remain bounded
    let (s_pos, s_neg) = detector.get_current_sums();
    assert!(s_pos < 10.0, "Positive CUSUM unbounded: {}", s_pos);
    assert!(s_neg < 10.0, "Negative CUSUM unbounded: {}", s_neg);
}

#[test]
fn test_cusum_mean_increase() {
    // Detect positive mean shift from 0 to +2σ
    let mut detector = CUSUMDetector::new(0.0, 1.0, 0.5, 4.0);

    // First 50 samples from N(0, 1)
    let mut rng = fastrand::Rng::with_seed(42);
    for _ in 0..50 {
        let value = rng.f64() * 2.0 - 1.0;
        let result = detector.update(value);
        assert!(result.is_none(), "Premature detection");
    }

    // Next 50 samples from N(+2, 1) - mean shift
    let mut break_detected = false;
    for _ in 0..50 {
        let value = (rng.f64() * 2.0 - 1.0) + 2.0; // Shift by +2σ
        if let Some(structural_break) = detector.update(value) {
            assert_eq!(structural_break.direction, "positive");
            assert!(structural_break.magnitude > 0.0);
            break_detected = true;
            break;
        }
    }

    assert!(break_detected, "Failed to detect positive mean shift");
}

#[test]
fn test_cusum_mean_decrease() {
    // Detect negative mean shift from 0 to -2σ
    let mut detector = CUSUMDetector::new(0.0, 1.0, 0.5, 4.0);

    // First 50 samples from N(0, 1)
    let mut rng = fastrand::Rng::with_seed(43);
    for _ in 0..50 {
        let value = rng.f64() * 2.0 - 1.0;
        let result = detector.update(value);
        assert!(result.is_none(), "Premature detection");
    }

    // Next 50 samples from N(-2, 1) - mean shift
    let mut break_detected = false;
    for _ in 0..50 {
        let value = (rng.f64() * 2.0 - 1.0) - 2.0; // Shift by -2σ
        if let Some(structural_break) = detector.update(value) {
            assert_eq!(structural_break.direction, "negative");
            assert!(structural_break.magnitude < 0.0);
            break_detected = true;
            break;
        }
    }

    assert!(break_detected, "Failed to detect negative mean shift");
}

#[test]
fn test_cusum_threshold_sensitivity() {
    // Lower threshold (h) should detect sooner
    let mut detector_low = CUSUMDetector::new(0.0, 1.0, 0.5, 3.0); // h=3
    let mut detector_high = CUSUMDetector::new(0.0, 1.0, 0.5, 5.0); // h=5

    let mut rng = fastrand::Rng::with_seed(44);
    let mut low_detected_at = None;
    let mut high_detected_at = None;

    for i in 0..100 {
        let value = (rng.f64() * 2.0 - 1.0) + 1.5; // Moderate shift +1.5σ

        if detector_low.update(value).is_some() && low_detected_at.is_none() {
            low_detected_at = Some(i);
        }

        if detector_high.update(value).is_some() && high_detected_at.is_none() {
            high_detected_at = Some(i);
        }
    }

    // Lower threshold should detect earlier (or at all)
    assert!(low_detected_at.is_some(), "Low threshold failed to detect");

    if let (Some(low), Some(high)) = (low_detected_at, high_detected_at) {
        assert!(low <= high, "Lower threshold should detect earlier");
    }
}

#[test]
fn test_cusum_drift_allowance() {
    // Higher drift allowance (k) makes detection more conservative
    let mut detector_low_k = CUSUMDetector::new(0.0, 1.0, 0.25, 4.0); // k=0.25
    let mut detector_high_k = CUSUMDetector::new(0.0, 1.0, 1.0, 4.0); // k=1.0

    let mut rng = fastrand::Rng::with_seed(45);
    let mut low_k_detected = false;
    let mut high_k_detected = false;

    for _ in 0..100 {
        let value = (rng.f64() * 2.0 - 1.0) + 1.2; // Small shift +1.2σ

        if detector_low_k.update(value).is_some() {
            low_k_detected = true;
        }

        if detector_high_k.update(value).is_some() {
            high_k_detected = true;
        }
    }

    // Lower k should be more sensitive (more likely to detect small shifts)
    assert!(low_k_detected, "Low k should detect small shifts");
}

#[test]
fn test_cusum_reset_after_detection() {
    // After detection, reset should clear CUSUM sums
    let mut detector = CUSUMDetector::new(0.0, 1.0, 0.5, 4.0);

    // Trigger detection
    let mut rng = fastrand::Rng::with_seed(46);
    for _ in 0..100 {
        let value = (rng.f64() * 2.0 - 1.0) + 2.0;
        if detector.update(value).is_some() {
            break;
        }
    }

    // Check sums before reset (should be high)
    let (s_pos_before, s_neg_before) = detector.get_current_sums();
    assert!(
        s_pos_before > 0.0 || s_neg_before > 0.0,
        "No CUSUM accumulation"
    );

    // Reset
    detector.reset();

    // Check sums after reset (should be zero)
    let (s_pos_after, s_neg_after) = detector.get_current_sums();
    assert_relative_eq!(s_pos_after, 0.0, epsilon = 1e-10);
    assert_relative_eq!(s_neg_after, 0.0, epsilon = 1e-10);
}

#[test]
fn test_cusum_false_positive_rate() {
    // False positive rate should be <5% on pure Gaussian noise
    let num_trials = 100;
    let samples_per_trial = 500;
    let mut false_positives = 0;

    for trial in 0..num_trials {
        let mut detector = CUSUMDetector::new(0.0, 1.0, 0.5, 5.0);
        let mut rng = fastrand::Rng::with_seed(100 + trial);

        let mut detected = false;
        for _ in 0..samples_per_trial {
            let value = rng.f64() * 2.0 - 1.0; // Uniform approx Gaussian
            if detector.update(value).is_some() {
                detected = true;
                break;
            }
        }

        if detected {
            false_positives += 1;
        }
    }

    let fpr = false_positives as f64 / num_trials as f64;
    assert!(
        fpr < 0.05,
        "False positive rate too high: {:.2}%",
        fpr * 100.0
    );
}

#[test]
fn test_cusum_detection_delay() {
    // Detection should occur within 5-10 bars after a 2σ shift
    let mut detector = CUSUMDetector::new(0.0, 1.0, 0.5, 4.0);

    // Stable period
    let mut rng = fastrand::Rng::with_seed(47);
    for _ in 0..50 {
        let value = rng.f64() * 2.0 - 1.0;
        detector.update(value);
    }

    // Mean shift and measure delay
    let mut detection_delay = None;
    for i in 0..20 {
        let value = (rng.f64() * 2.0 - 1.0) + 2.5; // Strong shift +2.5σ
        if detector.update(value).is_some() {
            detection_delay = Some(i);
            break;
        }
    }

    assert!(detection_delay.is_some(), "Failed to detect within 20 bars");
    let delay = detection_delay.unwrap();
    assert!(delay < 10, "Detection delay too high: {} bars", delay);
}

// ===== Performance Benchmarks =====

#[test]
fn test_cusum_performance_sub_50us() {
    // Each update should complete in <50μs
    let mut detector = CUSUMDetector::new(0.0, 1.0, 0.5, 5.0);
    let mut rng = fastrand::Rng::with_seed(48);

    let num_updates = 10_000;
    let start = Instant::now();

    for _ in 0..num_updates {
        let value = rng.f64() * 2.0 - 1.0;
        detector.update(value);
    }

    let elapsed = start.elapsed();
    let avg_latency_us = elapsed.as_micros() as f64 / num_updates as f64;

    println!("Average CUSUM update latency: {:.2}μs", avg_latency_us);
    assert!(
        avg_latency_us < 50.0,
        "Performance target not met: {:.2}μs",
        avg_latency_us
    );
}

// ===== Real Market Data Integration Tests =====

#[cfg(test)]
mod real_data_tests {
    use super::*;
    use dbn::decode::dbn::Decoder;
    use dbn::decode::DecodeRecord;
    use std::fs::File;
    use std::io::BufReader;

    fn load_dbn_file(path: &str) -> Vec<f64> {
        let file = File::open(path).expect("Failed to open DBN file");
        let reader = BufReader::new(file);
        let mut decoder = Decoder::new(reader).expect("Failed to create decoder");

        let mut prices = Vec::new();
        while let Some(record) = decoder
            .decode_record::<dbn::OhlcvMsg>()
            .expect("Failed to decode record")
        {
            // Use close price, convert from fixed-point (divide by 1e9)
            let close_price = record.close as f64 / 1_000_000_000.0;
            prices.push(close_price);
        }

        prices
    }

    fn compute_returns(prices: &[f64]) -> Vec<f64> {
        prices.windows(2).map(|w| (w[1] - w[0]) / w[0]).collect()
    }

    #[test]
    fn test_cusum_es_fut_real_data() {
        // ES.FUT (E-mini S&P 500) - test on real market data
        let path =
            "/home/jgrusewski/Work/foxhunt/test_data/real/databento/ES.FUT_ohlcv-1m_2024-01-02.dbn";

        if !std::path::Path::new(path).exists() {
            println!("Skipping ES.FUT test - file not found: {}", path);
            return;
        }

        let prices = load_dbn_file(path);
        assert!(!prices.is_empty(), "No data loaded from ES.FUT");
        println!("Loaded {} bars from ES.FUT", prices.len());

        // Compute returns
        let returns = compute_returns(&prices);

        // Estimate mean and std from first 100 bars
        let calibration_data = &returns[..100.min(returns.len())];
        let mean: f64 = calibration_data.iter().sum::<f64>() / calibration_data.len() as f64;
        let variance: f64 = calibration_data
            .iter()
            .map(|x| (x - mean).powi(2))
            .sum::<f64>()
            / calibration_data.len() as f64;
        let std_dev = variance.sqrt();

        println!(
            "ES.FUT return stats - mean: {:.6}, std: {:.6}",
            mean, std_dev
        );

        // Create detector
        let mut detector = CUSUMDetector::new(mean, std_dev, 0.5, 4.5);

        // Process remaining returns
        let mut breaks_detected = Vec::new();
        for (i, &ret) in returns.iter().enumerate().skip(100) {
            if let Some(structural_break) = detector.update(ret) {
                breaks_detected.push((i, structural_break));
                detector.reset(); // Reset after detection
            }
        }

        println!(
            "Detected {} structural breaks in ES.FUT",
            breaks_detected.len()
        );

        // Should detect at least some breaks in real market data
        assert!(
            breaks_detected.len() > 0,
            "Expected some structural breaks in real data"
        );
        assert!(
            breaks_detected.len() < returns.len() / 10,
            "Too many breaks detected"
        );
    }

    #[test]
    fn test_cusum_6e_fut_real_data() {
        // 6E.FUT (Euro FX) - test on currency futures
        let path = "/home/jgrusewski/Work/foxhunt/test_data/real/databento/ml_training_small/6E.FUT_ohlcv-1m_2024-01-02.dbn";

        if !std::path::Path::new(path).exists() {
            println!("Skipping 6E.FUT test - file not found: {}", path);
            return;
        }

        let prices = load_dbn_file(path);
        assert!(!prices.is_empty(), "No data loaded from 6E.FUT");
        println!("Loaded {} bars from 6E.FUT", prices.len());

        let returns = compute_returns(&prices);

        // Currency markets typically have different characteristics
        let calibration_data = &returns[..50.min(returns.len())];
        let mean: f64 = calibration_data.iter().sum::<f64>() / calibration_data.len() as f64;
        let variance: f64 = calibration_data
            .iter()
            .map(|x| (x - mean).powi(2))
            .sum::<f64>()
            / calibration_data.len() as f64;
        let std_dev = variance.sqrt();

        println!(
            "6E.FUT return stats - mean: {:.6}, std: {:.6}",
            mean, std_dev
        );

        let mut detector = CUSUMDetector::new(mean, std_dev, 0.5, 5.0);

        let mut breaks_detected = 0;
        for &ret in returns.iter().skip(50) {
            if detector.update(ret).is_some() {
                breaks_detected += 1;
                detector.reset();
            }
        }

        println!("Detected {} structural breaks in 6E.FUT", breaks_detected);
        assert!(
            breaks_detected < returns.len() / 5,
            "Too many breaks in currency data"
        );
    }

    #[test]
    fn test_cusum_multi_symbol_comparison() {
        // Compare break characteristics across different asset classes
        let symbols = vec![
            ("ES.FUT", "/home/jgrusewski/Work/foxhunt/test_data/real/databento/ES.FUT_ohlcv-1m_2024-01-02.dbn"),
            ("6E.FUT", "/home/jgrusewski/Work/foxhunt/test_data/real/databento/ml_training_small/6E.FUT_ohlcv-1m_2024-01-02.dbn"),
        ];

        for (symbol, path) in symbols {
            if !std::path::Path::new(path).exists() {
                println!("Skipping {} - file not found", symbol);
                continue;
            }

            let prices = load_dbn_file(path);
            let returns = compute_returns(&prices);

            if returns.len() < 100 {
                continue;
            }

            // Calibrate
            let calib = &returns[..50];
            let mean: f64 = calib.iter().sum::<f64>() / calib.len() as f64;
            let var: f64 =
                calib.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / calib.len() as f64;
            let std_dev = var.sqrt();

            let mut detector = CUSUMDetector::new(mean, std_dev, 0.5, 4.5);

            let mut positive_breaks = 0;
            let mut negative_breaks = 0;

            for &ret in returns.iter().skip(50) {
                if let Some(sb) = detector.update(ret) {
                    if sb.direction == "positive" {
                        positive_breaks += 1;
                    } else {
                        negative_breaks += 1;
                    }
                    detector.reset();
                }
            }

            println!(
                "{} - Positive: {}, Negative: {}",
                symbol, positive_breaks, negative_breaks
            );
        }
    }

    #[test]
    fn test_cusum_es_fut_integration_break_rate() {
        // ES.FUT integration test validating 5.5% break rate
        // Uses real market data from 2024-01-08 to validate CUSUM detection parameters
        let path = "/home/jgrusewski/Work/foxhunt/test_data/real/databento/ml_training/ES.FUT_ohlcv-1m_2024-01-08.dbn";

        if !std::path::Path::new(path).exists() {
            println!(
                "Skipping ES.FUT integration test - file not found: {}",
                path
            );
            println!("Expected path: {}", path);
            return;
        }

        // Load DBN data
        let prices = load_dbn_file(path);
        assert!(!prices.is_empty(), "No data loaded from ES.FUT");
        println!("Loaded {} bars from ES.FUT (2024-01-08)", prices.len());

        // Compute returns
        let returns = compute_returns(&prices);
        let total_bars = returns.len();
        println!("Total return bars: {}", total_bars);

        // Estimate mean and std from first 100 bars for calibration
        let calibration_size = 100.min(returns.len());
        let calibration_data = &returns[..calibration_size];

        let mean: f64 = calibration_data.iter().sum::<f64>() / calibration_data.len() as f64;
        let variance: f64 = calibration_data
            .iter()
            .map(|x| (x - mean).powi(2))
            .sum::<f64>()
            / calibration_data.len() as f64;
        let std_dev = variance.sqrt();

        println!("ES.FUT return statistics:");
        println!("  Mean: {:.6}", mean);
        println!("  Std Dev: {:.6}", std_dev);
        println!("  Calibration bars: {}", calibration_size);

        // Create CUSUM detector with standard parameters
        // k=0.5 (drift allowance), h=5.0 (detection threshold)
        let mut detector = CUSUMDetector::new(mean, std_dev, 0.5, 5.0);

        // Process remaining returns after calibration period
        let mut break_count = 0;
        let test_data = &returns[calibration_size..];

        for &ret in test_data.iter() {
            if let Some(structural_break) = detector.update(ret) {
                break_count += 1;
                println!(
                    "Break #{}: {} (magnitude: {:.3})",
                    break_count, structural_break.direction, structural_break.magnitude
                );
                detector.reset(); // Reset after detection
            }
        }

        // Calculate break rate as percentage of bars
        let break_rate = (break_count as f64 / test_data.len() as f64) * 100.0;

        println!("\nES.FUT Break Rate Analysis:");
        println!("  Breaks detected: {}", break_count);
        println!("  Test bars: {}", test_data.len());
        println!("  Break rate: {:.2}%", break_rate);

        // Validate break rate is within expected range [4.5%, 6.5%]
        // This validates that CUSUM parameters are properly tuned for ES.FUT
        assert!(
            break_rate >= 4.5 && break_rate <= 6.5,
            "ES.FUT break rate {:.2}% outside expected range [4.5%, 6.5%]. \
             Expected ~5.5% break rate for properly calibrated CUSUM detector.",
            break_rate
        );

        // Additional validation: ensure at least some breaks were detected
        assert!(
            break_count > 0,
            "No structural breaks detected - detector may be miscalibrated"
        );

        println!(
            "✓ ES.FUT integration test passed: break rate {:.2}% within [4.5%, 6.5%]",
            break_rate
        );
    }
}

// ===== Property-Based Tests =====

proptest! {
    #[test]
    fn test_cusum_invariant_nonnegative_sums(
        values in prop::collection::vec(-10.0..10.0f64, 10..100)
    ) {
        let mut detector = CUSUMDetector::new(0.0, 1.0, 0.5, 5.0);

        for value in values {
            detector.update(value);
            let (s_pos, s_neg) = detector.get_current_sums();

            // CUSUM sums must always be non-negative
            assert!(s_pos >= 0.0, "Positive CUSUM became negative: {}", s_pos);
            assert!(s_neg >= 0.0, "Negative CUSUM became negative: {}", s_neg);
        }
    }

    #[test]
    fn test_cusum_invariant_reset_clears_state(
        values in prop::collection::vec(-5.0..5.0f64, 10..50)
    ) {
        let mut detector = CUSUMDetector::new(0.0, 1.0, 0.5, 5.0);

        // Accumulate some state
        for value in values {
            detector.update(value);
        }

        // Reset
        detector.reset();

        // Verify zero state
        let (s_pos, s_neg) = detector.get_current_sums();
        assert_relative_eq!(s_pos, 0.0, epsilon = 1e-10);
        assert_relative_eq!(s_neg, 0.0, epsilon = 1e-10);
    }

    #[test]
    fn test_cusum_invariant_magnitude_bounds(
        shift in -5.0..5.0f64,
        std_dev in 0.1..2.0f64
    ) {
        let threshold = 4.0;
        let mut detector = CUSUMDetector::new(0.0, std_dev, 0.5, threshold);

        // Apply shift
        for _ in 0..100 {
            if let Some(sb) = detector.update(shift) {
                // When a break is detected, magnitude should exceed threshold
                assert!(sb.magnitude.abs() > threshold,
                        "Magnitude {} should exceed threshold {}", sb.magnitude.abs(), threshold);

                // Direction should match sign of shift
                if shift > 0.0 {
                    assert_eq!(sb.direction, "positive");
                    assert!(sb.magnitude > 0.0);
                } else if shift < 0.0 {
                    assert_eq!(sb.direction, "negative");
                    assert!(sb.magnitude < 0.0);
                }
                break;
            }
        }
    }
}

// ===== Edge Cases =====

#[test]
fn test_cusum_extreme_values() {
    let mut detector = CUSUMDetector::new(0.0, 1.0, 0.5, 5.0);

    // Should handle extreme values without panicking
    detector.update(f64::MAX / 1e10);
    detector.update(f64::MIN / 1e10);
    detector.update(0.0);

    let (s_pos, s_neg) = detector.get_current_sums();
    assert!(s_pos.is_finite());
    assert!(s_neg.is_finite());
}

#[test]
fn test_cusum_zero_variance() {
    // Zero variance should be handled gracefully
    let mut detector = CUSUMDetector::new(0.0, 0.0, 0.5, 5.0);

    // Should not panic
    detector.update(1.0);
    detector.update(2.0);

    let (s_pos, s_neg) = detector.get_current_sums();
    assert!(s_pos.is_finite());
    assert!(s_neg.is_finite());
}