Jump-and-Walk

benches/benchmark_utilities.rs

@@ -3,14 +3,16 @@ use core::fmt::{self, Display, Formatter};
 use criterion::{measurement::WallTime, BenchmarkGroup, BenchmarkId, Throughput};
 use rand::{distributions::uniform::SampleUniform, Rng, SeedableRng};
 use spade::{
-    DelaunayTriangulation, HierarchyHintGenerator, LastUsedVertexHintGenerator, Point2, SpadeNum,
-    Triangulation,
+    DelaunayTriangulation, HierarchyHintGenerator, JumpAndWalkHintGenerator,
+    LastUsedVertexHintGenerator, Point2, SpadeNum, Triangulation,
 };
 
 type LastUsedTriangulation =
     DelaunayTriangulation<Point2<f64>, (), (), (), LastUsedVertexHintGenerator>;
 type HierarchyTriangulation =
     DelaunayTriangulation<Point2<f64>, (), (), (), HierarchyHintGenerator<f64>>;
+type JumpAndWalkTriangulation =
+    DelaunayTriangulation<Point2<f64>, (), (), (), JumpAndWalkHintGenerator>;
 
 pub const SEED: &[u8; 32] = b"\xfb\xdc\x4e\xa0\x30\xde\x82\xba\x69\x97\x3c\x52\x49\x4d\x00\xca
 \x5c\x21\xa3\x8d\x5c\xf2\x34\x4e\x58\x7d\x80\x16\x66\x23\x30";
@@ -61,13 +63,15 @@ where
 pub enum CreationMethod {
     BulkLoad,
     Incremental,
+    JumpAndWalk,
 }
 
 impl Display for CreationMethod {
     fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
         f.write_str(match self {
             CreationMethod::BulkLoad => "bulk loading",
             CreationMethod::Incremental => "incremental loading",
+            CreationMethod::JumpAndWalk => "jump and walk",
         })
     }
 }
@@ -82,13 +86,15 @@ pub enum SampleSize {
 pub enum HintGeneratorType {
     LastUsedVertex,
     Hierarchy,
+    JumpAndWalk,
 }
 
 impl Display for HintGeneratorType {
     fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
         f.write_str(match self {
             HintGeneratorType::LastUsedVertex => "last used vertex heuristic",
             HintGeneratorType::Hierarchy => "hierarchy lookup",
+            HintGeneratorType::JumpAndWalk => "jump and walk",
         })
     }
 }
@@ -127,6 +133,7 @@ impl CreationBenchConfig {
                 let data = self.get_distribution().take(*size).collect::<Vec<_>>();
                 match self.creation_method {
                     CreationMethod::BulkLoad => b.iter(|| self.bulk_load(data.clone())),
+                    CreationMethod::JumpAndWalk => b.iter(|| self.bulk_load(data.clone())),
                     CreationMethod::Incremental => b.iter(|| self.sequential_load(data.clone())),
                 }
             });
@@ -165,6 +172,9 @@ impl CreationBenchConfig {
             HintGeneratorType::Hierarchy => {
                 HierarchyTriangulation::bulk_load(data).unwrap();
             }
+            HintGeneratorType::JumpAndWalk => {
+                HierarchyTriangulation::bulk_load(data).unwrap();
+            }
         };
     }
 
@@ -182,6 +192,12 @@ impl CreationBenchConfig {
                     triangulation.insert(point).unwrap();
                 }
             }
+            HintGeneratorType::JumpAndWalk => {
+                let mut triangulation = JumpAndWalkTriangulation::new();
+                for point in data {
+                    triangulation.insert(point).unwrap();
+                }
+            }
         }
     }
 }

benches/bulk_load_benchmark.rs

