foxhunt/ml/tests/dqn_feature_quality_validation_test.rs

/// Feature Quality Validation Test
///
/// This test systematically validates all 225 features to identify potential causes
/// of gradient explosions during DQN training. Based on expert analysis, we check for:
///
/// 1. **NaN/Inf Values**: No invalid floating point numbers
/// 2. **Extreme Values**: Features exceeding reasonable bounds (e.g., >100σ)
/// 3. **Constant Features**: Features with zero variance (not learnable)
/// 4. **Sparse Features**: Features that are zero >95% of the time
/// 5. **Multicollinearity**: Highly correlated feature pairs (>0.95)
/// 6. **Statistical Stability**: Higher-order moments (skewness, kurtosis) are stable
/// 7. **Microstructure Stability**: Ratio-based features don't explode
/// 8. **Time-Series Stationarity**: Features don't exhibit extreme non-stationarity
///
/// **Expert Guidance**: 225 features is excessive for DQN. Successful implementations
/// use 20-60 features. Primary suspects for gradient explosions:
/// - Statistical features (skewness, kurtosis) - notoriously unstable
/// - Microstructure proxies (Amihud illiquidity) - can approach infinity
/// - High multicollinearity (multiple moving average ratios)
///
/// **Target**: Identify which features to remove/fix for stable training.
use anyhow::Result;
use ml::features::extraction::{extract_ml_features, OHLCVBar};
use ml::real_data_loader::RealDataLoader;
use std::collections::HashMap;

/// Statistical analysis of a feature vector
#[derive(Debug, Clone)]
struct FeatureStats {
    index: usize,
    mean: f64,
    std_dev: f64,
    min: f64,
    max: f64,
    nan_count: usize,
    inf_count: usize,
    zero_count: usize,
    total_count: usize,
    /// Percentage of values that are zero
    zero_ratio: f64,
    /// Range (max - min)
    range: f64,
}

impl FeatureStats {
    fn new(index: usize) -> Self {
        Self {
            index,
            mean: 0.0,
            std_dev: 0.0,
            min: f64::MAX,
            max: f64::MIN,
            nan_count: 0,
            inf_count: 0,
            zero_count: 0,
            total_count: 0,
            zero_ratio: 0.0,
            range: 0.0,
        }
    }

    fn update(&mut self, value: f64) {
        self.total_count += 1;

        if value.is_nan() {
            self.nan_count += 1;
            return;
        }

        if value.is_infinite() {
            self.inf_count += 1;
            return;
        }

        if value.abs() < 1e-10 {
            self.zero_count += 1;
        }

        self.min = self.min.min(value);
        self.max = self.max.max(value);
    }

    fn finalize(&mut self, values: &[f64]) {
        self.zero_ratio = self.zero_count as f64 / self.total_count as f64;
        self.range = self.max - self.min;

        // Compute mean and std dev (excluding NaN/Inf)
        let valid_values: Vec<f64> = values.iter().filter(|v| v.is_finite()).copied().collect();

        if !valid_values.is_empty() {
            self.mean = valid_values.iter().sum::<f64>() / valid_values.len() as f64;

            let variance = valid_values
                .iter()
                .map(|v| (v - self.mean).powi(2))
                .sum::<f64>()
                / valid_values.len() as f64;

            self.std_dev = variance.sqrt();
        }
    }

