oxiproj-transformations 0.1.2

//! Batch 7-parameter Helmert transformation and least-squares parameter estimation.
//!
//! Provides:
//! - [`helmert7_batch_scalar`] — pure-Rust scalar implementation (always available)
//! - [`helmert7_batch_scalar_inv`] — inverse of the above
//! - `helmert7_batch_blas` — BLAS-accelerated variant (requires `blas` feature)
//! - [`estimate_helmert7_lsq`] — LSQ estimation of 7 Helmert parameters from
//!   paired control points via Cholesky normal equations (pure Rust, always available)

use oxiproj_core::ProjError;

/// The 7 parameters of a linearised Helmert transformation.
///
/// - `tx`, `ty`, `tz` — translations in metres
/// - `rx`, `ry`, `rz` — rotations in **radians** (coordinate-frame convention, linearised)
/// - `scale_ppm` — scale in parts-per-million
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Helmert7Params {
    /// Translation along X axis (metres).
    pub tx: f64,
    /// Translation along Y axis (metres).
    pub ty: f64,
    /// Translation along Z axis (metres).
    pub tz: f64,
    /// Rotation about X axis (radians, linearised).
    pub rx: f64,
    /// Rotation about Y axis (radians, linearised).
    pub ry: f64,
    /// Rotation about Z axis (radians, linearised).
    pub rz: f64,
    /// Scale factor in parts-per-million.
    pub scale_ppm: f64,
}

/// The result of a 7-parameter Helmert LSQ estimation: `(tx, ty, tz, rx, ry, rz, scale_ppm)`.
pub type Helmert7Est = (f64, f64, f64, f64, f64, f64, f64);

/// Apply the linearised 7-parameter Helmert transformation to a batch of 3-D points.
///
/// The linearised rotation matrix (coordinate-frame convention) is:
/// ```text
/// R ≈ [[1,  rz, -ry],
///      [-rz, 1,  rx],
///      [ ry,-rx,  1]]
/// ```
/// Transformed point: `X_out = (1 + scale_ppm * 1e-6) * R * X_in + T`
pub fn helmert7_batch_scalar(coords: &[[f64; 3]], p: &Helmert7Params) -> Vec<[f64; 3]> {
    let s = 1.0 + p.scale_ppm * 1e-6;
    coords
        .iter()
        .map(|&[x, y, z]| {
            let xo = s * (x + p.rz * y - p.ry * z) + p.tx;
            let yo = s * (-p.rz * x + y + p.rx * z) + p.ty;
            let zo = s * (p.ry * x - p.rx * y + z) + p.tz;
            [xo, yo, zo]
        })
        .collect()
}

/// Inverse of the linearised 7-parameter Helmert batch transform.
///
/// Uses the transpose of the linearised rotation matrix and inverts the scale
/// factor. Only accurate when the rotation angles and scale are small (the
/// linearised approximation).
pub fn helmert7_batch_scalar_inv(coords: &[[f64; 3]], p: &Helmert7Params) -> Vec<[f64; 3]> {
    let s = 1.0 + p.scale_ppm * 1e-6;
    coords
        .iter()
        .map(|&[x, y, z]| {
            let xi = (x - p.tx) / s;
            let yi = (y - p.ty) / s;
            let zi = (z - p.tz) / s;
            let xo = xi - p.rz * yi + p.ry * zi;
            let yo = p.rz * xi + yi - p.rx * zi;
            let zo = -p.ry * xi + p.rx * yi + zi;
            [xo, yo, zo]
        })
        .collect()
}