@@ -13,6 +13,8 @@ pub fn bulk_load_benchmark(c: &mut Criterion) {
         (LastUsedVertex, Uniform),
         (Hierarchy, Uniform),
         (Hierarchy, RandomWalk),
+        (JumpAndWalk, Uniform),
+        (JumpAndWalk, RandomWalk),
     ] {
         let config = CreationBenchConfig {
             creation_method: CreationMethod::BulkLoad,

benches/bulk_load_vs_incremental.rs

@@ -6,7 +6,7 @@ pub fn bulk_load_vs_incremental_benchmark(c: &mut Criterion) {
     {
         use CreationMethod::*;
         let mut group = c.benchmark_group("bulk vs incremental loading");
-        for creation_method in [BulkLoad, Incremental] {
+        for creation_method in [BulkLoad, Incremental, JumpAndWalk] {
             let config = CreationBenchConfig {
                 creation_method,
                 sample_size: SampleSize::SmallSampleSet,

benches/locate_benchmark.rs

@@ -3,8 +3,8 @@ use std::time::Duration;
 use criterion::{measurement::WallTime, BenchmarkGroup, Criterion};
 use rand::distributions::uniform::SampleUniform;
 use spade::{
-    DelaunayTriangulation, HierarchyHintGeneratorWithBranchFactor, HintGenerator, Point2, SpadeNum,
-    Triangulation,
+    DelaunayTriangulation, HierarchyHintGeneratorWithBranchFactor, HintGenerator,
+    JumpAndWalkHintGenerator, Point2, SpadeNum, Triangulation,
 };
 
 use crate::benchmark_utilities::{uniform_distribution, uniform_f64, SEED2};
@@ -56,6 +56,13 @@ pub fn locate_benchmark(c: &mut Criterion) {
         );
     }
 
+    single_locate_benchmark::<JumpAndWalkHintGenerator, _, _>(
+        &mut group,
+        "locate (jump and walk)".to_string(),
+        RANGE,
+        uniform_f64(),
+    );
+
     single_hierarchy::<2>(&mut group);
     single_hierarchy::<3>(&mut group);
     single_hierarchy::<4>(&mut group);

src/delaunay_core/hint_generator.rs

@@ -1,4 +1,7 @@
-use core::sync::atomic::{AtomicUsize, Ordering};
+use core::{
+    f64,
+    sync::atomic::{AtomicUsize, Ordering},
+};
 
 use crate::{
     DelaunayTriangulation, HasPosition, Point2, SpadeNum, Triangulation, TriangulationExt,
@@ -31,7 +34,11 @@ pub trait HintGenerator<S: SpadeNum>: Default {
     /// Returns a vertex handle that should be close to a given position.
     ///
     /// The returned vertex handle may be invalid.
-    fn get_hint(&self, position: Point2<S>) -> FixedVertexHandle;
+    fn get_hint<TR: Triangulation<Vertex = V>, V: HasPosition<Scalar = S>>(
+        &self,
+        position: Point2<S>,
+        triangulation: &TR,
+    ) -> FixedVertexHandle;
 
     /// Notifies the hint generator that an element was looked up
     fn notify_vertex_lookup(&self, vertex: FixedVertexHandle);
@@ -91,7 +98,11 @@ impl Clone for LastUsedVertexHintGenerator {
 }
 
 impl<S: SpadeNum> HintGenerator<S> for LastUsedVertexHintGenerator {
-    fn get_hint(&self, _: Point2<S>) -> FixedVertexHandle {
+    fn get_hint<TR: Triangulation<Vertex = V>, V: HasPosition<Scalar = S>>(
+        &self,
+        _: Point2<S>,
+        _triangulation: &TR,
+    ) -> FixedVertexHandle {
         FixedVertexHandle::new(self.index.load(Ordering::Relaxed))
     }
 
@@ -172,7 +183,11 @@ impl<S: SpadeNum, const BRANCH_FACTOR: u32> Default
 impl<S: SpadeNum, const BRANCH_FACTOR: u32> HintGenerator<S>
     for HierarchyHintGeneratorWithBranchFactor<S, BRANCH_FACTOR>
 {
-    fn get_hint(&self, position: Point2<S>) -> FixedVertexHandle {
+    fn get_hint<TR: Triangulation<Vertex = V>, V: HasPosition<Scalar = S>>(
+        &self,
+        position: Point2<S>,
+        _triangulation: &TR,
+    ) -> FixedVertexHandle {
         let mut nearest = FixedVertexHandle::new(0);
         for layer in self.hierarchy.iter().rev().skip(1) {
             nearest = layer.walk_to_nearest_neighbor(nearest, position).fix();
@@ -283,6 +298,80 @@ impl<S: SpadeNum, const BRANCH_FACTOR: u32> HintGenerator<S>
     }
 }
 
+/// Using the Jump-and-Walk algorithm to find points.
+#[derive(Default, Debug)]
+#[cfg_attr(
+    feature = "serde",
+    derive(Serialize, Deserialize),
+    serde(crate = "serde")
+)]
+pub struct JumpAndWalkHintGenerator;
+
+impl<S: SpadeNum> HintGenerator<S> for JumpAndWalkHintGenerator {
+    fn get_hint<TR: Triangulation<Vertex = V>, V: HasPosition<Scalar = S>>(
+        &self,
+        position: Point2<S>,
+        triangulation: &TR,
+    ) -> FixedVertexHandle {
+        let mut vertices = triangulation.s().vertices();
+        let Some(first) = vertices.next() else {
+            panic!("no vertices");
+        };
+        let len = vertices.len();
+        let step = (len / approx_cube_root(len)).max(1);
+
+        // Jump
+        let mut best_dist_2 = first.position().distance_2(position);
+        let mut best_vert = first.fix();
+        for vert in vertices.step_by(step) {
+            let dist_2 = vert.position().distance_2(position);
+            if dist_2 < best_dist_2 {
+                best_dist_2 = dist_2;
+                best_vert = vert.fix();
+            }
+        }
+
+        best_vert
+    }
+
+    fn notify_vertex_lookup(&self, _vertex: FixedVertexHandle) {}
+    fn notify_vertex_inserted(&mut self, _vertex: FixedVertexHandle, _vertex_position: Point2<S>) {}
+    fn notify_vertex_removed(
+        &mut self,
+        _swapped_in_point: Option<Point2<S>>,
+        _vertex: FixedVertexHandle,
+        _vertex_position: Point2<S>,
+    ) {
+    }
+
+    fn initialize_from_triangulation<TR, V>(_triangulation: &TR) -> Self
+    where
+        TR: Triangulation<Vertex = V>,
+        V: HasPosition<Scalar = S>,
+    {
+        Self
+    }
+}
+
+fn approx_cube_root(value: usize) -> usize {
+    const LEADING_ZEROS_TO_GUESS: [u32; 64] = [
+        1905390, 1512309, 1200320, 952695, 756155, 600160, 476348, 378078, 300080, 238174, 189039,
+        150040, 119087, 94520, 75020, 59544, 47260, 37510, 29772, 23630, 18755, 14886, 11815, 9378,
+        7443, 5908, 4689, 3722, 2954, 2345, 1861, 1477, 1173, 931, 739, 587, 466, 370, 294, 233,
+        185, 147, 117, 93, 74, 59, 47, 37, 30, 24, 19, 15, 12, 10, 8, 6, 5, 4, 3, 3, 2, 2, 2, 1,
+    ];
+    if value < 2 {
+        return value;
+    }
+    let index = value.leading_zeros() as usize;
+    let mut guess = LEADING_ZEROS_TO_GUESS[index] as usize;
+
+    // One iteration of Newton's method.
+    guess = (value / (guess * guess) + guess * 2) / 3;
+
+    guess
+}
+
 #[cfg(test)]
 mod test {
     use rand::{prelude::SliceRandom, RngCore, SeedableRng};
@@ -292,6 +381,8 @@ mod test {
         Triangulation, TriangulationExt,
     };
 
+    use super::approx_cube_root;
+
     use alloc::vec::Vec;
 
     const BRANCH_FACTOR: u32 = 3;
@@ -372,4 +463,17 @@ mod test {
         }
         Ok(())
     }
+
+    #[test]
+    fn test_approx_cube_root() {
+        for i in 0usize.. {
+            let Some(i3) = (i * i).checked_mul(i) else {
+                break;
+            };
+            let sut = approx_cube_root(i3);
+            // Tolerance of 0.015% of the input.
+            let tolerance = i3 / 6859;
+            assert!(sut.abs_diff(i) <= tolerance, "{sut} ~= {i}");
+        }
+    }
 }

src/delaunay_core/mod.rs

@@ -22,7 +22,7 @@ pub use triangulation_ext::{RemovalResult, TriangulationExt};
 pub use dcel::Dcel;
 pub use hint_generator::{
     HierarchyHintGenerator, HierarchyHintGeneratorWithBranchFactor, HintGenerator,
-    LastUsedVertexHintGenerator,
+    JumpAndWalkHintGenerator, LastUsedVertexHintGenerator,
 };
 
 pub use refinement::{AngleLimit, RefinementParameters, RefinementResult};

src/delaunay_core/triangulation_ext.rs

@@ -217,7 +217,7 @@ pub trait TriangulationExt: Triangulation {
         point: Point2<<Self::Vertex as HasPosition>::Scalar>,
         hint: Option<FixedVertexHandle>,
     ) -> PositionInTriangulation {
-        let start = hint.unwrap_or_else(|| self.hint_generator().get_hint(point));
+        let start = hint.unwrap_or_else(|| self.hint_generator().get_hint(point, self));
         self.locate_with_hint_fixed_core(point, start)
     }
 

src/delaunay_triangulation.rs

@@ -297,7 +297,7 @@ where
             return None;
         }
 
-        let hint = self.hint_generator().get_hint(position);
+        let hint = self.hint_generator().get_hint(position, self);
         let hint = self.validate_vertex_handle(hint);
 
         let vertex = self.walk_to_nearest_neighbor(hint, position);

src/lib.rs

@@ -52,7 +52,7 @@ pub use crate::delaunay_core::math::{
 
 pub use delaunay_core::{
     AngleLimit, HierarchyHintGenerator, HierarchyHintGeneratorWithBranchFactor, HintGenerator,
-    LastUsedVertexHintGenerator, RefinementParameters, RefinementResult,
+    JumpAndWalkHintGenerator, LastUsedVertexHintGenerator, RefinementParameters, RefinementResult,
 };
 
 pub use crate::delaunay_core::interpolation::{Barycentric, NaturalNeighbor};

src/triangulation.rs

@@ -403,7 +403,7 @@ pub trait Triangulation: Default {
         &self,
         point: Point2<<Self::Vertex as HasPosition>::Scalar>,
     ) -> PositionInTriangulation {
-        let hint = self.hint_generator().get_hint(point);
+        let hint = self.hint_generator().get_hint(point, self);
         self.locate_with_hint_option_core(point, Some(hint))
     }