//! Vector-tomography parity against tomopy (`libtomo/recon/vector.c`).
//!
//! Loads golden data generated by `tools/gen_tomopy_vector_golden.py`
//! (tomopy 1.15) and reconstructs the SAME inputs with tomoxide's port. Because
//! the Rust code is a line-for-line port of tomopy's C kernels (same mixed
//! float/double arithmetic, same `calc_quadrant` int trick, same
//! truncate-then-correct index rule), the reconstructions must match to the
//! f32 rounding floor.
//!
//! vector2/vector3 reconstruct a full 3-D vector field (`dx**3`), so their
//! golden uses a cube (`dy == dx`); the single-dataset `vector` uses a
//! rectangular stack (`dy != dx`) to exercise the general slice loop.

use ndarray::{Array1, Array3};
use ndarray_npy::read_npy;
use tomoxide::recon::vector::{vector, vector2, vector3};

const FIXTURES: &str = concat!(env!("CARGO_MANIFEST_DIR"), "/tests/fixtures");

fn load3(name: &str) -> Array3<f32> {
    read_npy(format!("{FIXTURES}/{name}")).unwrap_or_else(|e| panic!("load {name}: {e}"))
}

fn load1(name: &str) -> Array1<f32> {
    read_npy(format!("{FIXTURES}/{name}")).unwrap_or_else(|e| panic!("load {name}: {e}"))
}

/// Max absolute difference between two equally shaped arrays.
fn max_abs_diff(a: &Array3<f32>, b: &Array3<f32>) -> f32 {
    assert_eq!(a.dim(), b.dim(), "shape mismatch");
    a.iter()
        .zip(b.iter())
        .map(|(x, y)| (x - y).abs())
        .fold(0.0f32, f32::max)
}

// The port reproduces tomopy bit-for-bit (Δ = 0 on this machine — the same
// libm sin/cos/sqrt back both sides). A hair of slack guards against a stray
// ULP from a different system libm, without weakening the gate.
const TOL: f32 = 1e-6;

#[test]
fn vector_matches_tomopy() {
    let tomo = load3("vector_tomo1.npy"); // (dt, dy, dx) = (16, 2, 16)
    let theta = load1("vector_theta1.npy");
    let (r1, r2) = vector(tomo.view(), theta.as_slice().unwrap(), None, 4).unwrap();

    let g1 = load3("vector_out1.npy");
    let g2 = load3("vector_out2.npy");
    let d1 = max_abs_diff(&r1, &g1);
    let d2 = max_abs_diff(&r2, &g2);
    assert!(d1 < TOL && d2 < TOL, "vector mismatch: d1={d1}, d2={d2}");
}

#[test]
fn vector2_matches_tomopy() {
    let t1 = load3("vector2_tomo1.npy"); // cube (8, 8, 8)
    let t2 = load3("vector2_tomo2.npy");
    let theta = load1("vector2_theta1.npy");
    let th = theta.as_slice().unwrap();
    let (r1, r2, r3) = vector2(t1.view(), t2.view(), th, th, None, None, 4, 1, 2).unwrap();

    let d1 = max_abs_diff(&r1, &load3("vector2_out1.npy"));
    let d2 = max_abs_diff(&r2, &load3("vector2_out2.npy"));
    let d3 = max_abs_diff(&r3, &load3("vector2_out3.npy"));
    assert!(
        d1 < TOL && d2 < TOL && d3 < TOL,
        "vector2 mismatch: d1={d1}, d2={d2}, d3={d3}"
    );
}

#[test]
fn vector3_matches_tomopy() {
    let t1 = load3("vector3_tomo1.npy"); // cube (8, 8, 8)
    let t2 = load3("vector3_tomo2.npy");
    let t3 = load3("vector3_tomo3.npy");
    let theta = load1("vector3_theta1.npy");
    let th = theta.as_slice().unwrap();
    let (r1, r2, r3) = vector3(
        t1.view(),
        t2.view(),
        t3.view(),
        th,
        th,
        th,
        None,
        None,
        None,
        4,
        0,
        1,
        2,
    )
    .unwrap();

    let d1 = max_abs_diff(&r1, &load3("vector3_out1.npy"));
    let d2 = max_abs_diff(&r2, &load3("vector3_out2.npy"));
    let d3 = max_abs_diff(&r3, &load3("vector3_out3.npy"));
    assert!(
        d1 < TOL && d2 < TOL && d3 < TOL,
        "vector3 mismatch: d1={d1}, d2={d2}, d3={d3}"
    );
}

#[test]
fn vector2_requires_cube() {
    // dy != dx must be rejected (tomopy corrupts memory here).
    let t = Array3::<f32>::zeros((4, 2, 8));
    let theta = vec![0.0f32; 4];
    let err = vector2(t.view(), t.view(), &theta, &theta, None, None, 1, 1, 2);
    assert!(err.is_err(), "vector2 should reject non-cube dy != dx");
}