/// BLAS-accelerated 7-parameter Helmert batch transformation.
///
/// Builds the 3×3 rotation matrix R and applies it to all points via
/// `oxiblas::gemm`, treating the input as a 3×n matrix.
///
/// This is significantly faster than the scalar fallback for large batches
/// (n > ~100) due to SIMD-optimised GEMM kernels in oxiblas.
#[cfg(feature = "blas")]
pub fn helmert7_batch_blas(coords: &[[f64; 3]], p: &Helmert7Params) -> Vec<[f64; 3]> {
    use oxiblas::{gemm, Mat};

    let n = coords.len();
    if n == 0 {
        return Vec::new();
    }

    let s = 1.0 + p.scale_ppm * 1e-6;

    // Build the 3×3 rotation matrix R (coordinate-frame convention, linearised).
    // R = [[1, rz, -ry], [-rz, 1, rx], [ry, -rx, 1]]  (scaled by s)
    let mut r_mat: Mat<f64> = Mat::from_rows(&[
        &[1.0, p.rz, -p.ry],
        &[-p.rz, 1.0, p.rx],
        &[p.ry, -p.rx, 1.0],
    ]);
    // Scale all elements by s
    for col in 0..3 {
        for row in 0..3 {
            r_mat[(row, col)] *= s;
        }
    }

    // Build input matrix X: 3×n, column-major (each column = one point).
    // Mat::from_slice expects column-major data: element (i, j) at index i + j*3.
    let col_major_data: Vec<f64> = (0..n)
        .flat_map(|j| (0..3).map(move |i| coords[j][i]))
        .collect();
    let x_mat: Mat<f64> = Mat::from_slice(3, n, &col_major_data);

    // C = R * X  (3×3 · 3×n = 3×n)
    let mut c_mat: Mat<f64> = Mat::zeros(3, n);
    gemm(1.0, r_mat.as_ref(), x_mat.as_ref(), 0.0, c_mat.as_mut());

    // Add translation and collect
    (0..n)
        .map(|j| {
            [
                c_mat[(0, j)] + p.tx,
                c_mat[(1, j)] + p.ty,
                c_mat[(2, j)] + p.tz,
            ]
        })
        .collect()
}

/// GPU-accelerated 7-parameter Helmert batch transformation (feature `gpu`).
///
/// Applies the same linearised coordinate-frame transform as
/// [`helmert7_batch_scalar`] — `Xout = (1 + scale_ppm·1e-6)·R·Xin + T` — but with
/// one CUDA thread per 3-D point. The kernel is pure f64 arithmetic
/// (`fma.rn.f64` / `mul.f64` / `sub.f64`), so it reproduces the CPU result to
/// within a few ULP with no transcendental approximation.
///
/// This is the highest-throughput path for very large point batches (point
/// clouds, dense geodetic networks). When no CUDA device is reachable, or any
/// driver/launch error occurs, it transparently falls back to
/// [`helmert7_batch_scalar`], so callers never need to guard for GPU absence.
///
/// The transform runs on the GPU via the Pure-Rust `oxicuda` stack
/// (`libcuda.so` is loaded at runtime; no CUDA Toolkit is required to build).
#[cfg(feature = "gpu")]
pub fn helmert7_batch_gpu(coords: &[[f64; 3]], p: &Helmert7Params) -> Vec<[f64; 3]> {
    helmert7_batch_gpu_checked(coords, p).unwrap_or_else(|_| helmert7_batch_scalar(coords, p))
}

/// Fallible GPU path used by [`helmert7_batch_gpu`]. Returns `Err` — which
/// triggers the CPU fallback — when no CUDA device is present or a device-memory
/// / launch error occurs.
#[cfg(feature = "gpu")]
fn helmert7_batch_gpu_checked(
    coords: &[[f64; 3]],
    p: &Helmert7Params,
) -> Result<Vec<[f64; 3]>, oxiproj_core::gpu::CudaError> {
    use oxiproj_core::gpu::{self, ptx_math, DeviceBuffer};

    let n = coords.len();
    if n == 0 {
        return Ok(Vec::new());
    }
    let s = 1.0 + p.scale_ppm * 1e-6;
    // `[[f64;3]]` is already contiguous, but flatten explicitly so the host→device
    // copy works regardless of slice provenance.
    let flat: Vec<f64> = coords.iter().flatten().copied().collect();
    let mut out_host = vec![0.0f64; n * 3];

    gpu::with_kernel(
        "oxp_helmert7",
        "oxp_helmert7",
        |target| {
            let mut ptx = ptx_math::ptx_header(target);
            ptx.push_str(HELMERT7_PTX);
            ptx
        },
        |kernel, stream| {
            let d_in = gpu::upload(&flat, stream)?;
            let d_out = DeviceBuffer::<f64>::alloc(n * 3)?;
            // Argument order must match the kernel's `.param` declarations.
            let args = (
                d_in.as_device_ptr(),
                d_out.as_device_ptr(),
                n as u32,
                s,
                p.tx,
                p.ty,
                p.tz,
                p.rx,
                p.ry,
                p.rz,
            );
            gpu::launch_1d(kernel, stream, n as u32, &args)?;
            d_out.copy_to_host(&mut out_host)?;
            Ok(())
        },
    )?;

    Ok(out_host
        .chunks_exact(3)
        .map(|c| [c[0], c[1], c[2]])
        .collect())
}

