// SPDX-License-Identifier: LicenseRef-PolyForm-Noncommercial-1.0.0
//! Cross-solver validation tests.
//!
//! Loads reference solutions (generated by PYPOWER) and compares against our
//! Newton-Raphson solver. PYPOWER is a direct Python port of MATPOWER, so these
//! solutions are effectively MATPOWER-validated ground truth.

mod common;

use std::collections::HashMap;
use std::fs;

use serde::Deserialize;

#[derive(Deserialize)]
#[allow(dead_code)]
struct ReferenceSolution {
    solver: String,
    case: String,
    n_buses: usize,
    converged: bool,
    buses: Vec<BusSolution>,
}

#[derive(Deserialize)]
#[allow(dead_code)]
struct BusSolution {
    number: u32,
    vm_pu: f64,
    va_deg: f64,
    p_mw: f64,
    q_mvar: f64,
}

fn load_reference(solver: &str, case_name: &str) -> Option<ReferenceSolution> {
    let path =
        common::workspace_root().join(format!("benchmarks/references/{solver}/{case_name}.json"));
    if !path.exists() {
        return None;
    }
    let content = fs::read_to_string(&path)
        .unwrap_or_else(|e| panic!("failed to read {}: {e}", path.display()));
    Some(
        serde_json::from_str(&content)
            .unwrap_or_else(|e| panic!("failed to parse {}: {e}", path.display())),
    )
}

fn load_case(name: &str) -> surge_network::Network {
    let path = common::test_data_dir().join(format!("{name}.m"));
    surge_io::matpower::load(&path).unwrap_or_else(|e| panic!("failed to parse {name}: {e}"))
}

/// Validate our solver against a reference solution.
///
/// Tolerances are per the plan:
/// - Vm: 1e-4 p.u. for classic cases, looser for ACTIVS/Synthetic
/// - Va: 0.01 degrees for classic cases
/// - P/Q: 0.1 MW/Mvar (not checked for all cases due to different slack dispatch)
fn validate_case(case_name: &str, vm_tol: f64, va_tol_deg: f64) {
    let reference = match load_reference("pypower", case_name) {
        Some(r) => r,
        None => {
            eprintln!("SKIP {case_name}: no pypower reference found");
            return;
        }
    };

    assert!(
        reference.converged,
        "reference solution did not converge for {case_name}"
    );

    let network = load_case(case_name);
    // PYPOWER references were generated WITHOUT Q-limit enforcement
    // (PYPOWER default: ENFORCE_Q_LIMS = False).  Match that setting here.
    let options = surge_ac::AcPfOptions {
        enforce_q_limits: false,
        ..Default::default()
    };

    let solution = surge_ac::solve_ac_pf_kernel(&network, &options)
        .unwrap_or_else(|e| panic!("{case_name} NR-KLU failed: {e}"));

    assert_eq!(
        solution.status,
        surge_solution::SolveStatus::Converged,
        "{case_name}: solver did not converge"
    );

    // Build lookup: bus_number -> (vm, va_deg)
    let ref_buses: HashMap<u32, &BusSolution> =
        reference.buses.iter().map(|b| (b.number, b)).collect();

    assert_eq!(
        solution.bus_numbers.len(),
        reference.n_buses,
        "{case_name}: bus count mismatch (surge={}, ref={})",
        solution.bus_numbers.len(),
        reference.n_buses
    );

    let mut max_vm_diff = 0.0_f64;
    let mut max_va_diff = 0.0_f64;
    let mut worst_vm_bus = 0u32;
    let mut worst_va_bus = 0u32;

    for (i, &bus_num) in solution.bus_numbers.iter().enumerate() {
        let ref_bus = ref_buses
            .get(&bus_num)
            .unwrap_or_else(|| panic!("{case_name}: bus {bus_num} not found in reference"));

        let vm_diff = (solution.voltage_magnitude_pu[i] - ref_bus.vm_pu).abs();
        let va_deg = solution.voltage_angle_rad[i].to_degrees();
        let va_diff = (va_deg - ref_bus.va_deg).abs();

        if vm_diff > max_vm_diff {
            max_vm_diff = vm_diff;
            worst_vm_bus = bus_num;
        }
        if va_diff > max_va_diff {
            max_va_diff = va_diff;
            worst_va_bus = bus_num;
        }
    }

    eprintln!(
        "{case_name:25}: max|dVm|={max_vm_diff:.2e} (bus {worst_vm_bus}), \
         max|dVa|={max_va_diff:.4} deg (bus {worst_va_bus}), \
         {} iters, {:.3}ms",
        solution.iterations,
        solution.solve_time_secs * 1000.0
    );

    assert!(
        max_vm_diff < vm_tol,
        "{case_name}: Vm diff {max_vm_diff:.2e} > {vm_tol:.2e} at bus {worst_vm_bus}"
    );
    assert!(
        max_va_diff < va_tol_deg,
        "{case_name}: Va diff {max_va_diff:.4} deg > {va_tol_deg:.4} deg at bus {worst_va_bus}"
    );
}

// ── Classic MATPOWER cases ──
//
// Tolerances reflect the zero-impedance branch clamping fix (CRITICAL-02):
// branches with |r| < 1e-6 or |x| < 1e-6 use a 1e-6 substitution to prevent
// Y-bus singularity.  This introduces a small but intentional modeling
// difference vs PYPOWER's unclamped admittances, causing Vm differences up to
// ~1e-5 pu on large PEGASE cases.  Vm tolerance is set to 5e-5 pu (still 50x
// tighter than MATPOWER's 1e-3 pu convergence band) and Va to 0.05 deg.

#[test]
fn validate_case9() {
    validate_case("case9", 5e-5, 0.05);
}

#[test]
fn validate_case14() {
    validate_case("case14", 5e-5, 0.05);
}

#[test]
fn validate_case30() {
    validate_case("case30", 5e-5, 0.05);
}

#[test]
fn validate_case118() {
    validate_case("case118", 5e-5, 0.05);
}

#[test]
fn validate_case2383wp() {
    validate_case("case2383wp", 5e-5, 0.05);
}

#[test]
fn validate_case9241pegase() {
    validate_case("case9241pegase", 5e-5, 0.05);
}

#[test]
fn validate_case13659pegase() {
    validate_case("case13659pegase", 5e-5, 0.05);
}

// ── RTE case (moderate tolerances — slight modeling differences) ──

#[test]
fn validate_case6515rte() {
    // RTE case has small differences due to transformer modeling
    validate_case("case6515rte", 5e-4, 0.01);
}

// ── ACTIVS/Synthetic cases (loose tolerances — different reactive dispatch) ──

#[test]
fn validate_case_activsg2000() {
    validate_case("case_ACTIVSg2000", 0.05, 0.5);
}

#[test]
fn validate_case_activsg10k() {
    validate_case("case_ACTIVSg10k", 0.05, 0.5);
}

#[test]
fn validate_case_activsg25k() {
    validate_case("case_ACTIVSg25k", 0.05, 1.0);
}

#[test]
fn validate_case_activsg70k() {
    validate_case("case_ACTIVSg70k", 0.05, 2.0);
}

#[test]
fn validate_case_syntheticusa() {
    validate_case("case_SyntheticUSA", 0.1, 2.0);
}
