linreg-core 0.8.1

// ============================================================================
// Prediction Intervals Validation
// ============================================================================
//
// Validates prediction interval computation against R and Python reference
// values generated by:
//   R:      Rscript verification/scripts/r/core/test_prediction_intervals.R
//   Python: python verification/scripts/python/core/test_prediction_intervals.py
//
// Compares: predicted values, lower/upper bounds, SE_pred, leverage, df_residuals

use crate::common::{load_dataset_with_encoding, load_pi_result, CategoricalEncoding, PI_TOLERANCE};

use linreg_core::prediction_intervals::prediction_intervals;

/// Test datasets for prediction interval validation.
/// Matches OLS/WLS dataset lists with documented exclusions.
const PI_DATASETS: &[&str] = &[
    "bodyfat",
    "cars_stopping",
    "faithful",
    "iris",
    "lh",
    "longley",
    "mtcars",
    "prostate",
    // "synthetic_collinear",   // Excluded: perfect collinearity causes singular (X'X) matrix
    // "synthetic_high_vif",    // Excluded: high VIF (>5) causes numerical instability in leverage
    "synthetic_simple_linear",
    "synthetic_multiple",
    "synthetic_heteroscedastic",
    "synthetic_nonlinear",
    "synthetic_nonnormal",
    "synthetic_autocorrelated",
    "synthetic_outliers",
    "synthetic_small",
    "synthetic_interaction",
    "ToothGrowth",
];

/// Compute max absolute error and max percentage difference across two slices.
/// Returns (max_abs_err, max_pct_diff, worst_index).
fn compare_vectors(rust: &[f64], reference: &[f64]) -> (f64, f64, usize) {
    let mut max_abs = 0.0f64;
    let mut max_pct = 0.0f64;
    let mut worst_idx = 0;
    for (i, (r, ref_val)) in rust.iter().zip(reference.iter()).enumerate() {
        let abs_err = (r - ref_val).abs();
        if abs_err > max_abs {
            max_abs = abs_err;
            worst_idx = i;
        }
        if ref_val.abs() > 1e-15 {
            max_pct = max_pct.max(abs_err / ref_val.abs() * 100.0);
        }
    }
    (max_abs, max_pct, worst_idx)
}

/// Validate prediction intervals against a single reference (R or Python).
/// Uses the correct categorical encoding for each language:
///   R:      1-based (as.numeric(as.factor(...)))
///   Python: 0-based (pd.factorize(...))
fn validate_dataset_against_ref(
    dataset_name: &str,
    lang: &str,
    encoding: CategoricalEncoding,
    ref_dir: &std::path::Path,
    datasets_dir: &std::path::Path,
) -> Option<bool> {
    let csv_path = datasets_dir.join(format!("{}.csv", dataset_name));
    let ref_json = ref_dir.join(format!("{}_prediction_intervals.json", dataset_name));

    if !csv_path.exists() {
        return None;
    }

    let reference = match load_pi_result(&ref_json) {
        Some(r) => r,
        None => {
            println!("      {} reference not found", lang);
            return None;
        }
    };

    let dataset = match load_dataset_with_encoding(&csv_path, encoding) {
        Ok(d) => d,
        Err(e) => {
            println!("      Failed to load dataset: {}", e);
            return Some(false);
        }
    };

    let x_refs: Vec<&[f64]> = dataset.x_vars.iter().map(|v| v.as_slice()).collect();
    let result = match prediction_intervals(&dataset.y, &dataset.x_vars, &x_refs, reference.alpha) {
        Ok(r) => r,
        Err(e) => {
            println!("      Rust PI computation failed: {}", e);
            return Some(false);
        }
    };

    let (pred_abs, pred_pct, pred_idx) = compare_vectors(&result.predicted, &reference.train.predicted);
    let (lower_abs, lower_pct, lower_idx) = compare_vectors(&result.lower_bound, &reference.train.lower);
    let (upper_abs, upper_pct, upper_idx) = compare_vectors(&result.upper_bound, &reference.train.upper);
    let (se_abs, se_pct, se_idx) = compare_vectors(&result.se_pred, &reference.train.se_pred);

    let all_ok = pred_abs < PI_TOLERANCE
        && lower_abs < PI_TOLERANCE
        && upper_abs < PI_TOLERANCE
        && se_abs < PI_TOLERANCE;

    println!("    {} comparison (n={}):", lang, reference.train.predicted.len());
    println!("      predicted:   diff = {:.2e}, pct = {:.2e}% (worst idx={})", pred_abs, pred_pct, pred_idx);
    if pred_abs > PI_TOLERANCE {
        println!("                   Rust = {:.10}, {} = {:.10}", result.predicted[pred_idx], lang, reference.train.predicted[pred_idx]);
    }
    println!("      lower_bound: diff = {:.2e}, pct = {:.2e}% (worst idx={})", lower_abs, lower_pct, lower_idx);
    if lower_abs > PI_TOLERANCE {
        println!("                   Rust = {:.10}, {} = {:.10}", result.lower_bound[lower_idx], lang, reference.train.lower[lower_idx]);
    }
    println!("      upper_bound: diff = {:.2e}, pct = {:.2e}% (worst idx={})", upper_abs, upper_pct, upper_idx);
    if upper_abs > PI_TOLERANCE {
        println!("                   Rust = {:.10}, {} = {:.10}", result.upper_bound[upper_idx], lang, reference.train.upper[upper_idx]);
    }
    println!("      se_pred:     diff = {:.2e}, pct = {:.2e}% (worst idx={})", se_abs, se_pct, se_idx);
    if se_abs > PI_TOLERANCE {
        println!("                   Rust = {:.10}, {} = {:.10}", result.se_pred[se_idx], lang, reference.train.se_pred[se_idx]);
    }

    let status = if all_ok { "PASS" } else { "FAIL" };
    println!("      {} validation: {}", lang, status);

    Some(all_ok)
}