/// PTX `.entry` for [`helmert7_batch_gpu`]: one thread per 3-D point computing
/// `Xout = s·R·Xin + T` with the linearised coordinate-frame rotation
/// `R = [[1,rz,-ry],[-rz,1,rx],[ry,-rx,1]]`. Concatenated after a target header.
#[cfg(feature = "gpu")]
const HELMERT7_PTX: &str = r"
.visible .entry oxp_helmert7(
    .param .u64 p_in,
    .param .u64 p_out,
    .param .u32 p_n,
    .param .f64 p_s,
    .param .f64 p_tx,
    .param .f64 p_ty,
    .param .f64 p_tz,
    .param .f64 p_rx,
    .param .f64 p_ry,
    .param .f64 p_rz
)
{
    .reg .pred  %p<1>;
    .reg .b32   %r<6>;
    .reg .b64   %rd<6>;
    .reg .f64   %fd<16>;

    mov.u32      %r0, %ntid.x;
    mov.u32      %r1, %ctaid.x;
    mov.u32      %r2, %tid.x;
    mad.lo.u32   %r3, %r1, %r0, %r2;       // idx = blockIdx.x*blockDim.x + threadIdx.x
    ld.param.u32 %r4, [p_n];
    setp.ge.u32  %p0, %r3, %r4;
    @%p0 bra     DONE;

    mul.wide.u32 %rd0, %r3, 24;            // byte offset = idx * 3 * 8
    ld.param.u64 %rd1, [p_in];
    add.u64      %rd2, %rd1, %rd0;
    ld.param.u64 %rd3, [p_out];
    add.u64      %rd4, %rd3, %rd0;

    ld.global.f64 %fd0, [%rd2];            // x
    ld.global.f64 %fd1, [%rd2+8];          // y
    ld.global.f64 %fd2, [%rd2+16];         // z

    ld.param.f64 %fd3, [p_s];
    ld.param.f64 %fd4, [p_tx];
    ld.param.f64 %fd5, [p_ty];
    ld.param.f64 %fd6, [p_tz];
    ld.param.f64 %fd7, [p_rx];
    ld.param.f64 %fd8, [p_ry];
    ld.param.f64 %fd9, [p_rz];

    // xo = s*(x + rz*y - ry*z) + tx
    fma.rn.f64   %fd10, %fd9, %fd1, %fd0;
    mul.f64      %fd11, %fd8, %fd2;
    sub.f64      %fd10, %fd10, %fd11;
    fma.rn.f64   %fd12, %fd3, %fd10, %fd4;

    // yo = s*(y + rx*z - rz*x) + ty
    fma.rn.f64   %fd10, %fd7, %fd2, %fd1;
    mul.f64      %fd11, %fd9, %fd0;
    sub.f64      %fd10, %fd10, %fd11;
    fma.rn.f64   %fd13, %fd3, %fd10, %fd5;

    // zo = s*(z + ry*x - rx*y) + tz
    fma.rn.f64   %fd10, %fd8, %fd0, %fd2;
    mul.f64      %fd11, %fd7, %fd1;
    sub.f64      %fd10, %fd10, %fd11;
    fma.rn.f64   %fd14, %fd3, %fd10, %fd6;

    st.global.f64 [%rd4],    %fd12;
    st.global.f64 [%rd4+8],  %fd13;
    st.global.f64 [%rd4+16], %fd14;

DONE:
    ret;
}
";