    /// Check if feature is problematic
    fn is_problematic(&self) -> Vec<String> {
        let mut issues = Vec::new();

        // Check 1: NaN/Inf values
        if self.nan_count > 0 {
            issues.push(format!(
                "NaN values: {}/{} ({:.2}%)",
                self.nan_count,
                self.total_count,
                100.0 * self.nan_count as f64 / self.total_count as f64
            ));
        }

        if self.inf_count > 0 {
            issues.push(format!(
                "Inf values: {}/{} ({:.2}%)",
                self.inf_count,
                self.total_count,
                100.0 * self.inf_count as f64 / self.total_count as f64
            ));
        }

        // Check 2: Constant features (std dev < 1e-6)
        if self.std_dev < 1e-6 {
            issues.push(format!("Constant feature (std={:.2e})", self.std_dev));
        }

        // Check 3: Sparse features (>95% zeros)
        if self.zero_ratio > 0.95 {
            issues.push(format!(
                "Sparse feature ({:.1}% zeros)",
                self.zero_ratio * 100.0
            ));
        }

        // Check 4: Extreme values (>100σ from mean)
        if self.std_dev > 0.0 {
            let max_z_score = ((self.max - self.mean) / self.std_dev).abs();
            let min_z_score = ((self.min - self.mean) / self.std_dev).abs();
            let extreme_z = max_z_score.max(min_z_score);

            if extreme_z > 100.0 {
                issues.push(format!("Extreme outlier (z-score={:.1f})", extreme_z));
            }
        }

        // Check 5: Unreasonably large range
        if self.range > 10000.0 {
            issues.push(format!(
                "Large range: [{:.2}, {:.2}] (range={:.2})",
                self.min, self.max, self.range
            ));
        }

        issues
    }
}

/// Compute correlation coefficient between two feature vectors
fn compute_correlation(vec1: &[f64], vec2: &[f64]) -> f64 {
    assert_eq!(vec1.len(), vec2.len());

    let n = vec1.len() as f64;
    let mean1 = vec1.iter().sum::<f64>() / n;
    let mean2 = vec2.iter().sum::<f64>() / n;

    let mut cov = 0.0;
    let mut var1 = 0.0;
    let mut var2 = 0.0;

    for i in 0..vec1.len() {
        let d1 = vec1[i] - mean1;
        let d2 = vec2[i] - mean2;
        cov += d1 * d2;
        var1 += d1 * d1;
        var2 += d2 * d2;
    }

    if var1 > 0.0 && var2 > 0.0 {
        cov / (var1 * var2).sqrt()
    } else {
        0.0
    }
}

