feffit 0.1.0

//! Parity tests for the `feffdat` port.
//!
//! Two independent checks:
//!   1. The parser is verified against values read by hand from `feff0001.dat`
//!      (non-circular: literal numbers transcribed from the file).
//!   2. The EXAFS equation + linear interpolation is verified against the
//!      numpy reference in `tests/data/ref_*.txt` (generated by
//!      `scripts/ref_chi.py`, which replicates larch's `interp='lin'` branch).
//!
//! The cubic-spline path (larch's default) only gets a smoke test here: scipy
//! is unavailable, so its numerical parity is unverified.

use std::collections::HashMap;
use std::path::PathBuf;

use feffit::feffdat::{FeffDatFile, FeffPath, Interp, KGrid, PathParams, ff2chi};

fn data_dir() -> PathBuf {
    PathBuf::from(env!("CARGO_MANIFEST_DIR")).join("tests/data")
}

struct Reference {
    params: HashMap<String, f64>,
    scalars: HashMap<String, String>,
    k: Vec<f64>,
    chi_lin: Vec<f64>,
    /// Present only when the reference was generated with scipy (cubic column).
    chi_cubic: Option<Vec<f64>>,
}

fn load_ref(name: &str) -> Reference {
    let text = std::fs::read_to_string(data_dir().join(name)).unwrap();
    let mut params = HashMap::new();
    let mut scalars = HashMap::new();
    let mut k = Vec::new();
    let mut chi_lin = Vec::new();
    let mut chi_cubic = Vec::new();
    let mut in_data = false;
    for line in text.lines() {
        if line == "#data" {
            in_data = true;
            continue;
        }
        if let Some(rest) = line.strip_prefix("#param ") {
            let mut it = rest.split_whitespace();
            let key = it.next().unwrap().to_string();
            let val: f64 = it.next().unwrap().parse().unwrap();
            params.insert(key, val);
        } else if let Some(rest) = line.strip_prefix("#scalar ") {
            let mut it = rest.split_whitespace();
            let key = it.next().unwrap().to_string();
            scalars.insert(key, it.next().unwrap().to_string());
        } else if in_data && !line.starts_with('#') && !line.is_empty() {
            let mut it = line.split_whitespace();
            k.push(it.next().unwrap().parse().unwrap());
            chi_lin.push(it.next().unwrap().parse().unwrap());
            if let Some(tok) = it.next() {
                chi_cubic.push(tok.parse().unwrap());
            }
        }
    }
    let chi_cubic = if chi_cubic.len() == k.len() {
        Some(chi_cubic)
    } else {
        None
    };
    Reference {
        params,
        scalars,
        k,
        chi_lin,
        chi_cubic,
    }
}

fn params_from_ref(r: &Reference) -> PathParams {
    let g = |k: &str| *r.params.get(k).unwrap();
    PathParams {
        degen: g("degen"),
        s02: g("s02"),
        e0: g("e0"),
        ei: g("ei"),
        deltar: g("deltar"),
        sigma2: g("sigma2"),
        third: g("third"),
        fourth: g("fourth"),
    }
}

fn max_abs_diff(a: &[f64], b: &[f64]) -> f64 {
    assert_eq!(
        a.len(),
        b.len(),
        "length mismatch {} vs {}",
        a.len(),
        b.len()
    );
    a.iter()
        .zip(b)
        .map(|(x, y)| (x - y).abs())
        .fold(0.0, f64::max)
}

/// (1) Parser against values transcribed by hand from feff0001.dat.
#[test]
fn parser_matches_known_values() {
    let f = FeffDatFile::from_path(data_dir().join("feff0001.dat")).unwrap();

    assert_eq!(f.nleg, 2);
    assert_eq!(f.degen, 12.0);
    assert_eq!(f.reff, 2.5478);
    assert_eq!(f.rnorman, 2.6300);
    assert_eq!(f.edge, -5.53502);
    assert_eq!(f.shell, "K");
    assert_eq!(f.absorber.as_deref(), Some("Cu"));
    assert_eq!(f.gam_ch, 1.761e+00);

    // Mu line: Mu=-5.535 kf=1.806 Vint=-17.96 Rs_int=2.008
    assert_eq!(f.vmu, -5.535);
    assert_eq!(f.vint, -17.96);
    assert_eq!(f.rs_int, 2.008);

    // potentials: Abs (ipot 0) and Pot 1, both Z=29
    assert_eq!(f.potentials.len(), 2);
    assert_eq!(f.potentials[0].ipot, 0);
    assert_eq!(f.potentials[0].iz, 29);
    assert_eq!(f.potentials[1].ipot, 1);

    // geometry: absorber Cu at origin, scatterer Cu
    assert_eq!(f.geom.len(), 2);
    assert_eq!(f.geom[0].label, "Cu");
    assert_eq!(f.geom[0].ipot, 0);
    assert_eq!(f.geom[1].ipot, 1);
    assert_eq!(f.geom[1].y, -1.8016);
    assert_eq!(f.geom[1].z, 1.8016);

    // array block: 59 rows, k 0..20
    assert_eq!(f.k.len(), 59);
    assert_eq!(f.k[0], 0.0);
    assert_eq!(f.k[1], 0.1);
    assert_eq!(*f.k.last().unwrap(), 20.0);

    // derived columns from the first data row:
    // real_phc=3.5015, mag_feff=0.0, pha_feff=-5.1320, red_fact=0.9872
    // pha = real_phc + pha_feff ; amp = mag_feff * red_fact
    assert!((f.pha[0] - (3.5015 - 5.1320)).abs() < 1e-12);
    assert_eq!(f.amp[0], 0.0); // mag_feff[0] == 0
    assert!((f.real_phc[0] - 3.5015).abs() < 1e-12);
    assert!((f.lam[0] - 1.5637e+01).abs() < 1e-9);
    assert!((f.rep[0] - 1.8071e+00).abs() < 1e-9);
}