#[test]
fn validate_prediction_intervals_all_datasets() {
    println!("\n");
    println!("╔══════════════════════════════════════════════════════════════════════╗");
    println!("║  PREDICTION INTERVALS - COMPREHENSIVE MULTI-DATASET VALIDATION     ║");
    println!("╚══════════════════════════════════════════════════════════════════════╝");
    println!();

    let current_dir = std::env::current_dir().expect("Failed to get current dir");
    let datasets_dir = current_dir.join("verification/datasets/csv");
    let r_results_dir = current_dir.join("verification/results/r");
    let py_results_dir = current_dir.join("verification/results/python");

    let mut total_tests = 0;
    let mut passed_r = 0;
    let mut passed_python = 0;
    let mut failed_tests: Vec<(String, String)> = Vec::new();

    for dataset_name in PI_DATASETS {
        println!("  ┌─────────────────────────────────────────────────────────────────┐");
        println!("  │  Dataset: {:<52}│", dataset_name);
        println!("  └─────────────────────────────────────────────────────────────────┘");

        // Validate against R (1-based categorical encoding)
        if let Some(ok) = validate_dataset_against_ref(
            dataset_name, "R", CategoricalEncoding::OneBased,
            &r_results_dir, &datasets_dir,
        ) {
            total_tests += 1;
            if ok {
                passed_r += 1;
            } else {
                failed_tests.push((dataset_name.to_string(), "R mismatch".to_string()));
            }
        }

        // Validate against Python (0-based categorical encoding)
        if let Some(ok) = validate_dataset_against_ref(
            dataset_name, "Python", CategoricalEncoding::ZeroBased,
            &py_results_dir, &datasets_dir,
        ) {
            total_tests += 1;
            if ok {
                passed_python += 1;
            } else {
                failed_tests.push((dataset_name.to_string(), "Python mismatch".to_string()));
            }
        }

        println!();
    }

    // Summary
    println!("╔══════════════════════════════════════════════════════════════════════╗");
    println!("║  VALIDATION SUMMARY                                                ║");
    println!("╠══════════════════════════════════════════════════════════════════════╣");
    println!("║  Total tests run:            {:>35}  ║", total_tests);
    println!("║  R validations passed:       {:>35}  ║", passed_r);
    println!("║  Python validations passed:  {:>35}  ║", passed_python);
    println!("║  Failed tests:               {:>35}  ║", failed_tests.len());
    println!("║  Tolerance:                  {:>35}  ║", format!("{:.0e}", PI_TOLERANCE));
    println!("╚══════════════════════════════════════════════════════════════════════╝");

    if !failed_tests.is_empty() {
        println!();
        println!("Failed tests:");
        for (dataset, reason) in &failed_tests {
            println!("  - {}: {}", dataset, reason);
        }
    }

    assert!(total_tests > 0, "No validation tests were run. Check that reference files exist.");

    let pass_rate = (passed_r + passed_python) as f64 / total_tests as f64;
    assert!(
        pass_rate >= 0.9,
        "Validation pass rate ({:.1}%) is below 90% threshold.",
        pass_rate * 100.0
    );

    println!();
    println!(" Prediction intervals comprehensive validation passed!");
}