Add 1D triangulation support (#3)

* feat: add 1d rust triangulation support

* refactor: clarify 1d simplex filtering

* refactor: consolidate 1d triangulation tests

* docs: summarize issue-095 cleanup

* Remove cleanup reports and fix pre-commit

* fix: make 1D triangulation numerically robust

Three related fixes for failures found on small intervals away from the
origin (e.g. Triangulation(np.linspace(100, 100.001, 50).reshape(-1, 1))
crashed with an AssertionError from the volume-conservation check):

- circumsphere: solve in coordinates relative to the first vertex. The
  generic path computed |x_i|^2 - |x_0|^2, which cancels catastrophically
  when a simplex is small compared to its distance from the origin
  (center error ~1e-11 for a 1e-6 interval at x=0.5, vs the 5e-15
  circumcircle tolerance). 1D always uses this path; dims >= 4 had the
  same latent inaccuracy.
- point_in_circumcircle: short-circuit to true when the point is a vertex
  of the simplex, so rounding in the circumcenter cannot flip the answer
  bowyer_watson relies on after extend_hull.
- bowyer_watson: skip the circumcircle cascade in 1D. Splitting an
  interval never invalidates neighbouring intervals, and cascading on the
  (1 + eps) tolerance could delete a neighbour the point is strictly
  outside of, orphaning the shared vertex and leaving overlapping
  segments.

Also fixes clippy lints from newer toolchains and renames a test that now
asserts a deliberate divergence from the reference.

* perf: construct 1D triangulations from sorted adjacency

1D Delaunay triangulation is just the chain of adjacent points in sorted
order, so build it directly instead of running O(n^2) incremental
insertion (full-face scans in extend_hull made constructing 2000 sorted
points take ~860 ms; now 50k points build in ~17 ms). Exact duplicates
keep the first occurrence, matching the incremental path's behaviour.

Author: basnijholt

src/geometry.rs

@@ -288,34 +288,31 @@ pub fn circumsphere(pts: &[Vec<f64>]) -> Result<(Vec<f64>, f64), GeometryError>
         return Ok((center.into_iter().collect(), radius));
     }
 
+    // Solve in coordinates relative to the first vertex: computing
+    // |x_i|^2 - |x_0|^2 directly cancels catastrophically when the simplex
+    // is small compared to its distance from the origin.
     let x0 = &pts[0];
-    let x0_sq = squared_norm(x0);
     let mut matrix = vec![vec![0.0; dim]; dim];
     let mut rhs = vec![0.0; dim];
 
     for row in 0..dim {
-        let point = &pts[row + 1];
-        rhs[row] = squared_norm(point) - x0_sq;
-        for col in 0..dim {
-            matrix[row][col] = 2.0 * (point[col] - x0[col]);
+        let translated: Vec<f64> = pts[row + 1].iter().zip(x0).map(|(a, b)| a - b).collect();
+        rhs[row] = squared_norm(&translated);
+        for (col, value) in translated.iter().enumerate() {
+            matrix[row][col] = 2.0 * value;
         }
     }
 
-    let center = match solve_square(&matrix, &rhs) {
+    let relative_center = match solve_square(&matrix, &rhs) {
         Ok(center) => center,
         Err(GeometryError::SingularMatrix) => {
             return Ok((vec![f64::NAN; dim], f64::NAN));
         }
         Err(err) => return Err(err),
     };
 
-    let radius = fast_norm(
-        &center
-            .iter()
-            .zip(&pts[0])
-            .map(|(x, y)| x - y)
-            .collect::<Vec<_>>(),
-    );
+    let radius = fast_norm(&relative_center);
+    let center: Vec<f64> = relative_center.iter().zip(x0).map(|(c, o)| c + o).collect();
     Ok((center, radius))
 }
 
@@ -424,10 +421,12 @@ pub fn simplex_volume_in_embedding(vertices: &[Vec<f64>]) -> Result<f64, Geometr
         ));
     }
 
-    if vertices[0].len() == 2 && vertices.len() != 3 {
-        return Err(GeometryError::InvalidDimensions(
-            "Expected three 2D vertices".to_string(),
-        ));
+    if vertices.len() == 2 {
+        let length_sq = squared_distance(&vertices[0], &vertices[1]);
+        if length_sq == 0.0 {
+            return Err(GeometryError::DegenerateSimplex);
+        }
+        return Ok(length_sq.sqrt());
     }
 
     if vertices.len() == 3 {
@@ -465,3 +464,48 @@ pub fn simplex_volume_in_embedding(vertices: &[Vec<f64>]) -> Result<f64, Geometr
     }
     Ok(vol_square.sqrt())
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn simplex_volume_in_embedding_returns_segment_length() {
+        let length =
+            simplex_volume_in_embedding(&[vec![0.0, 0.0, 0.0], vec![3.0, 4.0, 0.0]]).unwrap();
+        assert!((length - 5.0).abs() < 1e-12);
+    }
+
+    #[test]
+    fn simplex_volume_in_embedding_rejects_identical_endpoints() {
+        let err = simplex_volume_in_embedding(&[vec![1.0, 2.0], vec![1.0, 2.0]]).unwrap_err();
+        assert!(matches!(err, GeometryError::DegenerateSimplex));
+    }
+
+    #[test]
+    fn circumsphere_supports_one_dimensional_segment() {
+        let (center, radius) = circumsphere(&[vec![0.0], vec![2.0]]).unwrap();
+        assert_eq!(center, vec![1.0]);
+        assert!((radius - 1.0).abs() < 1e-12);
+    }
+
+    #[test]
+    fn circumsphere_is_accurate_for_small_simplices_far_from_origin() {
+        let (center, radius) = circumsphere(&[vec![0.5], vec![0.5 + 1e-6]]).unwrap();
+        assert!((center[0] - (0.5 + 5e-7)).abs() < 1e-15);
+        assert!((radius - 5e-7).abs() < 1e-15);
+
+        let (center, radius) = circumsphere(&[
+            vec![100.0, 100.0, 100.0, 100.0],
+            vec![100.0 + 1e-5, 100.0, 100.0, 100.0],
+            vec![100.0, 100.0 + 1e-5, 100.0, 100.0],
+            vec![100.0, 100.0, 100.0 + 1e-5, 100.0],
+            vec![100.0, 100.0, 100.0, 100.0 + 1e-5],
+        ])
+        .unwrap();
+        for coord in &center {
+            assert!((coord - (100.0 + 5e-6)).abs() < 1e-12);
+        }
+        assert!((radius - 1e-5).abs() < 1e-12);
+    }
+}