diff --git a/benches/benchmark_utilities.rs b/benches/benchmark_utilities.rs index 23029a06..a458de0f 100644 --- a/benches/benchmark_utilities.rs +++ b/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, (), (), (), LastUsedVertexHintGenerator>; type HierarchyTriangulation = DelaunayTriangulation, (), (), (), HierarchyHintGenerator>; +type JumpAndWalkTriangulation = + DelaunayTriangulation, (), (), (), 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,6 +63,7 @@ where pub enum CreationMethod { BulkLoad, Incremental, + JumpAndWalk, } impl Display for CreationMethod { @@ -68,6 +71,7 @@ impl Display for CreationMethod { f.write_str(match self { CreationMethod::BulkLoad => "bulk loading", CreationMethod::Incremental => "incremental loading", + CreationMethod::JumpAndWalk => "jump and walk", }) } } @@ -82,6 +86,7 @@ pub enum SampleSize { pub enum HintGeneratorType { LastUsedVertex, Hierarchy, + JumpAndWalk, } impl Display for HintGeneratorType { @@ -89,6 +94,7 @@ impl Display for HintGeneratorType { 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::>(); 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(); + } + } } } } diff --git a/benches/bulk_load_benchmark.rs b/benches/bulk_load_benchmark.rs index 7f936291..01402fa7 100644 --- a/benches/bulk_load_benchmark.rs +++ b/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, diff --git a/benches/bulk_load_vs_incremental.rs b/benches/bulk_load_vs_incremental.rs index e121bbae..cdc4704e 100644 --- a/benches/bulk_load_vs_incremental.rs +++ b/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, diff --git a/benches/locate_benchmark.rs b/benches/locate_benchmark.rs index 59b71b2f..f23b182a 100644 --- a/benches/locate_benchmark.rs +++ b/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::( + &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); diff --git a/src/delaunay_core/hint_generator.rs b/src/delaunay_core/hint_generator.rs index 870b40de..ff84d5b9 100644 --- a/src/delaunay_core/hint_generator.rs +++ b/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: 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) -> FixedVertexHandle; + fn get_hint, V: HasPosition>( + &self, + position: Point2, + 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 HintGenerator for LastUsedVertexHintGenerator { - fn get_hint(&self, _: Point2) -> FixedVertexHandle { + fn get_hint, V: HasPosition>( + &self, + _: Point2, + _triangulation: &TR, + ) -> FixedVertexHandle { FixedVertexHandle::new(self.index.load(Ordering::Relaxed)) } @@ -172,7 +183,11 @@ impl Default impl HintGenerator for HierarchyHintGeneratorWithBranchFactor { - fn get_hint(&self, position: Point2) -> FixedVertexHandle { + fn get_hint, V: HasPosition>( + &self, + position: Point2, + _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 HintGenerator } } +/// 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 HintGenerator for JumpAndWalkHintGenerator { + fn get_hint, V: HasPosition>( + &self, + position: Point2, + 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) {} + fn notify_vertex_removed( + &mut self, + _swapped_in_point: Option>, + _vertex: FixedVertexHandle, + _vertex_position: Point2, + ) { + } + + fn initialize_from_triangulation(_triangulation: &TR) -> Self + where + TR: Triangulation, + V: HasPosition, + { + 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}"); + } + } } diff --git a/src/delaunay_core/mod.rs b/src/delaunay_core/mod.rs index 2509653d..bf09a8ca 100644 --- a/src/delaunay_core/mod.rs +++ b/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}; diff --git a/src/delaunay_core/triangulation_ext.rs b/src/delaunay_core/triangulation_ext.rs index d4d588fb..b0a13c15 100644 --- a/src/delaunay_core/triangulation_ext.rs +++ b/src/delaunay_core/triangulation_ext.rs @@ -217,7 +217,7 @@ pub trait TriangulationExt: Triangulation { point: Point2<::Scalar>, hint: Option, ) -> 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) } diff --git a/src/delaunay_triangulation.rs b/src/delaunay_triangulation.rs index 7b4f6343..967a1368 100644 --- a/src/delaunay_triangulation.rs +++ b/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); diff --git a/src/lib.rs b/src/lib.rs index 472a081a..f8fe33ee 100644 --- a/src/lib.rs +++ b/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}; diff --git a/src/triangulation.rs b/src/triangulation.rs index be7096de..7644942a 100644 --- a/src/triangulation.rs +++ b/src/triangulation.rs @@ -403,7 +403,7 @@ pub trait Triangulation: Default { &self, point: Point2<::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)) }