petektools 0.1.0

//! `gridding` — interpolate scattered `(x, y, z)` samples onto a regular
//! [`Lattice`], producing a dense `Array2<f64>` of node values (`NaN` where a
//! node is undefined).
//!
//! This is the GATE-0 contract surface. The [`grid`] dispatcher and
//! [`GridMethod`] are **locked**; the three kernels (`nearest`, `idw`,
//! `min_curvature`) are ported from petekio 0.2.0's proven Briggs/IDW/nearest
//! implementations, held at behaviour parity with their source.
//!
//! Method set is intentionally an `enum` (not yet a trait), mirroring petekio's
//! `GridMethod` so delegation stays a 1:1 map. A `Gridder` trait for pluggable
//! kriging/RBF backends is a documented future extension (see `SPEC.md`).

mod convergent;
mod idw;
mod min_curvature;
mod nearest;

pub use convergent::ConvergentGridder;

use crate::foundation::{AlgoError, Lattice, Result};
use ndarray::Array2;

/// Scattered-data → grid interpolation methods.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum GridMethod {
    /// Value of the single areally-closest sample (blocky, exact at data).
    Nearest,
    /// Inverse-distance weighting (global, p = 2). Exact at coincident samples.
    InverseDistance,
    /// Briggs minimum-curvature (biharmonic) — smooth, honours the samples.
    MinimumCurvature,
}

/// Grid `coords` (`[x, y, z]` rows) onto `lattice` with `method`, returning the
/// `(ncol × nrow)` node-value array. Errors only on empty input; per-node
/// undefined values are `NaN`.
pub fn grid(coords: &[[f64; 3]], lattice: &Lattice, method: GridMethod) -> Result<Array2<f64>> {
    if coords.is_empty() {
        return Err(AlgoError::EmptyInput("grid: no points to grid"));
    }
    Ok(match method {
        GridMethod::Nearest => nearest::grid_nearest(coords, lattice),
        GridMethod::InverseDistance => idw::grid_idw(coords, lattice),
        GridMethod::MinimumCurvature => min_curvature::grid_min_curvature(coords, lattice, None),
    })
}

/// Warm-start minimum-curvature gridding: like [`grid`] with
/// [`GridMethod::MinimumCurvature`], but relax the SOR from `seed` (a
/// lattice-shaped prior field) instead of the cold IDW seed.
///
/// For an incremental re-grid (control points nudged, a point added) this
/// converges in far fewer iterations while reaching the same field — `warm ==
/// cold` to tolerance. A `None` or wrong-shape seed falls back to the cold
/// behaviour, so this is a non-breaking superset of `grid(.., MinimumCurvature)`.
/// The whole field is re-solved either way; the seed only sets the start point.
///
/// Errors only on empty input. Held at parity with petekio's seeded
/// `grid_min_curvature`.
pub fn grid_min_curvature_seeded(
    coords: &[[f64; 3]],
    lattice: &Lattice,
    seed: Option<&Array2<f64>>,
) -> Result<Array2<f64>> {
    if coords.is_empty() {
        return Err(AlgoError::EmptyInput(
            "grid_min_curvature_seeded: no points to grid",
        ));
    }
    Ok(min_curvature::grid_min_curvature(coords, lattice, seed))
}

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

    #[test]
    fn empty_input_errors() {
        let g = Lattice::regular(0.0, 0.0, 1.0, 1.0, 4, 4);
        assert!(matches!(
            grid(&[], &g, GridMethod::Nearest),
            Err(AlgoError::EmptyInput(_))
        ));
    }

    #[test]
    fn seeded_empty_input_errors() {
        let g = Lattice::regular(0.0, 0.0, 1.0, 1.0, 4, 4);
        assert!(matches!(
            grid_min_curvature_seeded(&[], &g, None),
            Err(AlgoError::EmptyInput(_))
        ));
    }

    #[test]
    fn seeded_none_equals_cold_grid() {
        // grid_min_curvature_seeded(.., None) must equal grid(.., MinimumCurvature):
        // the seeded entry is a non-breaking superset of the cold dispatch.
        let g = Lattice::regular(0.0, 0.0, 1.0, 1.0, 8, 7);
        let coords = [[1.0, 1.0, 3.0], [6.0, 5.0, 12.0], [3.0, 4.0, 7.0]];
        let cold = grid(&coords, &g, GridMethod::MinimumCurvature).unwrap();
        let seeded_none = grid_min_curvature_seeded(&coords, &g, None).unwrap();
        assert_eq!(cold, seeded_none);
        // Warm-start from the cold field reproduces it to the solver tolerance
        // (both stop at max_delta < TOL·range; 1e-3 is above that, well below
        // the field magnitude).
        let warm = grid_min_curvature_seeded(&coords, &g, Some(&cold)).unwrap();
        for (w, c) in warm.iter().zip(cold.iter()) {
            assert!((w - c).abs() < 1e-3);
        }
    }
}