/// Estimate the 7-parameter Helmert transform from paired control points.
///
/// `source` and `target` are matched 3-D coordinate pairs. At least 3 pairs
/// (9 observations, 7 unknowns) are required; more pairs over-determine the
/// system and the result is a least-squares fit.
///
/// The design matrix follows the standard geodetic linearisation of the
/// 7-parameter transform:
/// ```text
/// δX   [1 0 0  0   -Z   Y   X·1e-6] [tx]
/// δY = [0 1 0  Z    0  -X   Y·1e-6] [ty]
/// δZ   [0 0 1 -Y    X   0   Z·1e-6] [tz]
///                                    [rx]
///                                    [ry]
///                                    [rz]
///                                    [scale_ppm]
/// ```
/// where δ = target − source, (X,Y,Z) are source coordinates.
///
/// Solved via normal equations `(AᵀA)·x = Aᵀb` with pure-Rust Cholesky
/// decomposition.
///
/// Returns `(tx, ty, tz, rx, ry, rz, scale_ppm)` as a [`Helmert7Est`].
pub fn estimate_helmert7_lsq(
    source: &[[f64; 3]],
    target: &[[f64; 3]],
) -> Result<Helmert7Est, ProjError> {
    if source.len() != target.len() {
        return Err(ProjError::InvalidCoord);
    }
    let n = source.len();
    if n < 3 {
        return Err(ProjError::InvalidCoord);
    }

    // Build A (3n × 7) and b (3n)
    let rows = 3 * n;
    const NCOLS: usize = 7;
    let mut a = vec![0.0f64; rows * NCOLS];
    let mut b_vec = vec![0.0f64; rows];

    for (i, (src, tgt)) in source.iter().zip(target.iter()).enumerate() {
        let [x, y, z] = *src;
        let base = i * 3;

        // Row base+0: δX equation  (from helmert7_batch_scalar: xo = s*(x + rz*y - ry*z) + tx)
        // δX ≈ tx + rz*y - ry*z + x*scale_ppm*1e-6
        a[base * NCOLS] = 1.0; // tx coefficient
                               // ty = 0, tz = 0, rx = 0
        a[base * NCOLS + 4] = -z; // ry coefficient
        a[base * NCOLS + 5] = y; // rz coefficient
        a[base * NCOLS + 6] = x * 1e-6; // scale_ppm coefficient

        // Row base+1: δY equation  (from helmert7_batch_scalar: yo = s*(-rz*x + y + rx*z) + ty)
        // δY ≈ ty + rx*z - rz*x + y*scale_ppm*1e-6
        a[(base + 1) * NCOLS + 1] = 1.0; // ty coefficient
        a[(base + 1) * NCOLS + 3] = z; // rx coefficient
                                       // ry = 0
        a[(base + 1) * NCOLS + 5] = -x; // rz coefficient
        a[(base + 1) * NCOLS + 6] = y * 1e-6; // scale_ppm coefficient

        // Row base+2: δZ equation  (from helmert7_batch_scalar: zo = s*(ry*x - rx*y + z) + tz)
        // δZ ≈ tz - rx*y + ry*x + z*scale_ppm*1e-6
        a[(base + 2) * NCOLS + 2] = 1.0; // tz coefficient
        a[(base + 2) * NCOLS + 3] = -y; // rx coefficient
        a[(base + 2) * NCOLS + 4] = x; // ry coefficient
                                       // rz = 0
        a[(base + 2) * NCOLS + 6] = z * 1e-6; // scale_ppm coefficient

        b_vec[base] = tgt[0] - src[0];
        b_vec[base + 1] = tgt[1] - src[1];
        b_vec[base + 2] = tgt[2] - src[2];
    }

    let x = solve_normal_equations(&a, &b_vec, rows, NCOLS)?;
    Ok((x[0], x[1], x[2], x[3], x[4], x[5], x[6]))
}