/// Cross-check the scalars the reference recorded against the parser.
fn assert_scalars_match(f: &FeffDatFile, r: &Reference) {
    let s = |k: &str| r.scalars.get(k).unwrap().clone();
    assert_eq!(f.reff, s("reff").parse::<f64>().unwrap());
    assert_eq!(f.nleg, s("nleg").parse::<usize>().unwrap());
    assert_eq!(f.degen, s("file_degen").parse::<f64>().unwrap());
    assert_eq!(f.k.len(), s("n_kgrid").parse::<usize>().unwrap());
    assert_eq!(f.k[0], s("dat_k_first").parse::<f64>().unwrap());
    assert_eq!(
        *f.k.last().unwrap(),
        s("dat_k_last").parse::<f64>().unwrap()
    );
    assert!((f.pha[0] - s("dat_pha_first").parse::<f64>().unwrap()).abs() < 1e-12);
    assert!((f.amp.last().unwrap() - s("dat_amp_last").parse::<f64>().unwrap()).abs() < 1e-12);
}

fn check_linear_against_ref(dat: &str, ref_name: &str, tol: f64) {
    let r = load_ref(ref_name);
    let f = FeffDatFile::from_path(data_dir().join(dat)).unwrap();
    assert_scalars_match(&f, &r);

    let mut path = FeffPath::new(f).with_params(params_from_ref(&r));
    path.calc_chi(&KGrid::default_step(), Interp::Linear);

    // k grid must match the reference exactly
    assert_eq!(path.k.len(), r.k.len());
    for (a, b) in path.k.iter().zip(&r.k) {
        assert_eq!(a, b, "k grid differs");
    }
    let diff = max_abs_diff(&path.chi, &r.chi_lin);
    println!("{ref_name}: max|chi - chi_lin| = {diff:.3e}");
    assert!(
        diff < tol,
        "{ref_name}: linear chi diff {diff:.3e} exceeds {tol:.0e}"
    );
}

/// (2) EXAFS equation + linear interp, default parameters.
#[test]
fn chi_linear_default_params() {
    check_linear_against_ref("feff0001.dat", "ref_cu_default.txt", 1e-9);
}

/// (2) EXAFS equation + linear interp, with e0 shift + sigma2 + deltar + third.
#[test]
fn chi_linear_shifted_params() {
    check_linear_against_ref("feff0001.dat", "ref_cu_shift.txt", 1e-9);
}

/// ff2chi: sum of three paths (default params) vs sum of the three references.
#[test]
fn ff2chi_sum_linear() {
    let refs = ["ref_cu_default.txt", "ref_cu_p2.txt", "ref_cu_p3.txt"];
    let dats = ["feff0001.dat", "feff0002.dat", "feff0003.dat"];

    let mut sum_ref: Vec<f64> = Vec::new();
    for name in refs {
        let r = load_ref(name);
        if sum_ref.is_empty() {
            sum_ref = vec![0.0; r.chi_lin.len()];
        }
        for (s, c) in sum_ref.iter_mut().zip(&r.chi_lin) {
            *s += c;
        }
    }

    let mut paths: Vec<FeffPath> = dats
        .iter()
        .map(|d| FeffPath::from_path(data_dir().join(d)).unwrap())
        .collect();
    let (_k, chi) = ff2chi(&mut paths, &KGrid::default_step(), Interp::Linear);

    let diff = max_abs_diff(&chi, &sum_ref);
    println!("ff2chi sum: max|chi - sum_ref| = {diff:.3e}");
    assert!(diff < 1e-9, "ff2chi sum diff {diff:.3e} exceeds 1e-9");
}

fn check_cubic_against_ref(dat: &str, ref_name: &str, tol: f64) {
    let r = load_ref(ref_name);
    let ref_cubic = r
        .chi_cubic
        .as_ref()
        .expect("reference has no chi_cubic column (regenerate with scipy installed)");
    let f = FeffDatFile::from_path(data_dir().join(dat)).unwrap();

    let mut path = FeffPath::new(f).with_params(params_from_ref(&r));
    path.calc_chi(&KGrid::default_step(), Interp::Cubic);

    let diff = max_abs_diff(&path.chi, ref_cubic);
    println!("{ref_name}: max|chi - chi_cubic(scipy)| = {diff:.3e}");
    assert!(
        diff < tol,
        "{ref_name}: cubic chi diff {diff:.3e} exceeds {tol:.0e}"
    );
}

/// EXAFS equation + not-a-knot cubic spline vs scipy `UnivariateSpline(s=0)`,
/// default parameters (q stays within the tabulated k range).
#[test]
fn chi_cubic_default_params() {
    check_cubic_against_ref("feff0001.dat", "ref_cu_default.txt", 1e-9);
}

/// Cubic parity with an e0 shift large enough to push low-k `q` negative, so
/// the spline is also exercised in its extrapolation region.
#[test]
fn chi_cubic_shifted_params() {
    check_cubic_against_ref("feff0001.dat", "ref_cu_shift.txt", 1e-9);
}