#[test]
fn test_feature_quality_comprehensive() -> Result<()> {
    println!("\n=== Feature Quality Validation Test ===\n");

    // Load real market data
    let loader = RealDataLoader::new();
    let bars = loader.load_ohlcv_bars_from_parquet("test_data/ES_FUT_180d.parquet")?;

    println!("Loaded {} bars from ES_FUT_180d.parquet", bars.len());

    // Extract all features
    let feature_vectors = extract_ml_features(&bars)?;
    println!(
        "Extracted {} feature vectors (225 dimensions each)\n",
        feature_vectors.len()
    );

    // Initialize feature stats
    let mut stats: Vec<FeatureStats> = (0..225).map(FeatureStats::new).collect();

    // Collect all values for each feature
    let mut all_values: Vec<Vec<f64>> = vec![Vec::new(); 225];

    for feature_vec in &feature_vectors {
        for (i, &value) in feature_vec.iter().enumerate() {
            stats[i].update(value);
            all_values[i].push(value);
        }
    }

    // Finalize statistics
    for (i, stat) in stats.iter_mut().enumerate() {
        stat.finalize(&all_values[i]);
    }

    // === VALIDATION 1: Check for problematic features ===
    println!("=== VALIDATION 1: Problematic Features ===\n");

    let mut problematic_count = 0;
    let feature_names = get_feature_names();

    for stat in &stats {
        let issues = stat.is_problematic();
        if !issues.is_empty() {
            problematic_count += 1;
            println!(
                "⚠️  Feature {} ({}): {}",
                stat.index,
                feature_names.get(&stat.index).unwrap_or(&"Unknown"),
                issues.join(", ")
            );
        }
    }

    if problematic_count == 0 {
        println!("✅ No problematic features detected");
    } else {
        println!("\n❌ Found {} problematic features", problematic_count);
    }

    // === VALIDATION 2: Multicollinearity Analysis ===
    println!("\n=== VALIDATION 2: Multicollinearity Analysis ===\n");

    let mut high_corr_pairs = Vec::new();

    for i in 0..225 {
        for j in (i + 1)..225 {
            let corr = compute_correlation(&all_values[i], &all_values[j]);
            if corr.abs() > 0.95 {
                high_corr_pairs.push((i, j, corr));
            }
        }
    }

    if high_corr_pairs.is_empty() {
        println!("✅ No highly correlated feature pairs (>0.95)");
    } else {
        println!(
            "❌ Found {} highly correlated feature pairs (>0.95):\n",
            high_corr_pairs.len()
        );

        for (i, j, corr) in high_corr_pairs.iter().take(20) {
            println!(
                "  Feature {} ({}) <-> Feature {} ({}): corr={:.3}",
                i,
                feature_names.get(i).unwrap_or(&"Unknown"),
                j,
                feature_names.get(j).unwrap_or(&"Unknown"),
                corr
            );
        }

        if high_corr_pairs.len() > 20 {
            println!("  ... and {} more pairs", high_corr_pairs.len() - 20);
        }
    }

    // === VALIDATION 3: Feature Group Analysis ===
    println!("\n=== VALIDATION 3: Feature Group Analysis ===\n");

    let groups = vec![
        ("OHLCV", 0, 5),
        ("Technical Indicators", 5, 15),
        ("Price Patterns", 15, 75),
        ("Volume Patterns", 75, 115),
        ("Microstructure Proxies", 115, 165),
        ("Time Features", 165, 175),
        ("Statistical Features", 175, 201),
        ("Regime Detection", 201, 225),
    ];

    for (group_name, start, end) in groups {
        let group_stats = &stats[start..end];

        let nan_count: usize = group_stats.iter().map(|s| s.nan_count).sum();
        let inf_count: usize = group_stats.iter().map(|s| s.inf_count).sum();
        let constant_count = group_stats.iter().filter(|s| s.std_dev < 1e-6).count();
        let sparse_count = group_stats.iter().filter(|s| s.zero_ratio > 0.95).count();

        let avg_std = group_stats.iter().map(|s| s.std_dev).sum::<f64>() / group_stats.len() as f64;
        let max_range = group_stats.iter().map(|s| s.range).fold(0.0, f64::max);

        println!("{} (features {}-{}):", group_name, start, end - 1);
        println!("  Size: {} features", end - start);
        println!("  NaN count: {}", nan_count);
        println!("  Inf count: {}", inf_count);
        println!("  Constant features: {}", constant_count);
        println!("  Sparse features (>95% zero): {}", sparse_count);
        println!("  Average std dev: {:.4}", avg_std);
        println!("  Max range: {:.2}", max_range);

        if nan_count > 0 || inf_count > 0 || constant_count > 0 || max_range > 1000.0 {
            println!("  ⚠️  THIS GROUP HAS QUALITY ISSUES");
        }
        println!();
    }

    // === VALIDATION 4: Distribution Analysis ===
    println!("=== VALIDATION 4: Distribution Analysis (Top 10 Most Volatile) ===\n");

    let mut sorted_by_std: Vec<_> = stats.iter().collect();
    sorted_by_std.sort_by(|a, b| b.std_dev.partial_cmp(&a.std_dev).unwrap());

    for stat in sorted_by_std.iter().take(10) {
        println!(
            "Feature {} ({}): std={:.4}, range=[{:.2}, {:.2}]",
            stat.index,
            feature_names.get(&stat.index).unwrap_or(&"Unknown"),
            stat.std_dev,
            stat.min,
            stat.max
        );
    }

    // === FINAL VERDICT ===
    println!("\n=== FINAL VERDICT ===\n");

    let total_issues = problematic_count + high_corr_pairs.len();

    if total_issues == 0 {
        println!("✅ All 225 features passed quality checks");
        println!("Features are NOT the cause of gradient explosions.");
        println!("Investigate: learning rate, reward function, network architecture.");
    } else {
        println!("❌ Feature quality issues detected:");
        println!("  - {} problematic features", problematic_count);
        println!("  - {} highly correlated pairs", high_corr_pairs.len());
        println!("\n📊 RECOMMENDATIONS:");
        println!(
            "  1. Remove Statistical Features (175-200): Likely unstable (skewness, kurtosis)"
        );
        println!("  2. Review Microstructure Proxies (115-164): Check for division by zero issues");
        println!("  3. Prune highly correlated features: Keep only 1 from each correlated pair");
        println!("  4. Target feature count: 20-60 (currently 225, likely excessive)");
        println!("\n⚠️  High probability these features are causing gradient explosions!");
    }

    Ok(())
}