/// Solve `AᵀA·x = Aᵀb` via Cholesky decomposition (pure Rust).
///
/// `a` is stored row-major with shape `(m, n)`. `n` must be ≤ `m`.
/// Returns the solution vector of length `n`.
fn solve_normal_equations(a: &[f64], b: &[f64], m: usize, n: usize) -> Result<Vec<f64>, ProjError> {
    // Compute AᵀA (n×n, symmetric positive definite)
    let mut ata = vec![0.0f64; n * n];
    for i in 0..n {
        for j in 0..=i {
            let mut s = 0.0f64;
            for k in 0..m {
                s += a[k * n + i] * a[k * n + j];
            }
            ata[i * n + j] = s;
            ata[j * n + i] = s;
        }
    }

    // Compute Aᵀb (n)
    let mut atb = vec![0.0f64; n];
    for i in 0..n {
        let mut s = 0.0f64;
        for k in 0..m {
            s += a[k * n + i] * b[k];
        }
        atb[i] = s;
    }

    // Cholesky decompose AᵀA → L (lower triangular, in-place)
    let mut l = vec![0.0f64; n * n];
    for i in 0..n {
        for j in 0..=i {
            let mut s = ata[i * n + j];
            for k in 0..j {
                s -= l[i * n + k] * l[j * n + k];
            }
            if i == j {
                if s <= 0.0 {
                    return Err(ProjError::NoConvergence);
                }
                l[i * n + i] = s.sqrt();
            } else {
                l[i * n + j] = s / l[j * n + j];
            }
        }
    }

    // Forward substitution: L·y = Aᵀb
    let mut y = vec![0.0f64; n];
    for i in 0..n {
        let mut s = atb[i];
        for j in 0..i {
            s -= l[i * n + j] * y[j];
        }
        y[i] = s / l[i * n + i];
    }

    // Back substitution: Lᵀ·x = y
    let mut x = vec![0.0f64; n];
    for i in (0..n).rev() {
        let mut s = y[i];
        for j in (i + 1)..n {
            s -= l[j * n + i] * x[j];
        }
        x[i] = s / l[i * n + i];
    }

    Ok(x)
}

#[cfg(test)]
mod tests {
    use super::*;

    fn close(a: f64, b: f64, tol: f64, label: &str) {
        assert!(
            (a - b).abs() < tol,
            "{label}: got {a}, expected {b}, diff {}",
            (a - b).abs()
        );
    }

    fn params(
        tx: f64,
        ty: f64,
        tz: f64,
        rx: f64,
        ry: f64,
        rz: f64,
        scale_ppm: f64,
    ) -> Helmert7Params {
        Helmert7Params {
            tx,
            ty,
            tz,
            rx,
            ry,
            rz,
            scale_ppm,
        }
    }

