diff --git a/src/geometry.rs b/src/geometry.rs index 26e1189..ae8f7c1 100644 --- a/src/geometry.rs +++ b/src/geometry.rs @@ -1,16 +1,30 @@ +//! Pure geometric primitives: norms, determinants, barycentric coordinates, +//! circumspheres, orientations, and simplex volumes. +//! +//! Nothing in this module touches Python or triangulation state, so every +//! function here is testable in isolation. Error messages mirror the Python +//! reference implementation in `adaptive.learner.triangulation`. + use nalgebra::{DMatrix, DVector}; use thiserror::Error; +use crate::tolerances::ORIENTATION_LOG_DET_CUTOFF; + +/// Errors produced by the geometric primitives. #[derive(Debug, Error)] pub enum GeometryError { + /// Input has the wrong shape (point/vertex counts or coordinate lengths). #[error("{0}")] InvalidDimensions(String), + /// The vertices span a lower-dimensional space than a simplex requires. #[error("Provided vertices do not form a simplex")] DegenerateSimplex, + /// A linear solve hit a (numerically) singular matrix. #[error("Singular matrix")] SingularMatrix, } +/// Euclidean norm, with unrolled 2D and 3D fast paths. #[inline] pub fn fast_norm(v: &[f64]) -> f64 { match v.len() { @@ -113,6 +127,9 @@ fn solve_square(matrix: &[Vec], rhs: &[f64]) -> Result, GeometryEr .ok_or(GeometryError::SingularMatrix) } +/// Rank of the matrix whose rows are `vectors`, counting singular values +/// above `tol`. A negative `tol` selects numpy's default tolerance +/// (`eps * max(rows, cols) * largest_singular_value`). pub fn matrix_rank(vectors: &[Vec], tol: f64) -> Result { if vectors.is_empty() { return Ok(0); @@ -139,10 +156,38 @@ pub fn matrix_rank(vectors: &[Vec], tol: f64) -> Result]) -> Result { matrix_rank(vectors, -1.0) } +/// Barycentric coordinates of `point` with respect to vertices `1..=dim` of +/// `simplex`; the coordinate of vertex 0 is `1 - sum(result)`. +/// +/// Accepts any slice of coordinate slices so callers can pass either owned +/// vertices or references into triangulation storage without cloning. +pub fn barycentric_coordinates>( + simplex: &[P], + point: &[f64], +) -> Result, GeometryError> { + let dim = point.len(); + let x0 = simplex[0].as_ref(); + let mut matrix = vec![vec![0.0; dim]; dim]; + let mut rhs = vec![0.0; dim]; + + for row in 0..dim { + rhs[row] = point[row] - x0[row]; + for col in 0..dim { + matrix[row][col] = simplex[col + 1].as_ref()[row] - x0[row]; + } + } + + solve_square(&matrix, &rhs) +} + +/// 2D point-in-triangle test via explicit barycentric formulas, with `eps` +/// slack on the barycentric coordinates (see +/// [`crate::tolerances::BARYCENTRIC_EPS`]). pub fn fast_2d_point_in_simplex( point: &[f64; 2], simplex: &[[f64; 2]; 3], @@ -165,6 +210,8 @@ pub fn fast_2d_point_in_simplex( Ok(t >= -eps && s + t <= 1.0 + eps) } +/// N-dimensional point-in-simplex test with `eps` slack on the barycentric +/// coordinates. Dispatches to [`fast_2d_point_in_simplex`] in 2D. pub fn point_in_simplex( point: &[f64], simplex: &[Vec], @@ -191,22 +238,12 @@ pub fn point_in_simplex( return fast_2d_point_in_simplex(&point, &simplex, eps); } - let x0 = &simplex[0]; - let dim = point.len(); - let mut matrix = vec![vec![0.0; dim]; dim]; - let mut rhs = vec![0.0; dim]; - - for row in 0..dim { - rhs[row] = point[row] - x0[row]; - for col in 0..dim { - matrix[row][col] = simplex[col + 1][row] - x0[row]; - } - } - - let alpha = solve_square(&matrix, &rhs)?; + let alpha = barycentric_coordinates(simplex, point)?; Ok(alpha.iter().all(|value| *value > -eps) && alpha.iter().sum::() < 1.0 + eps) } +/// Circumcircle (center, radius) of a triangle, computed in coordinates +/// relative to the first vertex for numerical stability. pub fn fast_2d_circumcircle(points: &[[f64; 2]; 3]) -> ([f64; 2], f64) { let [p0, p1, p2] = *points; let x1 = p1[0] - p0[0]; @@ -228,6 +265,8 @@ pub fn fast_2d_circumcircle(points: &[[f64; 2]; 3]) -> ([f64; 2], f64) { ([x + p0[0], y + p0[1]], radius) } +/// Circumsphere (center, radius) of a tetrahedron, computed in coordinates +/// relative to the first vertex for numerical stability. pub fn fast_3d_circumsphere(points: &[[f64; 3]; 4]) -> ([f64; 3], f64) { let [p0, p1, p2, p3] = *points; let x1 = p1[0] - p0[0]; @@ -258,6 +297,11 @@ pub fn fast_3d_circumsphere(points: &[[f64; 3]; 4]) -> ([f64; 3], f64) { ([cx + p0[0], cy + p0[1], cz + p0[2]], radius) } +/// Circumsphere (center, radius) of an N-dimensional simplex with `dim + 1` +/// vertices. Dispatches to the unrolled 2D/3D paths; the generic path solves +/// in coordinates relative to the first vertex (see PR #3 — the naive +/// `|x_i|^2 - |x_0|^2` form cancels catastrophically for small simplices far +/// from the origin). Returns NaNs for degenerate input, matching numpy. pub fn circumsphere(pts: &[Vec]) -> Result<(Vec, f64), GeometryError> { let dim = validate_points(pts)?; if pts.len() != dim + 1 { @@ -371,6 +415,10 @@ fn slogdet(matrix: &[Vec]) -> Result<(f64, f64), GeometryError> { Ok((sign, log_abs_det)) } +/// Sign (+1, -1, or 0) of the determinant of `face - origin`, i.e. which side +/// of the hyperplane through `face` the point `origin` lies on. Returns 0 +/// when the log-determinant falls below +/// [`ORIENTATION_LOG_DET_CUTOFF`]. pub fn orientation(face: &[Vec], origin: &[f64]) -> Result { let dim = validate_points(face)?; if face.len() != dim || origin.len() != dim { @@ -384,7 +432,7 @@ pub fn orientation(face: &[Vec], origin: &[f64]) -> Result], origin: &[f64]) -> Result]) -> Result { let dim = validate_points(vertices)?; if vertices.len() != dim + 1 { @@ -413,6 +463,10 @@ pub fn volume(vertices: &[Vec]) -> Result { Ok(determinant(&matrix)?.abs() / factorial(dim)) } +/// Volume of a `k`-simplex embedded in a higher-dimensional space (segment +/// length, triangle area via Heron's formula, Cayley-Menger determinant in +/// general). Unlike [`volume`], the number of vertices may be smaller than +/// `dim + 1`. pub fn simplex_volume_in_embedding(vertices: &[Vec]) -> Result { validate_points(vertices)?; if vertices.len() < 2 { diff --git a/src/lib.rs b/src/lib.rs index 7223c06..325a459 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,4 +1,10 @@ +//! Rust backend for `adaptive_triangulation`: a drop-in replacement for +//! `adaptive.learner.triangulation` (see [`triangulation`] for the class and +//! [`geometry`] for the free functions), with the numerical tolerance policy +//! documented in [`tolerances`]. + pub mod geometry; +pub mod tolerances; pub mod triangulation; use pyo3::exceptions::{PyValueError, PyZeroDivisionError}; @@ -7,21 +13,11 @@ use pyo3::types::{PyAny, PyModule, PyTuple}; use crate::geometry as geom; use crate::triangulation::{ - parse_point, parse_points_sized, point_tuple, PyFacesIter, PySimplicesProxy, PyTriangulation, - PyVertexToSimplicesIter, PyVertexToSimplicesProxy, PyVerticesIter, PyVerticesProxy, + numpy_linalg_error, parse_point, parse_points_sized, point_tuple, PyFacesIter, + PySimplicesProxy, PyTriangulation, PyVertexToSimplicesIter, PyVertexToSimplicesProxy, + PyVerticesIter, PyVerticesProxy, }; -fn numpy_linalg_error(py: Python<'_>, message: &str) -> PyErr { - PyModule::import(py, "numpy.linalg") - .and_then(|module| { - let error_type = module.getattr("LinAlgError")?; - let args = PyTuple::new(py, [message]).unwrap(); - let value = error_type.call1(args)?; - Ok(PyErr::from_value(value)) - }) - .unwrap_or_else(|_| PyValueError::new_err(message.to_string())) -} - fn point_in_simplex_error(py: Python<'_>, err: geom::GeometryError) -> PyErr { match err { geom::GeometryError::DegenerateSimplex => PyZeroDivisionError::new_err("division by zero"), @@ -126,7 +122,7 @@ fn py_point_in_simplex( ) -> PyResult { let point = parse_point(point)?; let simplex = parse_points_sized(simplex, "Please provide a 2-dimensional list of points")?; - geom::point_in_simplex(&point, &simplex, eps.unwrap_or(1e-8)) + geom::point_in_simplex(&point, &simplex, eps.unwrap_or(tolerances::BARYCENTRIC_EPS)) .map_err(|err| point_in_simplex_error(py, err)) } @@ -152,7 +148,7 @@ fn py_fast_2d_point_in_simplex( [simplex[1][0], simplex[1][1]], [simplex[2][0], simplex[2][1]], ], - eps.unwrap_or(1e-8), + eps.unwrap_or(tolerances::BARYCENTRIC_EPS), ) .map_err(|err| point_in_simplex_error(py, err)) } diff --git a/src/tolerances.rs b/src/tolerances.rs new file mode 100644 index 0000000..50b9e7b --- /dev/null +++ b/src/tolerances.rs @@ -0,0 +1,82 @@ +//! Numerical tolerance policy. +//! +//! Every fuzzy comparison in this crate goes through one of the constants +//! below. The values intentionally match the Python reference implementation +//! in `adaptive.learner.triangulation` so that both backends produce +//! identical triangulations from identical input; do not change one without +//! cross-validating against the reference (see `tests/test_triangulation.py`). +//! +//! When touching any of these, first check whether the comparison is +//! *scale-invariant* (relative to the simplex being tested) or *absolute* +//! (silently assumes coordinates of order 1). Mixing the two up is how small +//! simplices far from the origin get misclassified — see the regression tests +//! added for exactly that failure in PR #3. +//! +//! # Known limitation +//! +//! Point sets that mix widely separated coordinate scales in one +//! triangulation (e.g. a unit-sized cloud plus a 1e-5-sized cluster 100 +//! away) force sliver simplices with aspect ratios around 1e7 into the mesh. +//! Floating-point circumsphere tests on such slivers are unreliable, and the +//! Bowyer-Watson cascade can then assemble an inconsistent cavity; this is +//! caught loudly by the volume-conservation assertion rather than corrupting +//! state silently. The Python reference fails the same way (more often, in +//! fact). Fixing it would require exact geometric predicates, i.e. +//! adaptive-precision orientation/insphere tests, which is out of scope for +//! a tolerance constant. Homogeneous-scale inputs — at any magnitude and any +//! offset from the origin — are unaffected; see +//! `test_random_insertions_small_scale_far_from_origin`. + +/// Relative slack on barycentric coordinates when deciding whether a point +/// lies inside a simplex, on one of its faces, or coincides with a vertex. +/// +/// Scale-invariant: barycentric coordinates are already normalized by the +/// simplex size. Also exposed to Python as the default `eps` argument of +/// `point_in_simplex` and `get_reduced_simplex`. +pub const BARYCENTRIC_EPS: f64 = 1e-8; + +/// Relative slack on the circumsphere radius in the Bowyer-Watson +/// in-circumcircle test: a point is "inside" when +/// `distance < radius * (1.0 + CIRCUMCIRCLE_RTOL)`. +/// +/// Scale-invariant, but note it is *proportional to the neighbouring +/// simplex's size*, not the inserted point's distance to it — a point can be +/// within this slack of a large neighbour while being a meaningful fraction +/// of a small simplex away from it. Two guards in +/// `Triangulation::bowyer_watson` keep that asymmetry from corrupting the +/// triangulation: 1D insertion skips the circumcircle cascade entirely, and +/// in higher dimensions an insertion that would strand a cavity vertex is +/// rejected as a duplicate before anything is mutated. +pub const CIRCUMCIRCLE_RTOL: f64 = 1e-8; + +/// Threshold below which a candidate simplex is discarded as numerically +/// degenerate during Bowyer-Watson re-triangulation. A candidate is dropped +/// only when BOTH its raw volume (absolute; guarantees the omission cannot +/// leave a material hole) and its normalized volume (volume divided by +/// characteristic length to the dim-th power, i.e. scale-invariant +/// *flatness*) fall below this constant. +/// +/// The flatness condition is a deliberate divergence from the Python +/// reference, which uses the absolute test alone and therefore silently +/// deletes well-shaped simplices once a mesh is refined below ~1e-8 volume +/// (in 2D that is edge lengths of ~1e-4) — emptying whole regions of the +/// triangulation. See `Triangulation::simplex_is_numerically_degenerate`. +pub const DEGENERATE_VOLUME_EPS: f64 = 1e-8; + +/// Absolute and relative tolerances of the volume-conservation check after +/// Bowyer-Watson (the sum of deleted simplex volumes must equal the sum of +/// added ones). Same constants as `numpy.isclose`, but applied symmetrically +/// (scaled by `max(|a|, |b|)` rather than `|b|`) so the check cannot depend +/// on argument order. +pub const VOLUME_CONSERVATION_ATOL: f64 = 1e-8; +/// See [`VOLUME_CONSERVATION_ATOL`]. +pub const VOLUME_CONSERVATION_RTOL: f64 = 1e-5; + +/// `orientation` treats a face as coplanar with the origin point when the +/// log-determinant of the difference matrix is below this cutoff +/// (`|det| < exp(-50) ~= 2e-22`). +/// +/// NOT scale-invariant: the determinant scales with the dim-th power of the +/// face's edge lengths, so faces with edges below ~`exp(-50/dim)` are always +/// reported as degenerate. Matches the Python reference. +pub const ORIENTATION_LOG_DET_CUTOFF: f64 = -50.0; diff --git a/src/triangulation.rs b/src/triangulation.rs index 079bf2a..f2f3e9f 100644 --- a/src/triangulation.rs +++ b/src/triangulation.rs @@ -1,5 +1,16 @@ +//! N-dimensional Delaunay triangulation with incremental Bowyer-Watson +//! point insertion. +//! +//! The module has two layers: the pure-Rust [`Triangulation`] core, which +//! holds all state and algorithms and never touches Python, and the PyO3 +//! surface ([`PyTriangulation`] plus the proxy/iterator classes) which only +//! parses arguments, forwards to the core, and converts results and errors. +//! The Python API deliberately mirrors `adaptive.learner.triangulation` +//! (including its `point_in_cicumcircle` typo), so adaptive can use either +//! backend interchangeably. + use std::collections::VecDeque; -use std::sync::RwLock; +use std::sync::{PoisonError, RwLock}; use numpy::PyArray2; use pyo3::exceptions::{ @@ -12,22 +23,33 @@ use rustc_hash::{FxHashMap, FxHashSet}; use thiserror::Error; use crate::geometry::{self, GeometryError}; +use crate::tolerances::{ + BARYCENTRIC_EPS, CIRCUMCIRCLE_RTOL, DEGENERATE_VOLUME_EPS, VOLUME_CONSERVATION_ATOL, + VOLUME_CONSERVATION_RTOL, +}; +/// Index of a vertex in [`Triangulation::vertices`]. pub type PointIndex = usize; +/// A simplex as a sorted list of vertex indices. pub type Simplex = Vec; -const DEFAULT_EPS: f64 = 1e-8; - +/// Errors from triangulation operations. Each variant maps onto the Python +/// exception type the reference implementation raises in the same situation. #[derive(Debug, Error)] pub enum TriangulationError { + /// Maps to `ValueError`. #[error("{0}")] Value(String), + /// Maps to `IndexError`. #[error("{0}")] Index(String), + /// Maps to `RuntimeError`. #[error("{0}")] Runtime(String), + /// Maps to `AssertionError` (the reference uses bare `assert`). #[error("{0}")] Assertion(String), + /// Maps to `ValueError`. #[error(transparent)] Geometry(#[from] GeometryError), } @@ -39,10 +61,7 @@ impl TriangulationError { Self::Index(message) => PyIndexError::new_err(message), Self::Runtime(message) => PyRuntimeError::new_err(message), Self::Assertion(message) => PyAssertionError::new_err(message), - Self::Geometry(error) => match error { - GeometryError::SingularMatrix => PyValueError::new_err(error.to_string()), - _ => PyValueError::new_err(error.to_string()), - }, + Self::Geometry(error) => PyValueError::new_err(error.to_string()), } } } @@ -291,7 +310,10 @@ fn apply_transform(point: &[f64], transform: &[Vec]) -> Vec { result } -fn numpy_linalg_error(py: Python<'_>, message: &str) -> PyErr { +/// Raise `numpy.linalg.LinAlgError` (falling back to `ValueError` when numpy +/// is unavailable), matching where the reference implementation lets numpy +/// solves fail. +pub(crate) fn numpy_linalg_error(py: Python<'_>, message: &str) -> PyErr { PyModule::import(py, "numpy.linalg") .and_then(|module| { let error_type = module.getattr("LinAlgError")?; @@ -302,26 +324,26 @@ fn numpy_linalg_error(py: Python<'_>, message: &str) -> PyErr { .unwrap_or_else(|_| PyValueError::new_err(message.to_string())) } -fn barycentric_alpha(vertices: &[Vec], point: &[f64]) -> Result, TriangulationError> { - let dim = point.len(); - let x0 = &vertices[0]; - let mut matrix = vec![vec![0.0; dim]; dim]; - let mut rhs = vec![0.0; dim]; - - for row in 0..dim { - rhs[row] = point[row] - x0[row]; - for col in 0..dim { - matrix[row][col] = vertices[col + 1][row] - x0[row]; - } +/// Positions (`0..=dim`) of the simplex vertices kept after reducing against +/// the barycentric coordinates `alpha` (relative to vertex 0): vertices whose +/// coordinate exceeds `eps` survive, and an empty result means the point lies +/// outside the simplex. A single surviving position means the point coincides +/// with that vertex. +fn reduced_simplex_positions(alpha: &[f64], eps: f64) -> Vec { + let sum_alpha = alpha.iter().sum::(); + if alpha.iter().any(|value| *value < -eps) || sum_alpha > 1.0 + eps { + return Vec::new(); } - let flat: Vec = matrix.iter().flat_map(|row| row.iter().copied()).collect(); - let mat = nalgebra::DMatrix::from_row_slice(dim, dim, &flat); - let rhs = nalgebra::DVector::from_column_slice(&rhs); - mat.lu() - .solve(&rhs) - .map(|solution| solution.iter().copied().collect()) - .ok_or(TriangulationError::Geometry(GeometryError::SingularMatrix)) + let mut positions: Vec = alpha + .iter() + .enumerate() + .filter_map(|(idx, value)| (*value > eps).then_some(idx + 1)) + .collect(); + if sum_alpha < 1.0 - eps { + positions.insert(0, 0); + } + positions } fn combinations( @@ -380,15 +402,26 @@ fn is_close(a: f64, b: f64) -> bool { // Symmetric variant of numpy.isclose (which scales by |b| only), so the // volume-conservation check cannot depend on argument order. let scale = a.abs().max(b.abs()); - (a - b).abs() <= DEFAULT_EPS + 1e-5 * scale + (a - b).abs() <= VOLUME_CONSERVATION_ATOL + VOLUME_CONSERVATION_RTOL * scale } +/// An N-dimensional Delaunay triangulation supporting incremental point +/// insertion via the Bowyer-Watson algorithm. Pure Rust; the Python surface +/// lives on [`PyTriangulation`]. #[derive(Debug)] pub struct Triangulation { + /// Vertex coordinates, indexed by [`PointIndex`]. Vertices are never + /// removed; a vertex may belong to no simplex (e.g. rejected duplicates). pub vertices: Vec>, + /// All current simplices, each sorted ascending. pub simplices: FxHashSet, + /// For each vertex, the simplices it belongs to (inverse of + /// [`Self::simplices`]; see [`Self::reference_invariant`]). pub vertex_to_simplices: Vec>, + /// Dimensionality of the vertex coordinates. pub dim: usize, + // Cache of the simplex found by the last point location, used as the + // start of the next walk. Interior-mutable so lookups stay `&self`. last_simplex: RwLock>, } @@ -399,12 +432,28 @@ impl Clone for Triangulation { simplices: self.simplices.clone(), vertex_to_simplices: self.vertex_to_simplices.clone(), dim: self.dim, - last_simplex: RwLock::new(self.last_simplex.read().unwrap().clone()), + last_simplex: RwLock::new(self.last_simplex().clone()), } } } impl Triangulation { + // The location cache is only ever replaced wholesale, so a poisoned lock + // (a panic while it was held) cannot leave it inconsistent; recover the + // value rather than propagating the panic across the FFI boundary. + fn last_simplex(&self) -> std::sync::RwLockReadGuard<'_, Option> { + self.last_simplex + .read() + .unwrap_or_else(PoisonError::into_inner) + } + + fn set_last_simplex(&self, simplex: Option) { + *self + .last_simplex + .write() + .unwrap_or_else(PoisonError::into_inner) = simplex; + } + fn validate_coords(coords: &[Vec]) -> Result { if coords.is_empty() { return Err(TriangulationError::Value( @@ -481,6 +530,8 @@ impl Triangulation { )) } + /// Build a triangulation from known simplices (e.g. scipy's Delaunay + /// output) without running incremental insertion. pub fn from_simplices( coords: Vec>, simplices: impl IntoIterator, @@ -527,6 +578,9 @@ impl Triangulation { Self::from_simplices(coords, segments) } + /// Build a triangulation from scratch: sorted adjacency in 1D, otherwise + /// a seed simplex followed by incremental insertion of the remaining + /// points. Requires at least `dim + 1` points spanning full rank. pub fn new(coords: Vec>) -> Result { let dim = Self::validate_coords(&coords)?; if dim == 1 { @@ -560,7 +614,7 @@ impl Triangulation { } let reduced_simplex = - triangulation.get_reduced_simplex(&point, &actual_simplex, DEFAULT_EPS)?; + triangulation.get_reduced_simplex(&point, &actual_simplex, BARYCENTRIC_EPS)?; if reduced_simplex.is_empty() { return Err(TriangulationError::Value( "Point lies outside of the specified simplex.".to_string(), @@ -575,6 +629,8 @@ impl Triangulation { Ok(triangulation) } + /// Insert a simplex (sorted automatically), keeping + /// [`Self::vertex_to_simplices`] in sync. No geometric checks are done. pub fn add_simplex(&mut self, mut simplex: Simplex) -> Result<(), TriangulationError> { simplex.sort_unstable(); if simplex.len() != self.dim + 1 { @@ -592,6 +648,7 @@ impl Triangulation { Ok(()) } + /// Remove a simplex, keeping [`Self::vertex_to_simplices`] in sync. pub fn delete_simplex(&mut self, simplex: &[usize]) -> Result<(), TriangulationError> { let mut simplex = simplex.to_vec(); simplex.sort_unstable(); @@ -604,6 +661,7 @@ impl Triangulation { Ok(()) } + /// Coordinates of the given vertices, cloned in order. pub fn get_vertices( &self, indices: &[PointIndex], @@ -618,38 +676,26 @@ impl Triangulation { fn locate_point_scan(&self, point: &[f64]) -> Result, TriangulationError> { for simplex in &self.simplices { let vertices = self.get_vertices(simplex)?; - if geometry::point_in_simplex(point, &vertices, DEFAULT_EPS)? { - *self.last_simplex.write().unwrap() = Some(simplex.clone()); + if geometry::point_in_simplex(point, &vertices, BARYCENTRIC_EPS)? { + self.set_last_simplex(Some(simplex.clone())); return Ok(Some(simplex.clone())); } } Ok(None) } + /// Barycentric coordinates of `point` relative to vertices `1..=dim` of + /// `simplex`, without cloning the vertex coordinates. fn barycentric_alpha_for_simplex( &self, simplex: &[usize], point: &[f64], ) -> Result, TriangulationError> { - let dim = point.len(); - let x0 = &self.vertices[simplex[0]]; - let mut matrix = vec![vec![0.0; dim]; dim]; - let mut rhs = vec![0.0; dim]; - - for row in 0..dim { - rhs[row] = point[row] - x0[row]; - for col in 0..dim { - matrix[row][col] = self.vertices[simplex[col + 1]][row] - x0[row]; - } - } - - let flat: Vec = matrix.iter().flat_map(|row| row.iter().copied()).collect(); - let mat = nalgebra::DMatrix::from_row_slice(dim, dim, &flat); - let rhs = nalgebra::DVector::from_column_slice(&rhs); - mat.lu() - .solve(&rhs) - .map(|solution| solution.iter().copied().collect()) - .ok_or(TriangulationError::Geometry(GeometryError::SingularMatrix)) + let vertices: Vec<&[f64]> = simplex + .iter() + .map(|&index| self.vertices[index].as_slice()) + .collect(); + geometry::barycentric_coordinates(&vertices, point).map_err(TriangulationError::Geometry) } fn next_simplex_in_walk( @@ -669,8 +715,8 @@ impl Triangulation { } } - if worst_value >= -DEFAULT_EPS { - *self.last_simplex.write().unwrap() = Some(simplex.to_vec()); + if worst_value >= -BARYCENTRIC_EPS { + self.set_last_simplex(Some(simplex.to_vec())); return Ok(Some(simplex.to_vec())); } @@ -686,12 +732,13 @@ impl Triangulation { Ok(neighbours.into_iter().next()) } + /// Find a simplex containing `point` by walking from the last located + /// simplex, falling back to a linear scan if the walk hits a degenerate + /// solve. Returns `None` when the point is outside the hull. pub fn locate_point(&self, point: &[f64]) -> Result, TriangulationError> { self.validate_point_dim(point)?; let Some(mut current) = self - .last_simplex - .read() - .unwrap() + .last_simplex() .clone() .filter(|simplex| self.simplices.contains(simplex)) .or_else(|| self.simplices.iter().next().cloned()) @@ -713,6 +760,10 @@ impl Triangulation { self.locate_point_scan(point) } + /// Reduce `simplex` to the face of it that `point` actually lies on + /// (with `eps` barycentric slack): the full simplex for an interior + /// point, a lower-dimensional face for a point on the boundary, a single + /// vertex for a coincident point, or empty when the point lies outside. pub fn get_reduced_simplex( &self, point: &[f64], @@ -732,23 +783,14 @@ impl Triangulation { }; let alpha = self.barycentric_alpha_for_simplex(&simplex, point)?; - let sum_alpha = alpha.iter().sum::(); - - if alpha.iter().any(|value| *value < -eps) || sum_alpha > 1.0 + eps { - return Ok(Vec::new()); - } - - let mut result: Vec = alpha - .iter() - .enumerate() - .filter_map(|(idx, value)| (*value > eps).then_some(simplex[idx + 1])) - .collect(); - if sum_alpha < 1.0 - eps { - result.insert(0, simplex[0]); - } - Ok(result) + Ok(reduced_simplex_positions(&alpha, eps) + .into_iter() + .map(|position| simplex[position]) + .collect()) } + /// Whether `point` lies inside the simplex given by vertex `indices`, + /// with `eps` barycentric slack. pub fn point_in_simplex( &self, point: &[f64], @@ -759,6 +801,10 @@ impl Triangulation { Ok(geometry::point_in_simplex(point, &vertices, eps)?) } + /// All `dim`-vertex faces (default: (dim-1)-faces, i.e. facets) of either + /// the whole triangulation, the given `simplices`, or the simplices + /// touching the given `vertices` (at most one of the two filters). Faces + /// shared by several simplices are repeated. pub fn faces( &self, dim: Option, @@ -804,6 +850,23 @@ impl Triangulation { } } + /// Count, for every (dim-1)-face of the given simplex pool (or of the + /// whole triangulation when `None`), how many simplices it belongs to. + /// In a consistent triangulation a count of 1 means a boundary face and + /// 2 an interior face. + fn face_multiplicities( + &self, + simplices: Option<&FxHashSet>, + ) -> Result, TriangulationError> { + let faces = self.faces(None, simplices, None)?; + let mut multiplicities: FxHashMap = FxHashMap::default(); + for face in faces { + *multiplicities.entry(face).or_insert(0) += 1; + } + Ok(multiplicities) + } + + /// All simplices that contain every vertex of `face`. pub fn containing(&self, face: &[usize]) -> Result, TriangulationError> { if face.is_empty() { return Ok(FxHashSet::default()); @@ -848,6 +911,8 @@ impl Triangulation { Ok(geometry::circumsphere(&points)?) } + /// Circumsphere (center, radius) of a simplex, optionally after applying + /// a linear `transform` to its vertices. pub fn circumscribed_circle( &self, simplex: &[usize], @@ -875,9 +940,12 @@ impl Triangulation { .map(|(a, b)| a - b) .collect::>(), ); - Ok(distance < radius * (1.0 + DEFAULT_EPS)) + Ok(distance < radius * (1.0 + CIRCUMCIRCLE_RTOL)) } + /// Whether vertex `pt_index` lies inside the (optionally transformed) + /// circumsphere of `simplex`, with + /// [`CIRCUMCIRCLE_RTOL`] slack on the radius. pub fn point_in_circumcircle( &self, pt_index: usize, @@ -896,6 +964,12 @@ impl Triangulation { self.point_in_circumcircle_impl(&self.vertices[pt_index], simplex, transform.as_deref()) } + /// Re-triangulate around vertex `pt_index`: delete every simplex whose + /// circumsphere contains it (cascading through facet neighbours, except + /// in 1D) and fill the resulting hole with simplices through the vertex. + /// Returns `(deleted, created)` and verifies total volume is conserved. + /// Errors without mutating anything when the insertion would strand a + /// cavity vertex, i.e. the point is numerically a duplicate of it. pub fn bowyer_watson( &mut self, pt_index: usize, @@ -908,10 +982,12 @@ impl Triangulation { self.validate_simplex_indices(simplex)?; } + // Collect the cavity (every simplex whose circumsphere contains the + // point, cascading through facet neighbours) without mutating, so + // that pathological insertions can still be rejected cleanly below. let mut queue = VecDeque::new(); let mut queued = FxHashSet::default(); - let mut done_simplices = FxHashSet::default(); - let mut bad_triangles = FxHashSet::default(); + let mut bad_triangles: FxHashSet = FxHashSet::default(); if let Some(simplex) = containing_simplex { queued.insert(simplex.clone()); @@ -925,23 +1001,16 @@ impl Triangulation { } while let Some(simplex) = queue.pop_front() { - if !done_simplices.insert(simplex.clone()) { - continue; - } if !self.simplices.contains(&simplex) { continue; } if self.point_in_circumcircle(pt_index, &simplex, transform)? { - self.delete_simplex(&simplex)?; bad_triangles.insert(simplex.clone()); if self.dim == 1 { // Inserting a point into an interval never invalidates // neighbouring intervals, so there is no cascade in 1D. - // Cascading on the (1 + eps) circumcircle tolerance can - // delete a neighbour the point is strictly outside of, - // orphaning the shared vertex. continue; } @@ -952,9 +1021,6 @@ impl Triangulation { } for neighbour in neighbours { - if done_simplices.contains(&neighbour) { - continue; - } let shared = neighbour .iter() .filter(|vertex| simplex_vertices.contains(vertex)) @@ -966,16 +1032,13 @@ impl Triangulation { } } - let faces = self.faces(None, Some(&bad_triangles), None)?; - let mut multiplicities: FxHashMap = FxHashMap::default(); - for face in &faces { - *multiplicities.entry(face.clone()).or_insert(0) += 1; - } - let hole_faces: Vec = faces + let hole_faces: Vec = self + .face_multiplicities(Some(&bad_triangles))? .into_iter() - .filter(|face| multiplicities.get(face).copied().unwrap_or_default() < 2) + .filter_map(|(face, count)| (count == 1).then_some(face)) .collect(); + let mut candidates: Vec = Vec::new(); for face in hole_faces { if face.contains(&pt_index) { continue; @@ -987,6 +1050,33 @@ impl Triangulation { if self.simplex_is_numerically_degenerate(&simplex)? { continue; } + candidates.push(simplex); + } + + // A cavity vertex that would lose all of its simplices without + // being reconnected means the point is numerically indistinguishable + // from that vertex (it sits within the circumcircle slack of every + // simplex around it). Reject the insertion before mutating anything + // rather than orphaning the vertex and corrupting the triangulation. + let candidate_vertices: FxHashSet = candidates.iter().flatten().copied().collect(); + for &vertex in bad_triangles.iter().flatten() { + if vertex == pt_index || candidate_vertices.contains(&vertex) { + continue; + } + if self.vertex_to_simplices[vertex] + .iter() + .all(|simplex| bad_triangles.contains(simplex)) + { + return Err(TriangulationError::Value( + "Point already in triangulation.".to_string(), + )); + } + } + + for simplex in &bad_triangles { + self.delete_simplex(simplex)?; + } + for simplex in candidates { self.add_simplex(simplex)?; } @@ -1011,17 +1101,16 @@ impl Triangulation { Ok((deleted_simplices, new_simplices)) } + /// Connect vertex `pt_index`, which must lie outside the convex hull, to + /// every hull face that faces it. Returns the created simplices; errors + /// (and rolls back) if the vertex turns out to be inside the hull. pub fn extend_hull( &mut self, pt_index: usize, ) -> Result, TriangulationError> { self.validate_vertex_index(pt_index)?; - let faces = self.faces(None, None, None)?; - let mut multiplicities: FxHashMap = FxHashMap::default(); - for face in &faces { - *multiplicities.entry(face.clone()).or_insert(0) += 1; - } - let hull_faces: Vec = multiplicities + let hull_faces: Vec = self + .face_multiplicities(None)? .into_iter() .filter_map(|(face, count)| (count == 1).then_some(face)) .collect(); @@ -1072,6 +1161,10 @@ impl Triangulation { Ok(new_simplices) } + /// Insert a new point, locating it first unless a containing `simplex` + /// is supplied. Returns the `(deleted, created)` simplex sets. Errors + /// without modifying the triangulation when the point is already present + /// or lies outside the supplied simplex. pub fn add_point( &mut self, point: Vec, @@ -1100,7 +1193,20 @@ impl Triangulation { } }; let (deleted_simplices, added_simplices) = - self.bowyer_watson(pt_index, None, &transform)?; + match self.bowyer_watson(pt_index, None, &transform) { + Ok(result) => result, + // Value errors are raised before bowyer_watson mutates, + // so only the hull extension needs rolling back. + Err(err @ TriangulationError::Value(_)) => { + for simplex in self.vertex_to_simplices[pt_index].clone() { + self.delete_simplex(&simplex)?; + } + self.vertex_to_simplices.pop(); + self.vertices.pop(); + return Err(err); + } + Err(err) => return Err(err), + }; let deleted: FxHashSet = deleted_simplices .difference(&temporary_simplices) @@ -1113,7 +1219,7 @@ impl Triangulation { return Ok((deleted, added)); } - let reduced_simplex = self.get_reduced_simplex(&point, &simplex, DEFAULT_EPS)?; + let reduced_simplex = self.get_reduced_simplex(&point, &simplex, BARYCENTRIC_EPS)?; if reduced_simplex.is_empty() { self.vertex_to_simplices.pop(); return Err(TriangulationError::Value( @@ -1131,9 +1237,19 @@ impl Triangulation { let pt_index = self.vertices.len(); self.vertices.push(point); - self.bowyer_watson(pt_index, Some(actual_simplex), &transform) + match self.bowyer_watson(pt_index, Some(actual_simplex), &transform) { + // Value errors are raised before bowyer_watson mutates, so only + // the freshly pushed vertex needs rolling back. + Err(err @ TriangulationError::Value(_)) => { + self.vertex_to_simplices.pop(); + self.vertices.pop(); + Err(err) + } + other => other, + } } + /// Volume of the simplex with the given vertex indices. pub fn volume(&self, simplex: &[usize]) -> Result { Ok(geometry::volume(&self.get_vertices(simplex)?)?) } @@ -1165,13 +1281,19 @@ impl Triangulation { &self, simplex: &[usize], ) -> Result { - if self.dim == 1 { - // In 1D we only want to reject coincident endpoints, not tiny but valid intervals. - return Ok(self.normalized_volume(simplex)? < DEFAULT_EPS); - } - Ok(self.volume(simplex)? < DEFAULT_EPS) - } - + // Drop a candidate only when its volume is negligible in BOTH + // senses: absolutely tiny, so removing it cannot leave a material + // hole (this is the Python reference's criterion), AND flat relative + // to its own extent. Either test alone is wrong in some regime: the + // absolute cutoff alone empties finely refined or small-coordinate + // meshes whose well-shaped simplices all sit below it, while the + // flatness cutoff alone deletes large-but-flat simplices whose + // volume the cavity genuinely needs (breaking volume conservation). + Ok(self.volume(simplex)? < DEGENERATE_VOLUME_EPS + && self.normalized_volume(simplex)? < DEGENERATE_VOLUME_EPS) + } + + /// Whether the (order-insensitive) simplex is present. pub fn has_simplex(&self, simplex: &[usize]) -> Result { let mut simplex = simplex.to_vec(); simplex.sort_unstable(); @@ -1179,6 +1301,7 @@ impl Triangulation { Ok(self.simplices.contains(&simplex)) } + /// The simplices containing `vertex`, by reference. pub fn vertex_to_simplices_for( &self, vertex: usize, @@ -1187,6 +1310,7 @@ impl Triangulation { Ok(&self.vertex_to_simplices[vertex]) } + /// Volumes of all simplices, in iteration order of [`Self::simplices`]. pub fn volumes(&self) -> Result, TriangulationError> { self.simplices .iter() @@ -1194,6 +1318,8 @@ impl Triangulation { .collect() } + /// Whether [`Self::simplices`] and [`Self::vertex_to_simplices`] are + /// mutually consistent (a bookkeeping check, not a geometric one). pub fn reference_invariant(&self) -> bool { for vertex in 0..self.vertices.len() { if self.vertex_to_simplices[vertex] @@ -1214,18 +1340,16 @@ impl Triangulation { true } + /// Vertices on the convex hull (those belonging to a boundary face). + /// Errors when a face belongs to more than two simplices, which means + /// the triangulation is internally inconsistent. pub fn hull(&self) -> Result, TriangulationError> { - let faces = self.faces(None, None, None)?; - let mut counts: FxHashMap = FxHashMap::default(); - for face in faces { - let count = counts.entry(face).or_insert(0); - *count += 1; - if *count > 2 { - return Err(TriangulationError::Runtime( - "Broken triangulation, a (N-1)-dimensional appears in more than 2 simplices." - .to_string(), - )); - } + let counts = self.face_multiplicities(None)?; + if counts.values().any(|&count| count > 2) { + return Err(TriangulationError::Runtime( + "Broken triangulation, a (N-1)-dimensional appears in more than 2 simplices." + .to_string(), + )); } let mut hull = FxHashSet::default(); @@ -1238,11 +1362,16 @@ impl Triangulation { } } +/// Python-facing `Triangulation`, a thin argument-parsing wrapper around the +/// pure-Rust [`Triangulation`] core. Drop-in compatible with +/// `adaptive.learner.triangulation.Triangulation`. #[pyclass(name = "Triangulation")] pub struct PyTriangulation { pub core: Triangulation, } +/// Lazy, set-like view of the simplices (all of them, or those of one +/// vertex) that reads triangulation state on access instead of copying it. #[pyclass(name = "SimplicesProxy")] pub struct PySimplicesProxy { triangulation: Py, @@ -1265,6 +1394,8 @@ impl PySimplicesProxy { } } +/// Lazy, sequence-like view of the vertex coordinates (supports `__array__` +/// for zero-surprise numpy conversion). #[pyclass(name = "VerticesProxy")] pub struct PyVerticesProxy { triangulation: Py, @@ -1276,6 +1407,7 @@ impl PyVerticesProxy { } } +/// Lazy, sequence-like view mapping each vertex index to its simplex set. #[pyclass(name = "VertexToSimplicesProxy")] pub struct PyVertexToSimplicesProxy { triangulation: Py, @@ -1287,6 +1419,7 @@ impl PyVertexToSimplicesProxy { } } +/// Iterator over a snapshot of faces/simplices as Python tuples. #[pyclass] pub struct PyFacesIter { items: Vec, @@ -1464,8 +1597,8 @@ impl PyTriangulation { .map_err(TriangulationError::into_pyerr)?; let core = match scipy_delaunay_simplices(py, &parsed_coords, dim)? { - Some(initial) => Triangulation::from_simplices(parsed_coords.clone(), initial), - None => Triangulation::new(parsed_coords.clone()), + Some(initial) => Triangulation::from_simplices(parsed_coords, initial), + None => Triangulation::new(parsed_coords), } .map_err(TriangulationError::into_pyerr)?; Ok(Self { core }) @@ -1585,7 +1718,7 @@ impl PyTriangulation { } } - #[pyo3(signature = (point, simplex, eps=DEFAULT_EPS))] + #[pyo3(signature = (point, simplex, eps=BARYCENTRIC_EPS))] fn get_reduced_simplex( &self, py: Python<'_>, @@ -1595,37 +1728,25 @@ impl PyTriangulation { ) -> PyResult> { let point = parse_point(point)?; let simplex_signed = parse_signed_indices(simplex)?; + // For a full simplex the reference echoes back the caller's indices + // (including negative ones), so reduce positions rather than values. if simplex_signed.len() == self.core.dim + 1 { self.core .validate_point_dim(&point) .map_err(TriangulationError::into_pyerr)?; let simplex = normalize_indices(&simplex_signed, self.core.vertices.len()) .map_err(TriangulationError::into_pyerr)?; - let vertices = self - .core - .get_vertices(&simplex) - .map_err(TriangulationError::into_pyerr)?; - let alpha = match barycentric_alpha(&vertices, &point) { + let alpha = match self.core.barycentric_alpha_for_simplex(&simplex, &point) { Ok(alpha) => alpha, Err(TriangulationError::Geometry(GeometryError::SingularMatrix)) => { return Err(numpy_linalg_error(py, "Singular matrix")); } Err(other) => return Err(other.into_pyerr()), }; - let sum_alpha = alpha.iter().sum::(); - if alpha.iter().any(|value| *value < -eps) || sum_alpha > 1.0 + eps { - return signed_index_list_py(py, &[]); - } - - let mut reduced = Vec::new(); - for (idx, value) in alpha.iter().enumerate() { - if *value > eps { - reduced.push(simplex_signed[idx + 1]); - } - } - if sum_alpha < 1.0 - eps { - reduced.insert(0, simplex_signed[0]); - } + let reduced: Vec = reduced_simplex_positions(&alpha, eps) + .into_iter() + .map(|position| simplex_signed[position]) + .collect(); return signed_index_list_py(py, &reduced); } @@ -1739,7 +1860,7 @@ impl PyTriangulation { Ok(PyArray2::from_vec2(py, &identity)?.into()) } - #[pyo3(signature = (point, simplex, eps=DEFAULT_EPS))] + #[pyo3(signature = (point, simplex, eps=BARYCENTRIC_EPS))] fn point_in_simplex( &self, py: Python<'_>, @@ -1876,6 +1997,55 @@ mod tests { assert_eq!(added, FxHashSet::from_iter([vec![1, 2]])); } + #[test] + fn two_dimensional_small_scale_insertion_keeps_mesh() { + // Regression: the absolute degenerate-volume cutoff dropped every + // candidate simplex below 1e-8 volume, emptying the mesh. + let mut tri = Triangulation::from_simplices( + vec![ + vec![0.0, 0.0], + vec![1e-4, 0.0], + vec![0.0, 1e-4], + vec![1e-4, 1e-4], + ], + vec![vec![0, 1, 3], vec![0, 2, 3]], + ) + .unwrap(); + + let (deleted, added) = tri.add_point(vec![5e-5, 6e-5], None, None).unwrap(); + + assert_eq!(deleted.len(), 2); + assert_eq!(added.len(), 4); + assert_eq!(tri.simplices.len(), 4); + assert!((0..tri.vertices.len()).all(|v| !tri.vertex_to_simplices[v].is_empty())); + } + + #[test] + fn two_dimensional_near_duplicate_insertion_is_rejected_without_mutation() { + // Regression: a point within the circumcircle slack of every simplex + // around an existing vertex orphaned that vertex. + let mut tri = Triangulation::new(vec![ + vec![0.0, 0.0], + vec![1.0, 0.0], + vec![0.5, 1.0], + vec![0.5, 0.5], + vec![0.5 + 1e-6, 0.5], + vec![0.5, 0.5 + 1e-6], + ]) + .unwrap(); + let simplices_before = tri.simplices.clone(); + let n_vertices_before = tri.vertices.len(); + + let err = tri + .add_point(vec![0.5 + 2e-10, 0.5 + 2e-10], None, None) + .unwrap_err(); + + assert!(matches!(err, TriangulationError::Value(_))); + assert_eq!(tri.simplices, simplices_before); + assert_eq!(tri.vertices.len(), n_vertices_before); + assert!(tri.reference_invariant()); + } + #[test] fn one_dimensional_insertion_near_shared_vertex_keeps_vertex_connected() { // Regression: the circumcircle cascade deleted the long neighbouring diff --git a/tests/test_triangulation.py b/tests/test_triangulation.py index ab9ba48..f27c58a 100644 --- a/tests/test_triangulation.py +++ b/tests/test_triangulation.py @@ -277,6 +277,70 @@ def test_add_point_in_tiny_interval_far_from_origin_1d(): assert_1d_triangulation_state(tri) +def all_vertices_connected(tri) -> bool: + return all(len(tri.vertex_to_simplices[i]) > 0 for i in range(len(tri.vertices))) + + +def test_add_point_in_small_scale_mesh_2d_keeps_all_simplices(): + # Regression: the absolute degenerate-volume cutoff (matching the Python + # reference) deleted every candidate simplex below 1e-8 volume, silently + # emptying the mesh on insertion into a finely refined region. + tri = rust_tri.Triangulation([[0, 0], [1e-4, 0], [0, 1e-4], [1e-4, 1e-4]]) + + deleted, added = tri.add_point([5e-5, 6e-5]) + + assert len(deleted) == 2 + assert len(added) == 4 + assert len(tri.simplices) == 4 + assert all_vertices_connected(tri) + assert tri.reference_invariant() + + +def test_near_duplicate_insertion_2d_is_rejected_not_orphaned(): + # Regression: a point within the circumcircle slack of every simplex + # around an existing vertex deleted them all and orphaned that vertex. + # It must instead be rejected as a duplicate, leaving the state intact. + pts = [ + [0.0, 0.0], + [1.0, 0.0], + [0.5, 1.0], + [0.5, 0.5], + [0.5 + 1e-6, 0.5], + [0.5, 0.5 + 1e-6], + ] + tri = rust_tri.Triangulation(pts) + before = as_simplex_set(tri.simplices) + + with pytest.raises(ValueError, match="Point already in triangulation"): + tri.add_point([0.5 + 2e-10, 0.5 + 2e-10]) + + assert as_simplex_set(tri.simplices) == before + assert len(tri.vertices) == len(pts) + assert all_vertices_connected(tri) + assert tri.reference_invariant() + + +@pytest.mark.parametrize("dim", [2, 3]) +@pytest.mark.parametrize( + ("offset", "scale"), + [(100.0, 1e-5), (-1e6, 1e-2), (0.5, 1e-6)], + ids=["near-100", "near-minus-1e6", "near-half"], +) +def test_random_insertions_small_scale_far_from_origin(dim, offset, scale): + rng = np.random.default_rng(7) + pts = offset + scale * rng.random((30, dim)) + tri = rust_tri.Triangulation(pts[: dim + 2].tolist()) + + for p in pts[dim + 2 :]: + tri.add_point(p.tolist()) + + assert all_vertices_connected(tri) + assert tri.reference_invariant() + assert all(v > 0 for v in tri.volumes()) + interior = pts.mean(axis=0) + assert len(tri.locate_point(interior.tolist())) == dim + 1 + + def test_add_point_near_shared_vertex_keeps_vertex_connected_1d(): # Regression: the circumcircle cascade also deleted the long neighbouring # interval the point is strictly outside of, orphaning the shared vertex.