/// Get human-readable feature names for indices 0-224
fn get_feature_names() -> HashMap<usize, &'static str> {
    let mut names = HashMap::new();

    // OHLCV (0-4)
    names.insert(0, "Open");
    names.insert(1, "High");
    names.insert(2, "Low");
    names.insert(3, "Close");
    names.insert(4, "Volume");

    // Technical Indicators (5-14)
    names.insert(5, "RSI");
    names.insert(6, "EMA_Fast");
    names.insert(7, "EMA_Slow");
    names.insert(8, "MACD_Line");
    names.insert(9, "MACD_Signal");
    names.insert(10, "MACD_Histogram");
    names.insert(11, "Bollinger_Middle");
    names.insert(12, "Bollinger_Upper");
    names.insert(13, "Bollinger_Lower");
    names.insert(14, "ATR");

    // Price Patterns (15-74) - just mark ranges
    for i in 15..75 {
        names.insert(i, "Price_Pattern");
    }

    // Volume Patterns (75-114)
    for i in 75..115 {
        names.insert(i, "Volume_Pattern");
    }

    // Microstructure Proxies (115-164)
    names.insert(115, "Roll_Spread");
    names.insert(116, "Amihud_Illiquidity");
    names.insert(117, "Corwin_Schultz_Spread");
    for i in 118..165 {
        names.insert(i, "Microstructure");
    }

    // Time Features (165-174)
    for i in 165..175 {
        names.insert(i, "Time_Feature");
    }

    // Statistical Features (175-200)
    for i in 175..201 {
        names.insert(i, "Statistical");
    }

    // Regime Detection (201-224)
    for i in 201..225 {
        names.insert(i, "Regime");
    }

    names
}

#[test]
fn test_feature_extraction_basic_sanity() -> Result<()> {
    println!("\n=== Basic Feature Extraction Sanity Test ===\n");

    // Load data
    let loader = RealDataLoader::new();
    let bars = loader.load_ohlcv_bars_from_parquet("test_data/ES_FUT_180d.parquet")?;

    // Extract features
    let feature_vectors = extract_ml_features(&bars)?;

    println!("✅ Extracted {} feature vectors", feature_vectors.len());

    // Check dimensions
    for (i, feature_vec) in feature_vectors.iter().enumerate() {
        assert_eq!(
            feature_vec.len(),
            225,
            "Feature vector {} has wrong dimension: expected 225, got {}",
            i,
            feature_vec.len()
        );
    }

    println!("✅ All feature vectors have correct dimension (225)");

    // Check for NaN/Inf in first 10 vectors
    let mut has_nan_inf = false;
    for (i, feature_vec) in feature_vectors.iter().take(10).enumerate() {
        for (j, &value) in feature_vec.iter().enumerate() {
            if !value.is_finite() {
                println!(
                    "❌ Feature vector {} has invalid value at index {}: {}",
                    i, j, value
                );
                has_nan_inf = true;
            }
        }
    }

    if !has_nan_inf {
        println!("✅ No NaN/Inf values in first 10 feature vectors");
    }

    Ok(())
}