    #[test]
    fn helmert7_batch_scalar_identity() {
        let pts = vec![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [-1.0, 0.0, 100.0]];
        let p = params(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
        let result = helmert7_batch_scalar(&pts, &p);
        for (orig, out) in pts.iter().zip(result.iter()) {
            for j in 0..3 {
                close(orig[j], out[j], 1e-10, "identity");
            }
        }
    }

    #[test]
    fn helmert7_batch_scalar_translation_only() {
        let pts = vec![[100.0, 200.0, 300.0]];
        let p = params(10.0, 20.0, 30.0, 0.0, 0.0, 0.0, 0.0);
        let result = helmert7_batch_scalar(&pts, &p);
        close(result[0][0], 110.0, 1e-9, "tx");
        close(result[0][1], 220.0, 1e-9, "ty");
        close(result[0][2], 330.0, 1e-9, "tz");
    }

    #[test]
    fn helmert7_batch_scalar_round_trip() {
        let pts: Vec<[f64; 3]> = (0..10)
            .map(|i| {
                [
                    1_000_000.0 + i as f64 * 50_000.0,
                    2_000_000.0 + i as f64 * 30_000.0,
                    3_000_000.0 + i as f64 * 20_000.0,
                ]
            })
            .collect();
        let p = params(50.0, -30.0, 100.0, 1e-6, -2e-6, 3e-6, 2.0);
        let fwd = helmert7_batch_scalar(&pts, &p);
        let inv = helmert7_batch_scalar_inv(&fwd, &p);
        for (orig, back) in pts.iter().zip(inv.iter()) {
            for j in 0..3 {
                close(orig[j], back[j], 0.01, "round_trip");
            }
        }
    }

    #[test]
    fn helmert7_batch_scalar_empty() {
        let p = params(1.0, 2.0, 3.0, 0.0, 0.0, 0.0, 0.0);
        let result = helmert7_batch_scalar(&[], &p);
        assert!(result.is_empty());
    }

    #[test]
    fn estimate_helmert7_lsq_insufficient_points() {
        let src = vec![[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]];
        let tgt = vec![[1.0, 0.0, 0.0], [2.0, 0.0, 0.0]];
        let result = estimate_helmert7_lsq(&src, &tgt);
        assert!(result.is_err(), "should fail with < 3 points");
    }

    #[test]
    fn estimate_helmert7_lsq_mismatched_lengths() {
        let src = vec![[0.0, 0.0, 0.0]];
        let tgt: Vec<[f64; 3]> = vec![];
        let result = estimate_helmert7_lsq(&src, &tgt);
        assert!(result.is_err(), "mismatched lengths should fail");
    }

    #[test]
    fn estimate_helmert7_lsq_translation_only() {
        // Pure translation: 4 points, no rotation, no scale
        let sources: Vec<[f64; 3]> = vec![
            [1_000_000.0, 0.0, 0.0],
            [0.0, 2_000_000.0, 0.0],
            [0.0, 0.0, 3_000_000.0],
            [500_000.0, 500_000.0, 500_000.0],
        ];
        let tx = 100.0;
        let ty = 200.0;
        let tz = -50.0;
        let targets: Vec<[f64; 3]> = sources
            .iter()
            .map(|&[x, y, z]| [x + tx, y + ty, z + tz])
            .collect();
        let (tx_est, ty_est, tz_est, rx_est, ry_est, rz_est, scale_est) =
            estimate_helmert7_lsq(&sources, &targets).expect("LSQ should succeed");
        close(tx_est, tx, 0.001, "tx_est");
        close(ty_est, ty, 0.001, "ty_est");
        close(tz_est, tz, 0.001, "tz_est");
        close(rx_est, 0.0, 1e-10, "rx_est");
        close(ry_est, 0.0, 1e-10, "ry_est");
        close(rz_est, 0.0, 1e-10, "rz_est");
        close(scale_est, 0.0, 1e-6, "scale_est");
    }

    #[test]
    fn estimate_helmert7_lsq_recovers_params() {
        let p = params(100.0, 200.0, -50.0, 1e-6, 2e-6, -1e-6, 0.5);

        // Use spatially distributed points (non-collinear, spanning all 3 axes)
        // to ensure the 7-parameter normal equations are well-conditioned.
        let sources: Vec<[f64; 3]> = vec![
            [4_000_000.0, 500_000.0, 4_800_000.0],
            [-2_000_000.0, 3_500_000.0, 4_200_000.0],
            [1_500_000.0, -1_000_000.0, 6_100_000.0],
            [3_200_000.0, 2_800_000.0, -600_000.0],
            [-800_000.0, -2_400_000.0, 5_500_000.0],
            [2_100_000.0, 4_100_000.0, 3_000_000.0],
            [-3_000_000.0, 1_200_000.0, 2_700_000.0],
            [700_000.0, -3_300_000.0, 1_800_000.0],
        ];
        let targets = helmert7_batch_scalar(&sources, &p);

        let (tx_est, ty_est, tz_est, rx_est, ry_est, rz_est, scale_est) =
            estimate_helmert7_lsq(&sources, &targets).expect("LSQ should succeed");

        close(tx_est, p.tx, 0.01, "tx_recovery");
        close(ty_est, p.ty, 0.01, "ty_recovery");
        close(tz_est, p.tz, 0.01, "tz_recovery");
        close(rx_est, p.rx, 1e-9, "rx_recovery");
        close(ry_est, p.ry, 1e-9, "ry_recovery");
        close(rz_est, p.rz, 1e-9, "rz_recovery");
        close(scale_est, p.scale_ppm, 1e-6, "scale_recovery");
    }

    #[cfg(feature = "blas")]
    #[test]
    fn helmert7_batch_blas_matches_scalar() {
        let pts: Vec<[f64; 3]> = (0..20)
            .map(|i| {
                [
                    1_000_000.0 + i as f64 * 50_000.0,
                    2_000_000.0 + i as f64 * 30_000.0,
                    3_000_000.0 + i as f64 * 20_000.0,
                ]
            })
            .collect();
        let p = params(100.0, -50.0, 30.0, 1e-6, -2e-6, 3e-6, 1.5);

        let scalar_result = helmert7_batch_scalar(&pts, &p);
        let blas_result = helmert7_batch_blas(&pts, &p);

        assert_eq!(scalar_result.len(), blas_result.len());
        for (s, b) in scalar_result.iter().zip(blas_result.iter()) {
            for j in 0..3 {
                close(s[j], b[j], 1e-6, "blas_vs_scalar");
            }
        }
    }

    #[cfg(feature = "blas")]
    #[test]
    fn helmert7_batch_blas_empty() {
        let p = params(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
        let result = helmert7_batch_blas(&[], &p);
        assert!(result.is_empty());
    }

    // ── GPU (CUDA) batch — exercises the real kernel when a device is present,
    //    otherwise the transparent scalar fallback (still correct). ──────────────
    #[cfg(feature = "gpu")]
    #[test]
    fn helmert7_batch_gpu_matches_scalar() {
        // 10k points spans many thread blocks (block = 256).
        let pts: Vec<[f64; 3]> = (0..10_000)
            .map(|i| {
                let f = i as f64;
                [
                    4_000_000.0 + f * 137.0,
                    500_000.0 - f * 53.0,
                    4_800_000.0 + f * 31.0,
                ]
            })
            .collect();
        let p = params(100.0, -50.0, 30.0, 1e-6, -2e-6, 3e-6, 1.5);

        let scalar = helmert7_batch_scalar(&pts, &p);
        let gpu = helmert7_batch_gpu(&pts, &p);

        assert_eq!(scalar.len(), gpu.len());
        let mut max_abs = 0.0f64;
        for (s, g) in scalar.iter().zip(gpu.iter()) {
            for j in 0..3 {
                max_abs = max_abs.max((s[j] - g[j]).abs());
            }
        }
        // Coordinates are metre-scale (~1e6 m); 1e-6 m = 1 µm is a few ULP and
        // accounts for the kernel's fused multiply-adds vs the scalar mul+add.
        assert!(
            max_abs < 1e-6,
            "GPU vs scalar Helmert max abs diff {max_abs} m too large"
        );
    }

    #[cfg(feature = "gpu")]
    #[test]
    fn helmert7_batch_gpu_empty() {
        let p = params(1.0, 2.0, 3.0, 0.0, 0.0, 0.0, 0.0);
        assert!(helmert7_batch_gpu(&[], &p).is_empty());
    }

    #[cfg(feature = "gpu")]
    #[test]
    fn helmert7_batch_gpu_round_trip() {
        let pts: Vec<[f64; 3]> = (0..2_000)
            .map(|i| {
                let f = i as f64;
                [1e6 + f * 500.0, 2e6 + f * 300.0, 3e6 + f * 200.0]
            })
            .collect();
        let p = params(50.0, -30.0, 100.0, 1e-6, -2e-6, 3e-6, 2.0);
        let fwd = helmert7_batch_gpu(&pts, &p);
        let back = helmert7_batch_scalar_inv(&fwd, &p);
        for (orig, b) in pts.iter().zip(back.iter()) {
            for j in 0..3 {
                close(orig[j], b[j], 0.01, "gpu_round_trip");
            }
        }
    }
}