diff --git a/benches/allocation_hot_paths.rs b/benches/allocation_hot_paths.rs index 75a4dd4f..4fc117d1 100644 --- a/benches/allocation_hot_paths.rs +++ b/benches/allocation_hot_paths.rs @@ -24,9 +24,9 @@ mod allocation_contracts { use delaunay::prelude::construction::{ ConstructionOptions, DelaunayTriangulation, RetryPolicy, Vertex, vertex, }; - use delaunay::prelude::generators::generate_random_points_seeded; + use delaunay::prelude::generators::generate_random_points_in_range_seeded; use delaunay::prelude::geometry::{ - AdaptiveKernel, Coordinate, FastKernel, Point, simplex_volume, + AdaptiveKernel, Coordinate, CoordinateRange, FastKernel, Point, simplex_volume, }; use delaunay::prelude::query::measure_with_result; use delaunay::prelude::tds::{SimplexKey, TdsError, VertexKey, facet_key_from_vertices}; @@ -89,11 +89,16 @@ mod allocation_contracts { attempts } + fn benchmark_bounds() -> CoordinateRange { + bench_result( + CoordinateRange::try_new(-100.0_f64, 100.0), + "allocation benchmark bounds must be valid", + ) + } + fn canary_vertices(count: usize, seed: u64) -> Vec> { - let points = bench_result( - generate_random_points_seeded::(count, (-100.0, 100.0), seed), - format!("failed to generate {D}D allocation benchmark points"), - ); + let points = + generate_random_points_in_range_seeded::(count, benchmark_bounds(), seed); points.into_iter().map(|point| vertex!(point)).collect() } diff --git a/benches/boundary_uuid_iter.rs b/benches/boundary_uuid_iter.rs index 5a5875a9..c676c19a 100644 --- a/benches/boundary_uuid_iter.rs +++ b/benches/boundary_uuid_iter.rs @@ -8,7 +8,8 @@ use criterion::{BenchmarkId, Criterion, Throughput, criterion_group, criterion_main}; use delaunay::prelude::construction::{DelaunayTriangulation, Vertex, vertex}; -use delaunay::prelude::generators::generate_random_points_seeded; +use delaunay::prelude::generators::generate_random_points_in_range_seeded; +use delaunay::prelude::geometry::CoordinateRange; use delaunay::prelude::query::BoundaryAnalysis; use uuid::Uuid; @@ -20,19 +21,22 @@ use std::hint::black_box; pub mod bench_utils; use bench_utils::{bench_option, bench_result}; -const BOUNDS: (f64, f64) = (-100.0, 100.0); const BOUNDARY_COUNTS_3D: &[usize] = &[20, 40, 60, 80]; +fn benchmark_bounds() -> CoordinateRange { + bench_result( + CoordinateRange::try_new(-100.0_f64, 100.0), + "boundary benchmark bounds must be valid", + ) +} + fn boundary_triangulation_3d( requested_vertices: usize, ) -> DelaunayTriangulation, (), (), 3> { - let points = bench_result( - generate_random_points_seeded::( - requested_vertices, - BOUNDS, - 0xB0DA_FACE_0000_0000 ^ requested_vertices as u64, - ), - "failed to generate 3D boundary benchmark points", + let points = generate_random_points_in_range_seeded::( + requested_vertices, + benchmark_bounds(), + 0xB0DA_FACE_0000_0000 ^ requested_vertices as u64, ); let vertices = Vertex::from_points(&points); bench_result( diff --git a/benches/ci_performance_suite.rs b/benches/ci_performance_suite.rs index 345c61b4..be9e74a4 100644 --- a/benches/ci_performance_suite.rs +++ b/benches/ci_performance_suite.rs @@ -47,8 +47,8 @@ use delaunay::prelude::construction::{ ConstructionOptions, DelaunayTriangulation, RetryPolicy, Vertex, }; use delaunay::prelude::flips::{FacetHandle, RidgeHandle, SimplexKey}; -use delaunay::prelude::generators::generate_random_points_seeded; -use delaunay::prelude::geometry::{AdaptiveKernel, Coordinate, Point}; +use delaunay::prelude::generators::generate_random_points_in_range_seeded; +use delaunay::prelude::geometry::{AdaptiveKernel, Coordinate, CoordinateRange, Point}; use delaunay::prelude::query::ConvexHull; use delaunay::vertex; use std::{env, hint::black_box, num::NonZeroUsize, sync::Once}; @@ -321,7 +321,7 @@ fn print_manifest_once() { fn prepare_data( dim_seed: u64, count: usize, - bounds: (f64, f64), + bounds: CoordinateRange, attempts: NonZeroUsize, ) -> (u64, Vec>, Vec>) { // Fast path: use the pre-computed seed (single verification construction) @@ -340,7 +340,7 @@ fn prepare_data( find_seed_vertices::(base_seed, count, bounds, search_limit, attempts), format_args!( "No stable benchmark seed found for {D}D/{count}: \ - start_seed={base_seed}; search_limit={search_limit}; bounds={bounds:?}" + start_seed={base_seed}; search_limit={search_limit}; bounds={bounds}" ), ) } @@ -366,7 +366,10 @@ fn warn_known_seed_failed(seed: u64, count: usize, dataset: Data } fn prepare_dt(dim_seed: u64, count: usize) -> BenchTriangulation { - let bounds = (-100.0, 100.0); + let bounds = bench_result( + CoordinateRange::try_new(-100.0_f64, 100.0), + "well-conditioned benchmark bounds must be valid", + ); let attempts = retry_attempts(6); let (seed, _, vertices) = prepare_data::(dim_seed, count, bounds, attempts); let options = ConstructionOptions::default().with_retry_policy(RetryPolicy::Shuffled { @@ -404,9 +407,13 @@ fn prepare_inserts( seed ^= 0xA5A5_A5A5; } let points = match dataset { - Dataset::WellConditioned => bench_result( - generate_random_points_seeded::(count, (-50.0, 50.0), seed), - format!("insert point generation failed for {D}D"), + Dataset::WellConditioned => generate_random_points_in_range_seeded::( + count, + bench_result( + CoordinateRange::try_new(-50.0_f64, 50.0), + "insert benchmark bounds must be valid", + ), + seed, ), Dataset::Adversarial => generate_adv_points::(count, seed), }; @@ -416,16 +423,14 @@ fn prepare_inserts( fn find_seed_vertices( start_seed: u64, count: usize, - bounds: (f64, f64), + bounds: CoordinateRange, limit: usize, attempts: NonZeroUsize, ) -> SeedSearchResult { for offset in 0..limit { let candidate_seed = start_seed.wrapping_add(offset as u64); - let points = bench_result( - generate_random_points_seeded::(count, bounds, candidate_seed), - format!("generate_random_points_seeded failed for {D}D"), - ); + let points = + generate_random_points_in_range_seeded::(count, bounds, candidate_seed); let vertices = points.iter().map(|p| vertex!(*p)).collect::>(); let options = ConstructionOptions::default().with_retry_policy(RetryPolicy::Shuffled { @@ -496,9 +501,13 @@ fn prepare_adv_data( } fn generate_adv_points(count: usize, seed: u64) -> Vec> { - let base_points = bench_result( - generate_random_points_seeded::(count, (-1.0, 1.0), seed), - format!("generate_random_points_seeded failed for adversarial {D}D"), + let base_points = generate_random_points_in_range_seeded::( + count, + bench_result( + CoordinateRange::try_new(-1.0_f64, 1.0), + "adversarial benchmark bounds must be valid", + ), + seed, ); base_points @@ -900,7 +909,10 @@ macro_rules! benchmark_tds_new_dimension { // We avoid `std::process::exit` here so that destructors run and Criterion // can clean up state on both success and failure. if discover_seeds_enabled() { - let bounds = (-100.0, 100.0); + let bounds = bench_result( + CoordinateRange::try_new(-100.0_f64, 100.0), + "well-conditioned benchmark bounds must be valid", + ); let filters = criterion_filters(); let bench_id = format!("tds_new_{}d/tds_new/{count}", stringify!($dim)); @@ -945,7 +957,10 @@ macro_rules! benchmark_tds_new_dimension { format!("tds_new_{}d/tds_new_adversarial/{count}", stringify!($dim)); if benchmark_selected(&filters, &bench_id) { - let bounds = (-100.0, 100.0); + let bounds = bench_result( + CoordinateRange::try_new(-100.0_f64, 100.0), + "well-conditioned benchmark bounds must be valid", + ); let attempts = retry_attempts(6); let (seed, _, vertices) = prepare_data::<$dim>($seed, count, bounds, attempts); let options = ConstructionOptions::default().with_retry_policy( @@ -985,7 +1000,10 @@ macro_rules! benchmark_tds_new_dimension { group.bench_with_input(BenchmarkId::new("tds_new", count), &count, |b, &count| { // Reduce variance: pre-generate deterministic inputs outside the measured loop, // then benchmark only triangulation construction. - let bounds = (-100.0, 100.0); + let bounds = bench_result( + CoordinateRange::try_new(-100.0_f64, 100.0), + "well-conditioned benchmark bounds must be valid", + ); let attempts = retry_attempts(6); let (seed, points, vertices) = prepare_data::<$dim>($seed, count, bounds, attempts); diff --git a/benches/circumsphere_containment.rs b/benches/circumsphere_containment.rs index 6d81aa75..27aa51de 100644 --- a/benches/circumsphere_containment.rs +++ b/benches/circumsphere_containment.rs @@ -14,7 +14,8 @@ //! - Numerical consistency validation between all three algorithms use criterion::{Criterion, criterion_group, criterion_main}; -use delaunay::prelude::generators::generate_random_points_seeded; +use delaunay::prelude::generators::generate_random_points_in_range_seeded; +use delaunay::prelude::geometry::CoordinateRange; use delaunay::prelude::query::*; use std::hint::black_box; @@ -23,6 +24,10 @@ use std::hint::black_box; pub mod bench_utils; use bench_utils::{abort_benchmark, bench_option, bench_result}; +fn coordinate_range(min: f64, max: f64, context: &'static str) -> CoordinateRange { + bench_result(CoordinateRange::try_new(min, max), context) +} + /// Generate a standard D-dimensional simplex (D+1 vertices) /// /// Creates a simplex with vertices at: @@ -46,17 +51,19 @@ fn standard_simplex() -> Vec> { /// Generate a random 3D simplex (tetrahedron) for benchmarking using seeded generation fn generate_random_simplex_3d(seed: u64) -> Vec> { - bench_result( - generate_random_points_seeded(4, (-10.0, 10.0), seed), - "failed to generate random simplex points", + generate_random_points_in_range_seeded( + 4, + coordinate_range(-10.0, 10.0, "random simplex bounds must be valid"), + seed, ) } /// Generate a random 3D test point using seeded generation fn generate_random_test_point_3d(seed: u64) -> Point { - let points = bench_result( - generate_random_points_seeded(1, (-5.0, 5.0), seed), - "failed to generate random test point", + let points = generate_random_points_in_range_seeded( + 1, + coordinate_range(-5.0, 5.0, "random test point bounds must be valid"), + seed, ); bench_option(points.into_iter().next(), "expected exactly one test point") } @@ -67,9 +74,10 @@ fn benchmark_random_queries(c: &mut Criterion) { let simplex_points = generate_random_simplex_3d(42); // Generate many test points using seeded generation for reproducible results - let test_points = bench_result( - generate_random_points_seeded(1000, (-5.0, 5.0), 123), - "failed to generate random test points", + let test_points = generate_random_points_in_range_seeded( + 1000, + coordinate_range(-5.0, 5.0, "random query bounds must be valid"), + 123, ); c.bench_function("random/insphere_1000_queries", |b| { diff --git a/benches/cold_path_predicates.rs b/benches/cold_path_predicates.rs index 597c8718..9fe152d9 100644 --- a/benches/cold_path_predicates.rs +++ b/benches/cold_path_predicates.rs @@ -35,7 +35,8 @@ //! ``` use criterion::{BenchmarkId, Criterion, Throughput, criterion_group, criterion_main}; -use delaunay::prelude::generators::generate_random_points_seeded; +use delaunay::prelude::generators::generate_random_points_in_range_seeded; +use delaunay::prelude::geometry::CoordinateRange; use delaunay::prelude::query::*; use std::hint::black_box; @@ -44,6 +45,10 @@ use std::hint::black_box; pub mod bench_utils; use bench_utils::{abort_benchmark, bench_result}; +fn coordinate_range(min: f64, max: f64, context: &'static str) -> CoordinateRange { + bench_result(CoordinateRange::try_new(min, max), context) +} + /// Deterministic seed for query-point generation in the hot path. const HOT_SEED: u64 = 0xC01D_BEEF_0000_CAFE_u64; /// Deterministic seed for query-point generation in the near-boundary group. @@ -72,9 +77,10 @@ fn standard_simplex() -> Vec> { /// Uses the range `[-10, 10]` against a unit simplex so that the Shewchuk /// errbound comfortably resolves the sign in Stage 1. fn hot_queries() -> Vec> { - bench_result( - generate_random_points_seeded(HOT_QUERIES, (-10.0, 10.0), HOT_SEED), - "failed to generate hot-path query points", + generate_random_points_in_range_seeded( + HOT_QUERIES, + coordinate_range(-10.0, 10.0, "hot-path query bounds must be valid"), + HOT_SEED, ) } @@ -86,9 +92,10 @@ fn near_boundary_queries() -> Vec> { // Centered near the circumsphere radius of the standard simplex (~0.5 for // the D = 3 unit case); the exact value is unimportant — we just want a // high rate of errbound-ambiguous inputs. - bench_result( - generate_random_points_seeded(NEAR_BOUNDARY_QUERIES, (0.40, 0.60), NEAR_BOUNDARY_SEED), - "failed to generate near-boundary query points", + generate_random_points_in_range_seeded( + NEAR_BOUNDARY_QUERIES, + coordinate_range(0.40, 0.60, "near-boundary query bounds must be valid"), + NEAR_BOUNDARY_SEED, ) } diff --git a/benches/profiling_suite.rs b/benches/profiling_suite.rs index 3adae7ad..7059fd69 100644 --- a/benches/profiling_suite.rs +++ b/benches/profiling_suite.rs @@ -79,8 +79,8 @@ use delaunay::prelude::collections::SmallBuffer; use delaunay::prelude::construction::{ ConstructionOptions, DelaunayTriangulation, DelaunayTriangulationBuilder, RetryPolicy, Vertex, }; -use delaunay::prelude::generators::generate_random_points_seeded; -use delaunay::prelude::geometry::{AdaptiveKernel, Coordinate, Point}; +use delaunay::prelude::generators::generate_random_points_in_range_seeded; +use delaunay::prelude::geometry::{AdaptiveKernel, Coordinate, CoordinateRange, Point}; use delaunay::prelude::query::*; use delaunay::vertex; use sysinfo::{ProcessRefreshKind, ProcessesToUpdate, RefreshKind, System, get_current_pid}; @@ -97,6 +97,18 @@ fn retry_attempts(value: usize) -> NonZeroUsize { attempts } +fn coordinate_range(min: f64, max: f64, context: &'static str) -> CoordinateRange { + bench_result(CoordinateRange::try_new(min, max), context) +} + +fn wide_bounds() -> CoordinateRange { + coordinate_range(-100.0, 100.0, "wide benchmark bounds must be valid") +} + +fn adversarial_bounds() -> CoordinateRange { + coordinate_range(-1.0, 1.0, "adversarial benchmark bounds must be valid") +} + #[cfg(feature = "bench-logging")] fn init_tracing() { static INIT: Once = Once::new(); @@ -292,10 +304,7 @@ fn vertices_from_points(points: Vec>) -> Vec(n_points: usize, seed: u64) -> MemoryInfo { let mem_before = memory_usage_kib(); - let points = bench_result( - generate_random_points_seeded::(n_points, (-100.0, 100.0), seed), - "failed to generate points", - ); + let points = generate_random_points_in_range_seeded::(n_points, wide_bounds(), seed); let vertices = vertices_from_points(points); let mem_before_tds = memory_usage_kib(); @@ -349,17 +358,13 @@ fn gen_points( seed: u64, ) -> Vec> { match distribution { - PointDistribution::Random => bench_result( - generate_random_points_seeded(count, (-100.0, 100.0), seed), - "random point generation failed", - ), - PointDistribution::Adversarial => bench_result( - generate_random_points_seeded::( - count, - (-1.0, 1.0), - seed ^ 0xA5A5_A5A5_A5A5_A5A5, - ), - "adversarial base point generation failed", + PointDistribution::Random => { + generate_random_points_in_range_seeded(count, wide_bounds(), seed) + } + PointDistribution::Adversarial => generate_random_points_in_range_seeded::( + count, + adversarial_bounds(), + seed ^ 0xA5A5_A5A5_A5A5_A5A5, ) .iter() .enumerate() @@ -404,10 +409,8 @@ fn bench_construction(c: &mut Criterion, dimension_name: &str, n group.bench_function("construct", |b| { b.iter_batched( || { - let points = bench_result( - generate_random_points_seeded::(n_points, (-100.0, 100.0), seed), - "failed to generate points", - ); + let points = + generate_random_points_in_range_seeded::(n_points, wide_bounds(), seed); vertices_from_points(points) }, |vertices| { @@ -470,10 +473,7 @@ fn bench_validation(c: &mut Criterion, dimension_name: &str, n_p } let seed = seed_for_case::(n_points); - let points = bench_result( - generate_random_points_seeded::(n_points, (-100.0, 100.0), seed), - "failed to generate points", - ); + let points = generate_random_points_in_range_seeded::(n_points, wide_bounds(), seed); let vertices = vertices_from_points(points); let dt = construct_triangulation::(&vertices, seed); let tri = dt.as_triangulation(); @@ -509,10 +509,7 @@ fn bench_neighbor_queries( } let seed = seed_for_case::(n_points); - let points = bench_result( - generate_random_points_seeded::(n_points, (-100.0, 100.0), seed), - "failed to generate points", - ); + let points = generate_random_points_in_range_seeded::(n_points, wide_bounds(), seed); let vertices = vertices_from_points(points); let dt = construct_triangulation::(&vertices, seed); let tds = dt.tds(); @@ -549,10 +546,7 @@ fn bench_vertex_iteration( } let seed = seed_for_case::(n_points); - let points = bench_result( - generate_random_points_seeded::(n_points, (-100.0, 100.0), seed), - "failed to generate points", - ); + let points = generate_random_points_in_range_seeded::(n_points, wide_bounds(), seed); let vertices = vertices_from_points(points); let dt = construct_triangulation::(&vertices, seed); let tds = dt.tds(); @@ -587,10 +581,7 @@ fn bench_simplex_iteration( } let seed = seed_for_case::(n_points); - let points = bench_result( - generate_random_points_seeded::(n_points, (-100.0, 100.0), seed), - "failed to generate points", - ); + let points = generate_random_points_in_range_seeded::(n_points, wide_bounds(), seed); let vertices = vertices_from_points(points); let dt = construct_triangulation::(&vertices, seed); let tds = dt.tds(); diff --git a/benches/remove_vertex.rs b/benches/remove_vertex.rs index f7ec6f93..e3527d1b 100644 --- a/benches/remove_vertex.rs +++ b/benches/remove_vertex.rs @@ -17,8 +17,8 @@ use criterion::{BatchSize, BenchmarkId, Criterion, Throughput, criterion_group, criterion_main}; use delaunay::prelude::construction::{DelaunayTriangulation, Vertex}; -use delaunay::prelude::generators::generate_random_points_seeded; -use delaunay::prelude::geometry::{AdaptiveKernel, Coordinate, Point}; +use delaunay::prelude::generators::generate_random_points_in_range_seeded; +use delaunay::prelude::geometry::{AdaptiveKernel, Coordinate, CoordinateRange, Point}; use delaunay::prelude::tds::VertexKey; use std::hint::black_box; use std::time::Duration; @@ -28,7 +28,6 @@ use std::time::Duration; pub mod bench_utils; use bench_utils::{bench_option, bench_result}; -const INTERIOR_BOUNDS: (f64, f64) = (0.0, 1.0); const INTERIOR_RADIUS_MIN: f64 = 0.15; const INTERIOR_RADIUS_SPAN: f64 = 0.70; const NEAR_BOUNDARY_EPSILON: f64 = 1.0e-9; @@ -36,6 +35,13 @@ const NEAR_DEGENERATE_EPSILON: f64 = 1.0e-10; const COSPHERICAL_CENTER: f64 = 0.5; const COSPHERICAL_RADIUS: f64 = 0.25; const LARGE_COORDINATE_SCALE: f64 = 1.0e6; + +fn interior_bounds() -> CoordinateRange { + bench_result( + CoordinateRange::try_new(0.0_f64, 1.0), + "interior benchmark bounds must be valid", + ) +} const LARGE_COORDINATE_JITTER: f64 = 1.0e3; const SEED_SALT: u64 = 0x9E37_79B9_7F4A_7C15; const SEED_SEARCH_ATTEMPTS: usize = 64; @@ -139,10 +145,8 @@ fn generate_vertices( /// Generate well-conditioned interior points inside the canonical simplex. fn generate_interior_points(count: usize, seed: u64) -> Vec> { - let raw_points = bench_result( - generate_random_points_seeded::(count, INTERIOR_BOUNDS, seed), - format!("failed to generate {D}D interior benchmark points"), - ); + let raw_points = + generate_random_points_in_range_seeded::(count, interior_bounds(), seed); let mut points = Vec::with_capacity(count); for (index, raw_point) in raw_points.iter().enumerate() { @@ -160,10 +164,8 @@ fn generate_interior_points(count: usize, seed: u64) -> Vec(count: usize, seed: u64) -> Vec> { - let raw_points = bench_result( - generate_random_points_seeded::(count, INTERIOR_BOUNDS, seed), - format!("failed to generate {D}D near-boundary benchmark points"), - ); + let raw_points = + generate_random_points_in_range_seeded::(count, interior_bounds(), seed); let mut points = Vec::with_capacity(count); for (index, raw_point) in raw_points.iter().enumerate() { @@ -183,10 +185,8 @@ fn generate_near_boundary_points(count: usize, seed: u64) -> Vec /// Generate points on a shared sphere to stress cospherical predicates. fn generate_cospherical_points(count: usize, seed: u64) -> Vec> { - let raw_points = bench_result( - generate_random_points_seeded::(count, INTERIOR_BOUNDS, seed), - format!("failed to generate {D}D cospherical benchmark points"), - ); + let raw_points = + generate_random_points_in_range_seeded::(count, interior_bounds(), seed); let mut points = Vec::with_capacity(count); for raw_point in &raw_points { @@ -227,10 +227,8 @@ fn generate_near_degenerate_simplex(count: usize, seed: u64) -> /// Generate finite points with large coordinates to stress scale-sensitive paths. fn generate_large_coordinate_points(count: usize, seed: u64) -> Vec> { - let raw_points = bench_result( - generate_random_points_seeded::(count, INTERIOR_BOUNDS, seed), - format!("failed to generate {D}D large-coordinate benchmark points"), - ); + let raw_points = + generate_random_points_in_range_seeded::(count, interior_bounds(), seed); let mut points = Vec::with_capacity(count); for (index, raw_point) in raw_points.iter().enumerate() { diff --git a/benches/tds_clone.rs b/benches/tds_clone.rs index a2faaee7..ab890add 100644 --- a/benches/tds_clone.rs +++ b/benches/tds_clone.rs @@ -16,8 +16,8 @@ use criterion::{BenchmarkId, Criterion, Throughput, criterion_group, criterion_main}; use delaunay::prelude::construction::{DelaunayTriangulation, Vertex}; -use delaunay::prelude::generators::generate_random_points_seeded; -use delaunay::prelude::geometry::AdaptiveKernel; +use delaunay::prelude::generators::generate_random_points_in_range_seeded; +use delaunay::prelude::geometry::{AdaptiveKernel, CoordinateRange}; use delaunay::prelude::tds::Tds; use std::hint::black_box; use std::time::Duration; @@ -27,7 +27,6 @@ use std::time::Duration; pub mod bench_utils; use bench_utils::bench_result; -const BOUNDS: (f64, f64) = (-100.0, 100.0); const SEED_SALT: u64 = 0x9E37_79B9_7F4A_7C15; const SAMPLE_SIZE: usize = 10; const WARM_UP_TIME: Duration = Duration::from_millis(500); @@ -35,6 +34,13 @@ const MEASUREMENT_TIME: Duration = Duration::from_secs(2); type BenchTriangulation = DelaunayTriangulation, (), (), D>; +fn benchmark_bounds() -> CoordinateRange { + bench_result( + CoordinateRange::try_new(-100.0_f64, 100.0), + "clone benchmark bounds must be valid", + ) +} + struct CloneSource { vertex_count: usize, simplex_count: usize, @@ -56,9 +62,10 @@ fn generate_vertices( requested_vertices: usize, seed: u64, ) -> Vec> { - let points = bench_result( - generate_random_points_seeded::(requested_vertices, BOUNDS, seed), - format!("failed to generate {D}D benchmark points"), + let points = generate_random_points_in_range_seeded::( + requested_vertices, + benchmark_bounds(), + seed, ); Vertex::from_points(&points) } diff --git a/benches/topology_guarantee_construction.rs b/benches/topology_guarantee_construction.rs index 2145b5cb..8da08113 100644 --- a/benches/topology_guarantee_construction.rs +++ b/benches/topology_guarantee_construction.rs @@ -16,7 +16,8 @@ use criterion::{BenchmarkId, Criterion, Throughput, criterion_group, criterion_main}; use delaunay::prelude::construction::{DelaunayTriangulation, TopologyGuarantee}; -use delaunay::prelude::generators::generate_random_points_seeded; +use delaunay::prelude::generators::generate_random_points_in_range_seeded; +use delaunay::prelude::geometry::CoordinateRange; use delaunay::prelude::repair::DelaunayRepairPolicy; use delaunay::prelude::validation::ValidationPolicy; use delaunay::vertex; @@ -28,9 +29,15 @@ use std::time::Duration; pub mod bench_utils; use bench_utils::{abort_benchmark, bench_result}; -const BOUNDS: (f64, f64) = (-100.0, 100.0); const SEED_SALT: u64 = 0x9E37_79B9_7F4A_7C15; +fn benchmark_bounds() -> CoordinateRange { + bench_result( + CoordinateRange::try_new(-100.0_f64, 100.0), + "topology-guarantee benchmark bounds must be valid", + ) +} + fn bench_dimension( c: &mut Criterion, dim_label: &str, @@ -48,10 +55,8 @@ fn bench_dimension( // Deterministic input per (dimension, count). let seed = seed_base ^ (n_points as u64).wrapping_mul(SEED_SALT); - let points = bench_result( - generate_random_points_seeded::(n_points, BOUNDS, seed), - "failed to generate benchmark points", - ); + let points = + generate_random_points_in_range_seeded::(n_points, benchmark_bounds(), seed); let vertices = points.into_iter().map(|p| vertex!(p)).collect::>(); group.bench_with_input( diff --git a/examples/triangulation_and_hull.rs b/examples/triangulation_and_hull.rs index 8a7f066d..2cae74bb 100644 --- a/examples/triangulation_and_hull.rs +++ b/examples/triangulation_and_hull.rs @@ -18,8 +18,8 @@ use delaunay::prelude::construction::{ ConstructionOptions, DelaunayTriangulation, DelaunayTriangulationConstructionError, RetryPolicy, vertex, }; -use delaunay::prelude::generators::{RandomPointGenerationError, generate_random_points_seeded}; -use delaunay::prelude::geometry::AdaptiveKernel; +use delaunay::prelude::generators::generate_random_points_in_range_seeded; +use delaunay::prelude::geometry::{AdaptiveKernel, CoordinateRange, CoordinateRangeError}; use delaunay::prelude::query::{ AdjacencyIndexBuildError, ConvexHull, ConvexHullConstructionError, Coordinate, Point, QueryError, @@ -30,7 +30,7 @@ type WorkflowTriangulation = DelaunayTriangulation for WorkflowExampleError { } fn main() -> Result<(), WorkflowExampleError> { - run_case::<3>("3D", 750, 873, (-100.0, 100.0))?; + let bounds = CoordinateRange::try_new(-100.0_f64, 100.0)?; + run_case::<3>("3D", 750, 873, bounds)?; println!(); - run_case::<4>("4D", 75, 531, (-100.0, 100.0))?; + run_case::<4>("4D", 75, 531, bounds)?; Ok(()) } @@ -67,9 +68,9 @@ fn run_case( label: &str, point_count: usize, seed: u64, - bounds: (f64, f64), + bounds: CoordinateRange, ) -> Result<(), WorkflowExampleError> { - let points = generate_random_points_seeded::(point_count, bounds, seed)?; + let points = generate_random_points_in_range_seeded::(point_count, bounds, seed); let vertices = points .iter() .map(|point| vertex!(*point)) @@ -95,7 +96,7 @@ fn run_case( println!(" hull facets: {}", hull.number_of_facets()); let inside = centroid_point(&points)?; - let outside = Point::new([bounds.1 * 2.5; D]); + let outside = Point::new([bounds.max() * 2.5; D]); println!( " hull query: centroid outside? {}", diff --git a/src/core/util/hilbert.rs b/src/core/util/hilbert.rs index f2ee606e..f680f9fd 100644 --- a/src/core/util/hilbert.rs +++ b/src/core/util/hilbert.rs @@ -129,6 +129,39 @@ pub enum HilbertError { /// The out-of-range pre-quantized coordinate value. coordinate: u32, }, + + /// An internally constructed Hilbert sort permutation had the wrong length. + #[error( + "Hilbert sort permutation length mismatch: item count {item_count}, permutation count {permutation_count}" + )] + InvalidSortPermutationLength { + /// Number of items being sorted. + item_count: usize, + /// Number of indices in the permutation. + permutation_count: usize, + }, + + /// An internally constructed Hilbert sort permutation referenced an invalid item index. + #[error( + "Hilbert sort permutation index {permutation_index} has value {item_index}, which is outside item count {item_count}" + )] + InvalidSortPermutationIndex { + /// Position in the permutation whose value was invalid. + permutation_index: usize, + /// Invalid source item index. + item_index: usize, + /// Number of items being sorted. + item_count: usize, + }, + + /// An internally constructed Hilbert sort permutation referenced the same item twice. + #[error("Hilbert sort permutation index {permutation_index} repeats item index {item_index}")] + InvalidSortPermutationDuplicate { + /// Position in the permutation whose value repeated an earlier item. + permutation_index: usize, + /// Repeated source item index. + item_index: usize, + }, } /// Validated Hilbert bit depth per coordinate. @@ -221,6 +254,305 @@ impl fmt::Display for HilbertBitDepth { } } +/// Pre-quantized Hilbert coordinates proven to fit a selected bit-depth grid. +/// +/// This borrowed wrapper carries the validation evidence for a batch of +/// caller-supplied quantized coordinates. Use [`Self::try_new`] at the boundary, +/// then call [`Self::indices`] or [`hilbert_indices_for_quantized_batch`] for an +/// infallible mapping to Hilbert indices. +/// +/// # Examples +/// +/// ```rust +/// use delaunay::prelude::ordering::{ +/// HilbertBitDepth, HilbertError, HilbertQuantizedBatch, +/// }; +/// +/// let quantized = [[0_u32, 0], [3, 3]]; +/// let bits = HilbertBitDepth::try_new(2)?; +/// let batch = HilbertQuantizedBatch::try_new(&quantized, bits)?; +/// +/// let indices = batch.indices(); +/// assert_eq!(indices.len(), quantized.len()); +/// # Ok::<(), HilbertError>(()) +/// ``` +#[derive(Clone, Copy, Debug)] +#[must_use] +pub struct HilbertQuantizedBatch<'a, const D: usize> { + quantized: &'a [[u32; D]], + index_mode: HilbertIndexMode, +} + +impl<'a, const D: usize> HilbertQuantizedBatch<'a, D> { + /// Parses pre-quantized coordinates into a validated Hilbert batch. + /// + /// # Errors + /// + /// Returns [`HilbertError::IndexOverflow`] if `D * bits > 128` (index would not fit in `u128`). + /// Returns [`HilbertError::DimensionTooLarge`] if the dimension `D` exceeds `u32::MAX`. + /// Returns [`HilbertError::PrequantizedCoordinateOutOfRange`] if any pre-quantized + /// coordinate exceeds `2^bits - 1`. + /// + /// # Examples + /// + /// ```rust + /// use delaunay::prelude::ordering::{ + /// HilbertBitDepth, HilbertError, HilbertQuantizedBatch, + /// }; + /// + /// let quantized = [[0_u32, 0], [1, 2], [3, 3]]; + /// let bits = HilbertBitDepth::try_new(2)?; + /// + /// let batch = HilbertQuantizedBatch::try_new(&quantized, bits)?; + /// assert_eq!(batch.coordinates(), quantized.as_slice()); + /// # Ok::<(), HilbertError>(()) + /// ``` + pub fn try_new(quantized: &'a [[u32; D]], bits: HilbertBitDepth) -> Result { + let index_mode = HilbertIndexMode::try_new(bits)?; + + if D != 0 { + validate_prequantized_coordinates(quantized, bits)?; + } + + Ok(Self { + quantized, + index_mode, + }) + } + + /// Returns the validated pre-quantized coordinates. + /// + /// # Examples + /// + /// ```rust + /// use delaunay::prelude::ordering::{ + /// HilbertBitDepth, HilbertError, HilbertQuantizedBatch, + /// }; + /// + /// let quantized = [[0_u32, 0], [3, 3]]; + /// let bits = HilbertBitDepth::try_new(2)?; + /// let batch = HilbertQuantizedBatch::try_new(&quantized, bits)?; + /// + /// assert_eq!(batch.coordinates(), quantized.as_slice()); + /// # Ok::<(), HilbertError>(()) + /// ``` + #[must_use] + pub const fn coordinates(self) -> &'a [[u32; D]] { + self.quantized + } + + /// Returns the bit depth whose grid bounds were checked. + /// + /// # Examples + /// + /// ```rust + /// use delaunay::prelude::ordering::{ + /// HilbertBitDepth, HilbertError, HilbertQuantizedBatch, + /// }; + /// + /// let quantized = [[0_u32, 0], [3, 3]]; + /// let bits = HilbertBitDepth::try_new(2)?; + /// let batch = HilbertQuantizedBatch::try_new(&quantized, bits)?; + /// + /// assert_eq!(batch.bits(), bits); + /// # Ok::<(), HilbertError>(()) + /// ``` + pub const fn bits(self) -> HilbertBitDepth { + self.index_mode.bits() + } + + /// Computes Hilbert indices without revalidating the batch. + /// + /// # Examples + /// + /// ```rust + /// use delaunay::prelude::ordering::{ + /// HilbertBitDepth, HilbertError, HilbertQuantizedBatch, + /// }; + /// + /// let quantized = [[0_u32, 0], [3, 3]]; + /// let bits = HilbertBitDepth::try_new(2)?; + /// let batch = HilbertQuantizedBatch::try_new(&quantized, bits)?; + /// + /// let indices = batch.indices(); + /// assert_eq!(indices.len(), quantized.len()); + /// # Ok::<(), HilbertError>(()) + /// ``` + #[must_use] + pub fn indices(self) -> Vec { + hilbert_indices_for_quantized_batch(self) + } +} + +/// Owned, pre-quantized Hilbert coordinates proven in-grid by the quantizer +/// that produced them. +/// +/// Unlike [`HilbertQuantizedBatch`], which borrows caller-supplied coordinates +/// and revalidates them at the boundary, this type is only constructed by +/// [`hilbert_quantize_batch_in_range`]. That constructor clamps every +/// coordinate into the selected bit-depth grid, so the in-grid invariant is +/// carried structurally by the stored, validated index mode and +/// [`Self::indices`] is infallible — no second per-coordinate scan is needed. +/// +/// # Examples +/// +/// ```rust +/// use delaunay::prelude::geometry::CoordinateRange; +/// use delaunay::prelude::ordering::{ +/// HilbertBitDepth, HilbertError, hilbert_quantize_batch_in_range, +/// }; +/// +/// # fn main() -> Result<(), HilbertError> { +/// let points = [[0.1_f64, 0.2], [0.9, 0.8], [0.5, 0.5]]; +/// # let Ok(bounds) = CoordinateRange::try_new(0.0_f64, 1.0) else { +/// # return Ok(()); +/// # }; +/// let bits = HilbertBitDepth::try_new(8)?; +/// +/// let batch = hilbert_quantize_batch_in_range(&points, bounds, bits, |p| *p)?; +/// assert_eq!(batch.indices().len(), points.len()); +/// # Ok(()) +/// # } +/// ``` +#[derive(Clone, Debug, PartialEq, Eq)] +#[must_use] +pub struct HilbertQuantizedVec { + quantized: Vec<[u32; D]>, + index_mode: HilbertIndexMode, +} + +impl HilbertQuantizedVec { + /// Returns the validated pre-quantized coordinates. + #[must_use] + pub fn coordinates(&self) -> &[[u32; D]] { + &self.quantized + } + + /// Returns the bit depth whose grid the coordinates were quantized to. + pub const fn bits(&self) -> HilbertBitDepth { + self.index_mode.bits() + } + + /// Returns the number of quantized points in the batch. + #[must_use] + pub const fn len(&self) -> usize { + self.quantized.len() + } + + /// Returns `true` when the batch contains no points. + #[must_use] + pub const fn is_empty(&self) -> bool { + self.quantized.is_empty() + } + + /// Computes Hilbert indices without revalidating the batch. + /// + /// This is infallible: the constructor already validated the index width + /// and clamped every coordinate into the selected bit-depth grid. + #[must_use] + pub fn indices(&self) -> Vec { + indices_for_mode(&self.quantized, self.index_mode) + } + + /// Consumes the batch, returning its Hilbert indices alongside the owned + /// quantized coordinates. + /// + /// This fuses the two products batch callers usually need — the per-point + /// Hilbert index and the quantized cell used as a sort tie-break — without + /// an extra allocation or a per-coordinate revalidation pass. + #[must_use] + pub fn into_indices_and_coordinates(self) -> (Vec, Vec<[u32; D]>) { + let indices = self.indices(); + (indices, self.quantized) + } + + /// Consumes the batch, returning only the owned quantized coordinates. + #[must_use] + pub fn into_coordinates(self) -> Vec<[u32; D]> { + self.quantized + } +} + +/// Quantizes a batch of items into an owned, proof-bearing Hilbert batch. +/// +/// Coordinates are extracted with `coords_of`, normalized against `bounds` +/// (already parsed at an upstream boundary), and clamped into the +/// `0..=2^bits - 1` grid. The index width and quantization scale are validated +/// once for the whole batch, after which [`HilbertQuantizedVec::indices`] is +/// infallible. +/// +/// This is the single-pass bulk entry point preferred by construction +/// preprocessing: it avoids both the per-item bound parsing of +/// [`try_hilbert_quantize`] and the per-coordinate revalidation that +/// [`hilbert_indices_prequantized`] performs on caller-supplied grids. +/// +/// # Errors +/// +/// Returns [`HilbertError::IndexOverflow`] if `D * bits > 128` (index would not fit in `u128`). +/// Returns [`HilbertError::DimensionTooLarge`] if the dimension `D` exceeds `u32::MAX`. +/// Returns [`HilbertError::QuantizationScaleConversionFailed`] if the quantization +/// grid maximum cannot be represented by the coordinate scalar type. Returns +/// [`HilbertError::NonFiniteBoundsExtent`], [`HilbertError::NonFiniteCoordinate`], +/// [`HilbertError::NonFiniteNormalizedCoordinate`], or +/// [`HilbertError::QuantizedCoordinateConversionFailed`] if coordinate +/// quantization fails. +/// +/// # Examples +/// +/// ```rust +/// use delaunay::prelude::geometry::CoordinateRange; +/// use delaunay::prelude::ordering::{ +/// HilbertBitDepth, HilbertError, hilbert_quantize_batch_in_range, +/// }; +/// +/// # fn main() -> Result<(), HilbertError> { +/// let points = [[0.0_f64, 0.0], [1.0, 1.0]]; +/// # let Ok(bounds) = CoordinateRange::try_new(0.0_f64, 1.0) else { +/// # return Ok(()); +/// # }; +/// let batch = +/// hilbert_quantize_batch_in_range(&points, bounds, HilbertBitDepth::try_new(4)?, |p| *p)?; +/// +/// let (indices, quantized) = batch.into_indices_and_coordinates(); +/// assert_eq!(indices.len(), 2); +/// assert_eq!(quantized.len(), 2); +/// # Ok(()) +/// # } +/// ``` +pub fn hilbert_quantize_batch_in_range( + items: &[Item], + bounds: CoordinateRange, + bits: HilbertBitDepth, + mut coords_of: impl FnMut(&Item) -> [T; D], +) -> Result, HilbertError> +where + T: CoordinateScalar, +{ + let index_mode = HilbertIndexMode::try_new(bits)?; + + if D == 0 { + return Ok(HilbertQuantizedVec { + quantized: vec![[0_u32; D]; items.len()], + index_mode, + }); + } + + let (max_val_u32, max_val_t) = quantization_scale::(bits)?; + + let quantized = items + .iter() + .map(|item| { + let coords = coords_of(item); + quantize_with_scale(&coords, bounds, bits, max_val_u32, max_val_t) + }) + .collect::, HilbertError>>()?; + + Ok(HilbertQuantizedVec { + quantized, + index_mode, + }) +} + /// Converts a validated bit depth into a scalar grid maximum so Hilbert /// ordering cannot silently collapse when a coordinate type cannot represent it. fn quantization_scale( @@ -251,8 +583,45 @@ fn total_bits(bits: HilbertBitDepth) -> Result { + bits: HilbertBitDepth, +} + +impl HilbertIndexParams { + const fn bits(self) -> u32 { + self.bits.get() + } +} + +/// Validated Hilbert indexing mode, including the zero-dimensional special case. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +enum HilbertIndexMode { + ZeroDimensional { bits: HilbertBitDepth }, + Positive(HilbertIndexParams), +} + +impl HilbertIndexMode { + fn try_new(bits: HilbertBitDepth) -> Result { + validate_index_width::(bits)?; + if D == 0 { + Ok(Self::ZeroDimensional { bits }) + } else { + Ok(Self::Positive(HilbertIndexParams { bits })) + } + } + + const fn bits(self) -> HilbertBitDepth { + match self { + Self::ZeroDimensional { bits } => bits, + Self::Positive(params) => params.bits, + } + } +} + /// Centralizes index-width validation shared by indexing and ordering APIs. -fn validate_index_params(bits: HilbertBitDepth) -> Result<(), HilbertError> { +fn validate_index_width(bits: HilbertBitDepth) -> Result<(), HilbertError> { let total_bits = total_bits::(bits)?; if total_bits > 128 { return Err(HilbertError::IndexOverflow { @@ -304,28 +673,72 @@ fn parse_hilbert_bounds( /// # Examples /// /// ```rust -/// use delaunay::prelude::ordering::{HilbertBitDepth, HilbertError, hilbert_quantize}; +/// use delaunay::prelude::ordering::{HilbertBitDepth, HilbertError, try_hilbert_quantize}; /// /// let coords = [0.5_f64, 0.25]; -/// let q = hilbert_quantize(&coords, (0.0, 1.0), HilbertBitDepth::try_new(2)?)?; +/// let q = try_hilbert_quantize(&coords, (0.0, 1.0), HilbertBitDepth::try_new(2)?)?; /// assert!(q[0] <= 3 && q[1] <= 3); /// # Ok::<(), HilbertError>(()) /// ``` -pub fn hilbert_quantize( +pub fn try_hilbert_quantize( coords: &[T; D], bounds: (T, T), bits: HilbertBitDepth, ) -> Result<[u32; D], HilbertError> { + let bounds = parse_hilbert_bounds(bounds)?; + if D == 0 { return Ok([0_u32; D]); } let (max_val_u32, max_val_t) = quantization_scale::(bits)?; - let bounds = parse_hilbert_bounds(bounds)?; quantize_with_scale(coords, bounds, bits, max_val_u32, max_val_t) } +/// Quantizes coordinates against bounds already parsed by an upstream boundary. +/// +/// # Errors +/// +/// Returns [`HilbertError::QuantizationScaleConversionFailed`] if the +/// quantization grid maximum cannot be represented by the coordinate scalar +/// type. Returns [`HilbertError::NonFiniteBoundsExtent`] if the validated bounds +/// produce a non-finite extent. Returns [`HilbertError::NonFiniteCoordinate`] or +/// [`HilbertError::NonFiniteNormalizedCoordinate`] if a coordinate or its +/// normalized value is non-finite. Returns +/// [`HilbertError::QuantizedCoordinateConversionFailed`] if a rounded scaled +/// coordinate cannot be represented as `u32`. +/// +/// # Examples +/// +/// ```rust +/// use delaunay::prelude::geometry::CoordinateRange; +/// use delaunay::prelude::ordering::{ +/// HilbertBitDepth, HilbertError, hilbert_quantize_in_range, +/// }; +/// +/// # fn main() -> Result<(), HilbertError> { +/// # let Ok(bounds) = CoordinateRange::try_new(0.0_f64, 1.0) else { +/// # return Ok(()); +/// # }; +/// let q = hilbert_quantize_in_range(&[0.5_f64, 0.25], bounds, HilbertBitDepth::try_new(2)?)?; +/// assert!(q[0] <= 3 && q[1] <= 3); +/// # Ok(()) +/// # } +/// ``` +pub fn hilbert_quantize_in_range( + coords: &[T; D], + bounds: CoordinateRange, + bits: HilbertBitDepth, +) -> Result<[u32; D], HilbertError> { + if D == 0 { + return Ok([0_u32; D]); + } + + let (max_val_u32, max_val_t) = quantization_scale::(bits)?; + quantize_with_scale(coords, bounds, bits, max_val_u32, max_val_t) +} + /// Quantizes coordinates with a precomputed scalar grid maximum so hot callers /// can validate conversion once before sorting or batch index generation. #[inline] @@ -378,21 +791,60 @@ fn quantize_with_scale( /// Applies a prevalidated permutation after key construction succeeds so sort /// helpers never partially reorder items before returning a Hilbert error. -fn apply_order(items: &mut [Item], order: Vec) { - debug_assert_eq!(items.len(), order.len()); +fn apply_order( + items: &mut [Item], + order: impl ExactSizeIterator, +) -> Result<(), HilbertError> { + let item_len = items.len(); + let permutation_count = order.len(); + if item_len != permutation_count { + return Err(HilbertError::InvalidSortPermutationLength { + item_count: item_len, + permutation_count, + }); + } - let mut ranks = vec![0_usize; order.len()]; + let mut ranks = vec![usize::MAX; item_len]; + let mut observed_count = 0_usize; for (new_index, old_index) in order.into_iter().enumerate() { + observed_count = new_index + 1; + if new_index >= item_len { + return Err(HilbertError::InvalidSortPermutationLength { + item_count: item_len, + permutation_count: observed_count, + }); + } + if old_index >= item_len { + return Err(HilbertError::InvalidSortPermutationIndex { + permutation_index: new_index, + item_index: old_index, + item_count: item_len, + }); + } + if ranks[old_index] != usize::MAX { + return Err(HilbertError::InvalidSortPermutationDuplicate { + permutation_index: new_index, + item_index: old_index, + }); + } ranks[old_index] = new_index; } + if observed_count != item_len { + return Err(HilbertError::InvalidSortPermutationLength { + item_count: item_len, + permutation_count: observed_count, + }); + } - for index in 0..items.len() { + for index in 0..item_len { while ranks[index] != index { let target = ranks[index]; items.swap(index, target); ranks.swap(index, target); } } + + Ok(()) } /// Compute the Hilbert curve index for a point in D-dimensional space. @@ -423,25 +875,61 @@ fn apply_order(items: &mut [Item], order: Vec) { /// # Examples /// /// ```rust -/// use delaunay::prelude::ordering::{HilbertBitDepth, HilbertError, hilbert_index}; +/// use delaunay::prelude::ordering::{HilbertBitDepth, HilbertError, try_hilbert_index}; /// -/// let idx = hilbert_index(&[0.0_f64, 0.0], (0.0, 1.0), HilbertBitDepth::try_new(4)?)?; +/// let idx = try_hilbert_index(&[0.0_f64, 0.0], (0.0, 1.0), HilbertBitDepth::try_new(4)?)?; /// assert_eq!(idx, 0); /// # Ok::<(), HilbertError>(()) /// ``` -pub fn hilbert_index( +pub fn try_hilbert_index( coords: &[T; D], bounds: (T, T), bits: HilbertBitDepth, ) -> Result { - validate_index_params::(bits)?; + let bounds = parse_hilbert_bounds(bounds)?; + hilbert_index_in_range(coords, bounds, bits) +} - if D == 0 { +/// Computes a Hilbert index against bounds already parsed by an upstream boundary. +/// +/// # Errors +/// +/// Returns [`HilbertError::IndexOverflow`] if `D * bits > 128` (index would not fit in `u128`). +/// Returns [`HilbertError::DimensionTooLarge`] if the dimension `D` exceeds `u32::MAX`. +/// Returns [`HilbertError::QuantizationScaleConversionFailed`] if the quantization +/// grid maximum cannot be represented by the coordinate scalar type. Returns +/// [`HilbertError::NonFiniteBoundsExtent`], [`HilbertError::NonFiniteCoordinate`], +/// [`HilbertError::NonFiniteNormalizedCoordinate`], or +/// [`HilbertError::QuantizedCoordinateConversionFailed`] if coordinate +/// quantization fails. +/// +/// # Examples +/// +/// ```rust +/// use delaunay::prelude::geometry::CoordinateRange; +/// use delaunay::prelude::ordering::{HilbertBitDepth, HilbertError, hilbert_index_in_range}; +/// +/// # fn main() -> Result<(), HilbertError> { +/// # let Ok(bounds) = CoordinateRange::try_new(0.0_f64, 1.0) else { +/// # return Ok(()); +/// # }; +/// let idx = hilbert_index_in_range(&[0.0_f64, 0.0], bounds, HilbertBitDepth::try_new(4)?)?; +/// assert_eq!(idx, 0); +/// # Ok(()) +/// # } +/// ``` +pub fn hilbert_index_in_range( + coords: &[T; D], + bounds: CoordinateRange, + bits: HilbertBitDepth, +) -> Result { + let index_mode = HilbertIndexMode::try_new(bits)?; + let HilbertIndexMode::Positive(index_params) = index_mode else { return Ok(0); - } + }; - let q = hilbert_quantize(coords, bounds, bits)?; - Ok(index_from_quantized(&q, bits)) + let q = hilbert_quantize_in_range(coords, bounds, bits)?; + Ok(index_from_quantized(&q, index_params)) } /// Compute Hilbert index from pre-quantized integer coordinates. @@ -452,17 +940,8 @@ pub fn hilbert_index( /// The resulting ordering is continuous on the integer grid (successive indices move to /// adjacent cells). #[must_use] -fn index_from_quantized(coords: &[u32; D], bits: HilbertBitDepth) -> u128 { - let bits = bits.get(); - debug_assert!(D > 0, "caller should handle D==0"); - debug_assert!( - bits > 0 && bits <= 31, - "bits must be in range [1, 31], got {bits}" - ); - debug_assert!( - (D as u128) * u128::from(bits) <= 128, - "Hilbert index would overflow u128 for D={D} and bits={bits}" - ); +fn index_from_quantized(coords: &[u32; D], params: HilbertIndexParams) -> u128 { + let bits = params.bits(); // Work on a local copy in "transposed" form. let mut transposed = *coords; @@ -554,30 +1033,82 @@ fn index_from_quantized(coords: &[u32; D], bits: HilbertBitDepth /// Returns [`HilbertError::QuantizedCoordinateConversionFailed`] if a rounded /// scaled coordinate cannot be represented as `u32`. /// +/// Returns [`HilbertError::InvalidSortPermutationLength`], +/// [`HilbertError::InvalidSortPermutationIndex`], or +/// [`HilbertError::InvalidSortPermutationDuplicate`] if an internally +/// constructed permutation is inconsistent. +/// /// # Examples /// /// ```rust -/// use delaunay::prelude::ordering::{HilbertBitDepth, HilbertError, hilbert_sort_by_stable}; +/// use delaunay::prelude::ordering::{HilbertBitDepth, HilbertError, try_hilbert_sort_by_stable}; /// /// let mut points = vec![[0.9_f64, 0.9], [0.1, 0.1], [0.5, 0.5]]; -/// hilbert_sort_by_stable(&mut points, (0.0, 1.0), HilbertBitDepth::try_new(8)?, |p| *p)?; +/// try_hilbert_sort_by_stable(&mut points, (0.0, 1.0), HilbertBitDepth::try_new(8)?, |p| *p)?; /// assert_eq!(points[0], [0.1, 0.1]); /// # Ok::<(), HilbertError>(()) /// ``` -pub fn hilbert_sort_by_stable( +pub fn try_hilbert_sort_by_stable( items: &mut [Item], bounds: (T, T), bits: HilbertBitDepth, - mut coords_of: impl FnMut(&Item) -> [T; D], + coords_of: impl FnMut(&Item) -> [T; D], ) -> Result<(), HilbertError> { - validate_index_params::(bits)?; + let bounds = parse_hilbert_bounds(bounds)?; + hilbert_sort_by_stable_in_range(items, bounds, bits, coords_of) +} - if D == 0 { +/// Stable sort helper using bounds already parsed by an upstream boundary. +/// +/// This is equivalent to [`try_hilbert_sort_by_stable`], but accepts a +/// [`CoordinateRange`] so callers can carry range validation evidence inward. +/// +/// # Errors +/// +/// Returns [`HilbertError::IndexOverflow`] if `D * bits > 128` (index would not fit in `u128`). +/// Returns [`HilbertError::DimensionTooLarge`] if the dimension `D` exceeds `u32::MAX`. +/// Returns [`HilbertError::QuantizationScaleConversionFailed`] if the quantization +/// grid maximum cannot be represented by the coordinate scalar type. Returns +/// [`HilbertError::NonFiniteBoundsExtent`], [`HilbertError::NonFiniteCoordinate`], +/// [`HilbertError::NonFiniteNormalizedCoordinate`], or +/// [`HilbertError::QuantizedCoordinateConversionFailed`] if coordinate +/// quantization fails. +/// +/// Returns [`HilbertError::InvalidSortPermutationLength`], +/// [`HilbertError::InvalidSortPermutationIndex`], or +/// [`HilbertError::InvalidSortPermutationDuplicate`] if an internally +/// constructed permutation is inconsistent. +/// +/// # Examples +/// +/// ```rust +/// use delaunay::prelude::geometry::CoordinateRange; +/// use delaunay::prelude::ordering::{ +/// HilbertBitDepth, HilbertError, hilbert_sort_by_stable_in_range, +/// }; +/// +/// # fn main() -> Result<(), HilbertError> { +/// let mut points = vec![[0.9_f64, 0.9], [0.1, 0.1], [0.5, 0.5]]; +/// # let Ok(bounds) = CoordinateRange::try_new(0.0_f64, 1.0) else { +/// # return Ok(()); +/// # }; +/// hilbert_sort_by_stable_in_range(&mut points, bounds, HilbertBitDepth::try_new(8)?, |p| *p)?; +/// assert_eq!(points[0], [0.1, 0.1]); +/// # Ok(()) +/// # } +/// ``` +pub fn hilbert_sort_by_stable_in_range( + items: &mut [Item], + bounds: CoordinateRange, + bits: HilbertBitDepth, + mut coords_of: impl FnMut(&Item) -> [T; D], +) -> Result<(), HilbertError> { + let index_mode = HilbertIndexMode::try_new(bits)?; + let HilbertIndexMode::Positive(index_params) = index_mode else { return Ok(()); - } + }; let (max_val_u32, max_val_t) = quantization_scale::(bits)?; - let bounds = parse_hilbert_bounds(bounds)?; let mut keyed: Vec<((u128, [u32; D]), usize)> = items .iter() @@ -585,13 +1116,13 @@ pub fn hilbert_sort_by_stable( .map(|(i, item)| { let c = coords_of(item); let q = quantize_with_scale(&c, bounds, bits, max_val_u32, max_val_t)?; - let idx = index_from_quantized(&q, bits); + let idx = index_from_quantized(&q, index_params); Ok(((idx, q), i)) }) .collect::>()?; keyed.sort_by_key(|(key, _)| *key); - apply_order(items, keyed.into_iter().map(|(_, i)| i).collect()); + apply_order(items, keyed.into_iter().map(|(_, i)| i))?; Ok(()) } @@ -599,7 +1130,7 @@ pub fn hilbert_sort_by_stable( /// Unstable sort helper: sort items by Hilbert index + quantized-coordinate tie-break. /// /// This precomputes fallible Hilbert keys once, then applies an unstable ordering. -/// Prefer [`hilbert_sort_by_stable`] when equal-key items must preserve their +/// Prefer [`try_hilbert_sort_by_stable`] when equal-key items must preserve their /// original relative order. /// /// When `D == 0`, all items are considered equivalent (index 0) and the sort order is @@ -608,47 +1139,99 @@ pub fn hilbert_sort_by_stable( /// # Errors /// /// Returns [`HilbertError::IndexOverflow`] if `D * bits > 128` (index would not fit in `u128`). -/// -/// Returns [`HilbertError::DimensionTooLarge`] if the dimension `D` exceeds `u32::MAX` -/// (extremely unlikely in practice). -/// +/// +/// Returns [`HilbertError::DimensionTooLarge`] if the dimension `D` exceeds `u32::MAX` +/// (extremely unlikely in practice). +/// +/// Returns [`HilbertError::QuantizationScaleConversionFailed`] if the quantization +/// grid maximum cannot be represented by the coordinate scalar type. +/// +/// Returns [`HilbertError::NonFiniteBounds`], +/// [`HilbertError::NonIncreasingBounds`], +/// [`HilbertError::NonFiniteBoundsExtent`], +/// [`HilbertError::NonFiniteCoordinate`], or +/// [`HilbertError::NonFiniteNormalizedCoordinate`] if quantization input or +/// normalization arithmetic is non-finite. +/// +/// Returns [`HilbertError::QuantizedCoordinateConversionFailed`] if a rounded +/// scaled coordinate cannot be represented as `u32`. +/// +/// Returns [`HilbertError::InvalidSortPermutationLength`], +/// [`HilbertError::InvalidSortPermutationIndex`], or +/// [`HilbertError::InvalidSortPermutationDuplicate`] if an internally +/// constructed permutation is inconsistent. +/// +/// # Examples +/// +/// ```rust +/// use delaunay::prelude::ordering::{HilbertBitDepth, HilbertError, try_hilbert_sort_by_unstable}; +/// +/// let mut points = vec![[0.9_f64, 0.9], [0.1, 0.1], [0.5, 0.5]]; +/// try_hilbert_sort_by_unstable(&mut points, (0.0, 1.0), HilbertBitDepth::try_new(8)?, |p| *p)?; +/// assert_eq!(points[0], [0.1, 0.1]); +/// # Ok::<(), HilbertError>(()) +/// ``` +pub fn try_hilbert_sort_by_unstable( + items: &mut [Item], + bounds: (T, T), + bits: HilbertBitDepth, + coords_of: impl FnMut(&Item) -> [T; D], +) -> Result<(), HilbertError> { + let bounds = parse_hilbert_bounds(bounds)?; + hilbert_sort_by_unstable_in_range(items, bounds, bits, coords_of) +} + +/// Unstable sort helper using bounds already parsed by an upstream boundary. +/// +/// This is equivalent to [`try_hilbert_sort_by_unstable`], but accepts a +/// [`CoordinateRange`] so callers can carry range validation evidence inward. +/// +/// # Errors +/// +/// Returns [`HilbertError::IndexOverflow`] if `D * bits > 128` (index would not fit in `u128`). +/// Returns [`HilbertError::DimensionTooLarge`] if the dimension `D` exceeds `u32::MAX`. /// Returns [`HilbertError::QuantizationScaleConversionFailed`] if the quantization -/// grid maximum cannot be represented by the coordinate scalar type. -/// -/// Returns [`HilbertError::NonFiniteBounds`], -/// [`HilbertError::NonIncreasingBounds`], -/// [`HilbertError::NonFiniteBoundsExtent`], -/// [`HilbertError::NonFiniteCoordinate`], or -/// [`HilbertError::NonFiniteNormalizedCoordinate`] if quantization input or -/// normalization arithmetic is non-finite. +/// grid maximum cannot be represented by the coordinate scalar type. Returns +/// [`HilbertError::NonFiniteBoundsExtent`], [`HilbertError::NonFiniteCoordinate`], +/// [`HilbertError::NonFiniteNormalizedCoordinate`], or +/// [`HilbertError::QuantizedCoordinateConversionFailed`] if coordinate +/// quantization fails. /// -/// Returns [`HilbertError::QuantizedCoordinateConversionFailed`] if a rounded -/// scaled coordinate cannot be represented as `u32`. +/// Returns [`HilbertError::InvalidSortPermutationLength`], +/// [`HilbertError::InvalidSortPermutationIndex`], or +/// [`HilbertError::InvalidSortPermutationDuplicate`] if an internally +/// constructed permutation is inconsistent. /// /// # Examples /// /// ```rust -/// use delaunay::prelude::ordering::{HilbertBitDepth, HilbertError, hilbert_sort_by_unstable}; +/// use delaunay::prelude::geometry::CoordinateRange; +/// use delaunay::prelude::ordering::{ +/// HilbertBitDepth, HilbertError, hilbert_sort_by_unstable_in_range, +/// }; /// +/// # fn main() -> Result<(), HilbertError> { /// let mut points = vec![[0.9_f64, 0.9], [0.1, 0.1], [0.5, 0.5]]; -/// hilbert_sort_by_unstable(&mut points, (0.0, 1.0), HilbertBitDepth::try_new(8)?, |p| *p)?; +/// # let Ok(bounds) = CoordinateRange::try_new(0.0_f64, 1.0) else { +/// # return Ok(()); +/// # }; +/// hilbert_sort_by_unstable_in_range(&mut points, bounds, HilbertBitDepth::try_new(8)?, |p| *p)?; /// assert_eq!(points[0], [0.1, 0.1]); -/// # Ok::<(), HilbertError>(()) +/// # Ok(()) +/// # } /// ``` -pub fn hilbert_sort_by_unstable( +pub fn hilbert_sort_by_unstable_in_range( items: &mut [Item], - bounds: (T, T), + bounds: CoordinateRange, bits: HilbertBitDepth, mut coords_of: impl FnMut(&Item) -> [T; D], ) -> Result<(), HilbertError> { - validate_index_params::(bits)?; - - if D == 0 { + let index_mode = HilbertIndexMode::try_new(bits)?; + let HilbertIndexMode::Positive(index_params) = index_mode else { return Ok(()); - } + }; let (max_val_u32, max_val_t) = quantization_scale::(bits)?; - let bounds = parse_hilbert_bounds(bounds)?; let mut keyed: Vec<((u128, [u32; D]), usize)> = items .iter() @@ -656,13 +1239,13 @@ pub fn hilbert_sort_by_unstable( .map(|(i, item)| { let c = coords_of(item); let q = quantize_with_scale(&c, bounds, bits, max_val_u32, max_val_t)?; - let idx = index_from_quantized(&q, bits); + let idx = index_from_quantized(&q, index_params); Ok(((idx, q), i)) }) .collect::>()?; keyed.sort_unstable_by_key(|(key, _)| *key); - apply_order(items, keyed.into_iter().map(|(_, i)| i).collect()); + apply_order(items, keyed.into_iter().map(|(_, i)| i))?; Ok(()) } @@ -698,7 +1281,7 @@ fn validate_prequantized_coordinates( /// /// This is a bulk API that avoids recomputing quantization parameters for large /// insertion batches. When inserting many points, quantize them once using -/// [`hilbert_quantize`] and then call this function to compute all indices in bulk. +/// [`try_hilbert_quantize`] and then call this function to compute all indices in bulk. /// Pre-quantized coordinates must be in the inclusive range `0..=2^bits - 1`; /// values outside that grid are rejected instead of being truncated. /// @@ -706,8 +1289,10 @@ fn validate_prequantized_coordinates( /// /// This function validates index width and pre-quantized coordinate ranges, then /// maps each quantized coordinate through the internal Hilbert index computation. -/// For large batches, this is significantly faster than calling [`hilbert_index`] +/// For large batches, this is significantly faster than calling [`try_hilbert_index`] /// individually for each point. +/// If the same pre-quantized batch is reused, construct a [`HilbertQuantizedBatch`] +/// once and call [`HilbertQuantizedBatch::indices`] to avoid repeated validation. /// /// # Errors /// @@ -723,7 +1308,7 @@ fn validate_prequantized_coordinates( /// /// ```rust /// use delaunay::prelude::ordering::{ -/// HilbertBitDepth, HilbertError, hilbert_indices_prequantized, hilbert_quantize, +/// HilbertBitDepth, HilbertError, hilbert_indices_prequantized, try_hilbert_quantize, /// }; /// /// let coords = vec![[0.1_f64, 0.2], [0.5, 0.5], [0.9, 0.8]]; @@ -733,7 +1318,7 @@ fn validate_prequantized_coordinates( /// // Quantize once /// let quantized: Vec<[u32; 2]> = coords /// .iter() -/// .map(|c| hilbert_quantize(c, bounds, bits)) +/// .map(|c| try_hilbert_quantize(c, bounds, bits)) /// .collect::>()?; /// /// // Compute all indices in bulk @@ -777,19 +1362,55 @@ pub fn hilbert_indices_prequantized( quantized: &[[u32; D]], bits: HilbertBitDepth, ) -> Result, HilbertError> { - validate_index_params::(bits)?; - - // Handle D == 0 case: zero-dimensional space has only one point, all map to index 0 - if D == 0 { - return Ok(vec![0_u128; quantized.len()]); - } + Ok(HilbertQuantizedBatch::try_new(quantized, bits)?.indices()) +} - validate_prequantized_coordinates(quantized, bits)?; +/// Computes Hilbert indices from a validated pre-quantized batch. +/// +/// This is the infallible companion to [`hilbert_indices_prequantized`]. The +/// [`HilbertQuantizedBatch`] constructor has already checked both the Hilbert +/// index width and every coordinate against the selected bit-depth grid. +/// +/// # Examples +/// +/// ```rust +/// use delaunay::prelude::ordering::{ +/// HilbertBitDepth, HilbertError, HilbertQuantizedBatch, +/// hilbert_indices_for_quantized_batch, +/// }; +/// +/// let quantized = [[0_u32, 0], [3, 3]]; +/// let bits = HilbertBitDepth::try_new(2)?; +/// let batch = HilbertQuantizedBatch::try_new(&quantized, bits)?; +/// +/// let indices = hilbert_indices_for_quantized_batch(batch); +/// assert_eq!(indices.len(), quantized.len()); +/// # Ok::<(), HilbertError>(()) +/// ``` +#[must_use] +pub fn hilbert_indices_for_quantized_batch( + batch: HilbertQuantizedBatch<'_, D>, +) -> Vec { + indices_for_mode(batch.quantized, batch.index_mode) +} - Ok(quantized - .iter() - .map(|q| index_from_quantized(q, bits)) - .collect()) +/// Maps validated quantized coordinates to Hilbert indices for a known index +/// mode. +/// +/// Both [`HilbertQuantizedBatch`] and [`HilbertQuantizedVec`] carry a validated +/// [`HilbertIndexMode`], so neither revalidates coordinates before computing +/// indices. +fn indices_for_mode( + quantized: &[[u32; D]], + index_mode: HilbertIndexMode, +) -> Vec { + match index_mode { + HilbertIndexMode::ZeroDimensional { .. } => vec![0_u128; quantized.len()], + HilbertIndexMode::Positive(index_params) => quantized + .iter() + .map(|q| index_from_quantized(q, index_params)) + .collect(), + } } /// Return the indices that would sort `coords` by Hilbert order. @@ -820,33 +1441,74 @@ pub fn hilbert_indices_prequantized( /// # Examples /// /// ```rust -/// use delaunay::prelude::ordering::{HilbertBitDepth, HilbertError, hilbert_sorted_indices}; +/// use delaunay::prelude::ordering::{HilbertBitDepth, HilbertError, try_hilbert_sorted_indices}; /// /// let coords = vec![[0.9_f64, 0.9], [0.1, 0.1], [0.5, 0.5]]; -/// let order = hilbert_sorted_indices(&coords, (0.0, 1.0), HilbertBitDepth::try_new(8)?)?; +/// let order = try_hilbert_sorted_indices(&coords, (0.0, 1.0), HilbertBitDepth::try_new(8)?)?; /// assert_eq!(order.len(), coords.len()); /// # Ok::<(), HilbertError>(()) /// ``` -pub fn hilbert_sorted_indices( +pub fn try_hilbert_sorted_indices( coords: &[[T; D]], bounds: (T, T), bits: HilbertBitDepth, ) -> Result, HilbertError> { - validate_index_params::(bits)?; + let bounds = parse_hilbert_bounds(bounds)?; + hilbert_sorted_indices_in_range(coords, bounds, bits) +} - if D == 0 { +/// Return the indices that would sort `coords` by Hilbert order using validated bounds. +/// +/// This is equivalent to [`try_hilbert_sorted_indices`], but accepts a +/// [`CoordinateRange`] so callers can carry range validation evidence inward. +/// +/// # Errors +/// +/// Returns [`HilbertError::IndexOverflow`] if `D * bits > 128` (index would not fit in `u128`). +/// Returns [`HilbertError::DimensionTooLarge`] if the dimension `D` exceeds `u32::MAX`. +/// Returns [`HilbertError::QuantizationScaleConversionFailed`] if the quantization +/// grid maximum cannot be represented by the coordinate scalar type. Returns +/// [`HilbertError::NonFiniteBoundsExtent`], [`HilbertError::NonFiniteCoordinate`], +/// [`HilbertError::NonFiniteNormalizedCoordinate`], or +/// [`HilbertError::QuantizedCoordinateConversionFailed`] if coordinate +/// quantization fails. +/// +/// # Examples +/// +/// ```rust +/// use delaunay::prelude::geometry::CoordinateRange; +/// use delaunay::prelude::ordering::{ +/// HilbertBitDepth, HilbertError, hilbert_sorted_indices_in_range, +/// }; +/// +/// # fn main() -> Result<(), HilbertError> { +/// let coords = vec![[0.9_f64, 0.9], [0.1, 0.1], [0.5, 0.5]]; +/// # let Ok(bounds) = CoordinateRange::try_new(0.0_f64, 1.0) else { +/// # return Ok(()); +/// # }; +/// let order = hilbert_sorted_indices_in_range(&coords, bounds, HilbertBitDepth::try_new(8)?)?; +/// assert_eq!(order.len(), coords.len()); +/// # Ok(()) +/// # } +/// ``` +pub fn hilbert_sorted_indices_in_range( + coords: &[[T; D]], + bounds: CoordinateRange, + bits: HilbertBitDepth, +) -> Result, HilbertError> { + let index_mode = HilbertIndexMode::try_new(bits)?; + let HilbertIndexMode::Positive(index_params) = index_mode else { return Ok((0..coords.len()).collect()); - } + }; let (max_val_u32, max_val_t) = quantization_scale::(bits)?; - let bounds = parse_hilbert_bounds(bounds)?; let mut keyed: Vec<((u128, [u32; D]), usize)> = coords .iter() .enumerate() .map(|(i, c)| { let q = quantize_with_scale(c, bounds, bits, max_val_u32, max_val_t)?; - let idx = index_from_quantized(&q, bits); + let idx = index_from_quantized(&q, index_params); Ok(((idx, q), i)) }) .collect::>()?; @@ -867,6 +1529,43 @@ mod tests { HilbertBitDepth::try_new(value).expect("test bit depth must be valid") } + fn positive_index_params(bits: HilbertBitDepth) -> HilbertIndexParams { + let HilbertIndexMode::Positive(params) = + HilbertIndexMode::::try_new(bits).expect("test index parameters must be valid") + else { + panic!("test index parameters must be positive-dimensional"); + }; + params + } + + struct LyingOrder { + values: std::vec::IntoIter, + reported_len: usize, + } + + impl LyingOrder { + fn new(values: Vec, reported_len: usize) -> Self { + Self { + values: values.into_iter(), + reported_len, + } + } + } + + impl Iterator for LyingOrder { + type Item = usize; + + fn next(&mut self) -> Option { + self.values.next() + } + } + + impl ExactSizeIterator for LyingOrder { + fn len(&self) -> usize { + self.reported_len + } + } + /// Asserts that the bulk pre-quantized API matches per-point Hilbert indexing. fn assert_prequantized_matches_hilbert_index( coords: &[[f64; D]], @@ -875,13 +1574,13 @@ mod tests { ) { let quantized: Vec<[u32; D]> = coords .iter() - .map(|c| hilbert_quantize(c, bounds, bits).unwrap()) + .map(|c| try_hilbert_quantize(c, bounds, bits).unwrap()) .collect(); let indices_bulk = hilbert_indices_prequantized(&quantized, bits).expect("valid quantized points"); let indices_individual: Vec = coords .iter() - .map(|c| hilbert_index(c, bounds, bits).unwrap()) + .map(|c| try_hilbert_index(c, bounds, bits).unwrap()) .collect(); assert_eq!(indices_bulk, indices_individual); @@ -918,9 +1617,9 @@ mod tests { #[test] fn test_hilbert_index_2d() { let bits = bit_depth(4); - let origin = hilbert_index(&[0.0_f64, 0.0], (0.0, 1.0), bits).unwrap(); - let corner = hilbert_index(&[1.0_f64, 1.0], (0.0, 1.0), bits).unwrap(); - let center = hilbert_index(&[0.5_f64, 0.5], (0.0, 1.0), bits).unwrap(); + let origin = try_hilbert_index(&[0.0_f64, 0.0], (0.0, 1.0), bits).unwrap(); + let corner = try_hilbert_index(&[1.0_f64, 1.0], (0.0, 1.0), bits).unwrap(); + let center = try_hilbert_index(&[0.5_f64, 0.5], (0.0, 1.0), bits).unwrap(); assert_eq!(origin, 0); assert_ne!(origin, center); @@ -930,8 +1629,8 @@ mod tests { #[test] fn test_hilbert_index_3d() { let bits = bit_depth(8); - let origin = hilbert_index(&[0.0_f64, 0.0, 0.0], (-1.0, 1.0), bits).unwrap(); - let corner = hilbert_index(&[1.0_f64, 1.0, 1.0], (-1.0, 1.0), bits).unwrap(); + let origin = try_hilbert_index(&[0.0_f64, 0.0, 0.0], (-1.0, 1.0), bits).unwrap(); + let corner = try_hilbert_index(&[1.0_f64, 1.0, 1.0], (-1.0, 1.0), bits).unwrap(); assert_ne!(origin, corner); } @@ -996,16 +1695,16 @@ mod tests { let coords: Vec<[f64; 2]> = vec![[0.9, 0.9], [0.1, 0.1], [0.5, 0.5], [0.1, 0.9], [0.9, 0.1]]; let bits = bit_depth(16); - let order = hilbert_sorted_indices(&coords, (0.0, 1.0), bits).unwrap(); + let order = try_hilbert_sorted_indices(&coords, (0.0, 1.0), bits).unwrap(); assert_eq!(order.len(), coords.len()); // Apply the ordering to a parallel payload. let mut payload: Vec = (0..coords.len()).collect(); - hilbert_sort_by_stable(&mut payload, (0.0_f64, 1.0), bits, |&i| coords[i]).unwrap(); + try_hilbert_sort_by_stable(&mut payload, (0.0_f64, 1.0), bits, |&i| coords[i]).unwrap(); // Sorting by stable helper should be deterministic. let mut payload2: Vec = (0..coords.len()).collect(); - hilbert_sort_by_stable(&mut payload2, (0.0_f64, 1.0), bits, |&i| coords[i]).unwrap(); + try_hilbert_sort_by_stable(&mut payload2, (0.0_f64, 1.0), bits, |&i| coords[i]).unwrap(); assert_eq!(payload, order); assert_eq!(payload, payload2); } @@ -1014,11 +1713,11 @@ mod tests { fn test_sort_helpers_accept_stateful_coordinate_closures() { let coords: Vec<[f64; 2]> = vec![[0.9, 0.9], [0.1, 0.1], [0.5, 0.5]]; let bits = bit_depth(8); - let expected_order = hilbert_sorted_indices(&coords, (0.0, 1.0), bits).unwrap(); + let expected_order = try_hilbert_sorted_indices(&coords, (0.0, 1.0), bits).unwrap(); let mut stable_calls = 0_usize; let mut stable_payload: Vec = (0..coords.len()).collect(); - hilbert_sort_by_stable(&mut stable_payload, (0.0_f64, 1.0), bits, |&i| { + try_hilbert_sort_by_stable(&mut stable_payload, (0.0_f64, 1.0), bits, |&i| { stable_calls += 1; coords[i] }) @@ -1028,7 +1727,7 @@ mod tests { let mut unstable_calls = 0_usize; let mut unstable_payload: Vec = (0..coords.len()).collect(); - hilbert_sort_by_unstable(&mut unstable_payload, (0.0_f64, 1.0), bits, |&i| { + try_hilbert_sort_by_unstable(&mut unstable_payload, (0.0_f64, 1.0), bits, |&i| { unstable_calls += 1; coords[i] }) @@ -1042,10 +1741,10 @@ mod tests { let coords: Vec<[f64; 2]> = vec![[0.9, 0.9], [0.1, 0.1], [0.5, 0.5], [0.1, 0.9], [0.9, 0.1]]; let bits = bit_depth(16); - let expected_order = hilbert_sorted_indices(&coords, (0.0, 1.0), bits).unwrap(); + let expected_order = try_hilbert_sorted_indices(&coords, (0.0, 1.0), bits).unwrap(); let mut payload: Vec = (0..coords.len()).collect(); - hilbert_sort_by_unstable(&mut payload, (0.0_f64, 1.0), bits, |&i| coords[i]).unwrap(); + try_hilbert_sort_by_unstable(&mut payload, (0.0_f64, 1.0), bits, |&i| coords[i]).unwrap(); assert_eq!(payload, expected_order); } @@ -1054,18 +1753,53 @@ mod tests { fn test_zero_dim_sort_helpers_noop() { let coords: Vec<[f64; 0]> = vec![[], [], []]; let bits = bit_depth(8); - let order = hilbert_sorted_indices(&coords, (0.0, 1.0), bits).unwrap(); + let order = try_hilbert_sorted_indices(&coords, (0.0, 1.0), bits).unwrap(); assert_eq!(order, vec![0, 1, 2]); let mut stable_payload = vec![3, 2, 1]; - hilbert_sort_by_stable(&mut stable_payload, (0.0_f64, 1.0), bits, |_| []).unwrap(); + try_hilbert_sort_by_stable(&mut stable_payload, (0.0_f64, 1.0), bits, |_| []).unwrap(); assert_eq!(stable_payload, vec![3, 2, 1]); let mut unstable_payload = vec![3, 2, 1]; - hilbert_sort_by_unstable(&mut unstable_payload, (0.0_f64, 1.0), bits, |_| []).unwrap(); + try_hilbert_sort_by_unstable(&mut unstable_payload, (0.0_f64, 1.0), bits, |_| []).unwrap(); assert_eq!(unstable_payload, vec![3, 2, 1]); } + #[test] + fn test_zero_dim_raw_bound_apis_reject_invalid_bounds() { + let bits = bit_depth(8); + let coords = [0.0_f64; 0]; + let coordinate_batch = vec![coords; 3]; + let mut payload = vec![3, 2, 1]; + + assert_matches!( + try_hilbert_quantize(&coords, (f64::NAN, 1.0), bits), + Err(HilbertError::NonFiniteBounds { + lower_bound_finite: false, + upper_bound_finite: true + }) + ); + assert_matches!( + try_hilbert_index(&coords, (1.0, 1.0), bits), + Err(HilbertError::NonIncreasingBounds { + ordering: CoordinateRangeOrdering::Equal + }) + ); + assert_matches!( + try_hilbert_sorted_indices(&coordinate_batch, (1.0, 0.0), bits), + Err(HilbertError::NonIncreasingBounds { + ordering: CoordinateRangeOrdering::Decreasing + }) + ); + assert_matches!( + try_hilbert_sort_by_stable(&mut payload, (f64::NEG_INFINITY, 1.0), bits, |_| coords), + Err(HilbertError::NonFiniteBounds { + lower_bound_finite: false, + upper_bound_finite: true + }) + ); + } + #[test] fn test_scaled_quantize_reports_conversion_error() { let result = quantize_with_scale( @@ -1086,9 +1820,196 @@ mod tests { ); } + macro_rules! gen_in_range_quantization_tests { + ($dim:literal, $points:expr, $sample:expr) => { + pastey::paste! { + #[test] + fn []() { + let bits = bit_depth(8); + let range = CoordinateRange::try_new(-2.0_f64, 3.0).unwrap(); + let coords: [f64; $dim] = $sample; + + let parsed = try_hilbert_quantize(&coords, range.bounds(), bits).unwrap(); + let prevalidated = hilbert_quantize_in_range(&coords, range, bits).unwrap(); + + assert_eq!(prevalidated, parsed); + } + + #[test] + fn []() { + let bits = bit_depth(8); + let range = CoordinateRange::try_new(-2.0_f64, 3.0).unwrap(); + let points: [[f64; $dim]; 4] = $points; + let quantized: Vec<[u32; $dim]> = points + .iter() + .map(|point| hilbert_quantize_in_range(point, range, bits).unwrap()) + .collect(); + let batch = HilbertQuantizedBatch::try_new(&quantized, bits).unwrap(); + + assert_eq!(batch.coordinates(), quantized.as_slice()); + assert_eq!(batch.bits(), bits); + + let checked = hilbert_indices_prequantized(&quantized, bits).unwrap(); + assert_eq!(batch.indices(), checked); + assert_eq!(hilbert_indices_for_quantized_batch(batch), checked); + } + + #[test] + fn []() { + let bits = bit_depth(8); + let bounds = CoordinateRange::try_new(-2.0_f64, 3.0).unwrap(); + let points: [[f64; $dim]; 4] = $points; + + let two_step: Vec<[u32; $dim]> = points + .iter() + .map(|point| hilbert_quantize_in_range(point, bounds, bits).unwrap()) + .collect(); + let two_step_indices = hilbert_indices_prequantized(&two_step, bits).unwrap(); + + let batch = hilbert_quantize_batch_in_range(&points, bounds, bits, |point| *point) + .unwrap(); + assert_eq!(batch.coordinates(), two_step.as_slice()); + assert_eq!(batch.bits(), bits); + assert_eq!(batch.len(), points.len()); + assert!(!batch.is_empty()); + + let (indices, quantized) = batch.into_indices_and_coordinates(); + assert_eq!(quantized, two_step); + assert_eq!(indices, two_step_indices); + } + } + }; + } + + gen_in_range_quantization_tests!( + 2, + [[-2.0_f64, -1.0], [-1.5, 0.25], [0.1, -0.7], [3.0, 3.0]], + [0.25_f64, 0.75] + ); + gen_in_range_quantization_tests!( + 3, + [ + [-2.0_f64, -1.0, 0.0], + [-1.5, 0.25, 1.75], + [0.1, -0.7, 2.2], + [3.0, 3.0, -2.0] + ], + [0.25_f64, 0.75, -1.0] + ); + gen_in_range_quantization_tests!( + 4, + [ + [-2.0_f64, -1.0, 0.0, 1.0], + [-1.5, 0.25, 1.75, 2.5], + [0.1, -0.7, 2.2, -1.8], + [3.0, 3.0, -2.0, -2.0], + ], + [0.25_f64, 0.75, -1.0, 2.5] + ); + gen_in_range_quantization_tests!( + 5, + [ + [-2.0_f64, -1.0, 0.0, 1.0, 2.0], + [-1.5, 0.25, 1.75, 2.5, -0.5], + [0.1, -0.7, 2.2, -1.8, 1.4], + [3.0, 3.0, -2.0, -2.0, 0.5], + ], + [0.25_f64, 0.75, -1.0, 2.5, 0.0] + ); + + #[test] + fn test_quantized_batch_rejects_out_of_range_coordinate() { + let bits = bit_depth(2); + let quantized = [[0_u32, 0], [4, 1]]; + let result = HilbertQuantizedBatch::try_new(&quantized, bits); + + assert_matches!( + result, + Err(HilbertError::PrequantizedCoordinateOutOfRange { + bits: 2, + max_grid_value: 3, + point_index: 1, + coordinate_index: 0, + coordinate: 4 + }) + ); + } + + #[test] + fn test_quantized_batch_rejects_index_overflow() { + let quantized = [[1_u32, 2, 3, 4, 5]]; + let result = HilbertQuantizedBatch::try_new(&quantized, bit_depth(26)); + + assert_matches!( + result, + Err(HilbertError::IndexOverflow { + dimension: 5, + bits: 26, + total_bits: 130 + }) + ); + } + + #[test] + fn test_quantized_batch_handles_zero_dimension() { + let bits = bit_depth(8); + let quantized = [[], [], []]; + let batch = HilbertQuantizedBatch::try_new(&quantized, bits).unwrap(); + + assert_eq!(batch.coordinates(), quantized.as_slice()); + assert_eq!(batch.bits(), bits); + assert_eq!(batch.indices(), vec![0_u128, 0_u128, 0_u128]); + assert_eq!( + hilbert_indices_for_quantized_batch(batch), + vec![0_u128, 0_u128, 0_u128] + ); + } + + #[test] + fn test_quantize_batch_in_range_handles_zero_dimension() { + let bits = bit_depth(8); + let bounds = CoordinateRange::try_new(0.0_f64, 1.0).unwrap(); + let items = [(), (), ()]; + + let batch = + hilbert_quantize_batch_in_range(&items, bounds, bits, |()| [0.0_f64; 0]).unwrap(); + assert_eq!(batch.len(), 3); + assert!(!batch.is_empty()); + assert_eq!(batch.indices(), vec![0_u128, 0_u128, 0_u128]); + assert_eq!(batch.into_coordinates(), vec![[0_u32; 0]; 3]); + } + + #[test] + fn test_quantize_in_range_handles_zero_dimension() { + let bits = bit_depth(8); + let bounds = CoordinateRange::try_new(0.0_f64, 1.0).unwrap(); + let coords = [0.0_f64; 0]; + + assert_eq!( + hilbert_quantize_in_range(&coords, bounds, bits).unwrap(), + [0_u32; 0] + ); + } + + #[test] + fn test_quantize_batch_in_range_rejects_index_overflow() { + let bits = bit_depth(26); + let bounds = CoordinateRange::try_new(0.0_f64, 1.0).unwrap(); + let items = [[0.0_f64; 5]]; + + assert_matches!( + hilbert_quantize_batch_in_range(&items, bounds, bits, |p| *p), + Err(HilbertError::IndexOverflow { + dimension: 5, + bits: 26, + total_bits: 130 + }) + ); + } + #[test] fn test_quantize_rejects_nonfinite_bounds() { - let result = hilbert_quantize(&[0.5_f64], (f64::NAN, 1.0), bit_depth(8)); + let result = try_hilbert_quantize(&[0.5_f64], (f64::NAN, 1.0), bit_depth(8)); assert_eq!( result, @@ -1098,7 +2019,8 @@ mod tests { }) ); - let both_non_finite = hilbert_quantize(&[0.5_f64], (f64::NAN, f64::INFINITY), bit_depth(8)); + let both_non_finite = + try_hilbert_quantize(&[0.5_f64], (f64::NAN, f64::INFINITY), bit_depth(8)); assert_eq!( both_non_finite, Err(HilbertError::NonFiniteBounds { @@ -1110,14 +2032,14 @@ mod tests { #[test] fn test_quantize_rejects_nonfinite_extent() { - let result = hilbert_quantize(&[0.0_f64], (-f64::MAX, f64::MAX), bit_depth(8)); + let result = try_hilbert_quantize(&[0.0_f64], (-f64::MAX, f64::MAX), bit_depth(8)); assert_eq!(result, Err(HilbertError::NonFiniteBoundsExtent {})); } #[test] fn test_quantize_rejects_nonfinite_coordinate() { - let result = hilbert_quantize(&[0.25_f64, f64::INFINITY], (0.0, 1.0), bit_depth(8)); + let result = try_hilbert_quantize(&[0.25_f64, f64::INFINITY], (0.0, 1.0), bit_depth(8)); assert_eq!( result, @@ -1129,7 +2051,8 @@ mod tests { #[test] fn test_quantize_rejects_nonfinite_normalized() { - let result = hilbert_quantize(&[f64::MAX], (-f64::MAX / 2.0, f64::MAX / 2.0), bit_depth(8)); + let result = + try_hilbert_quantize(&[f64::MAX], (-f64::MAX / 2.0, f64::MAX / 2.0), bit_depth(8)); assert_eq!( result, @@ -1144,7 +2067,8 @@ mod tests { let coords = [[0.5_f64], [f64::NAN], [0.25]]; let mut payload = vec![0_usize, 1, 2]; - let result = hilbert_sort_by_stable(&mut payload, (0.0, 1.0), bit_depth(8), |&i| coords[i]); + let result = + try_hilbert_sort_by_stable(&mut payload, (0.0, 1.0), bit_depth(8), |&i| coords[i]); assert_eq!( result, @@ -1155,9 +2079,85 @@ mod tests { assert_eq!(payload, vec![0, 1, 2]); } + #[test] + fn test_apply_order_rejects_length_mismatch_without_reordering() { + let mut payload = vec![10, 20, 30]; + let result = apply_order(&mut payload, [0_usize, 1].into_iter()); + + assert_eq!(payload, vec![10, 20, 30]); + assert_matches!( + result, + Err(HilbertError::InvalidSortPermutationLength { + item_count: 3, + permutation_count: 2 + }) + ); + } + + #[test] + fn test_apply_order_rejects_out_of_range_index_without_reordering() { + let mut payload = vec![10, 20, 30]; + let result = apply_order(&mut payload, [0_usize, 3, 1].into_iter()); + + assert_eq!(payload, vec![10, 20, 30]); + assert_matches!( + result, + Err(HilbertError::InvalidSortPermutationIndex { + permutation_index: 1, + item_index: 3, + item_count: 3 + }) + ); + } + + #[test] + fn test_apply_order_rejects_duplicate_index_without_reordering() { + let mut payload = vec![10, 20, 30]; + let result = apply_order(&mut payload, [0_usize, 1, 1].into_iter()); + + assert_eq!(payload, vec![10, 20, 30]); + assert_matches!( + result, + Err(HilbertError::InvalidSortPermutationDuplicate { + permutation_index: 2, + item_index: 1 + }) + ); + } + + #[test] + fn test_apply_order_rejects_short_iterator_length_lie_without_reordering() { + let mut payload = vec![10, 20, 30]; + let result = apply_order(&mut payload, LyingOrder::new(vec![0, 1], 3)); + + assert_eq!(payload, vec![10, 20, 30]); + assert_matches!( + result, + Err(HilbertError::InvalidSortPermutationLength { + item_count: 3, + permutation_count: 2 + }) + ); + } + + #[test] + fn test_apply_order_rejects_long_iterator_length_lie_without_reordering() { + let mut payload = vec![10, 20, 30]; + let result = apply_order(&mut payload, LyingOrder::new(vec![0, 1, 2, 0], 3)); + + assert_eq!(payload, vec![10, 20, 30]); + assert_matches!( + result, + Err(HilbertError::InvalidSortPermutationLength { + item_count: 3, + permutation_count: 4 + }) + ); + } + #[test] fn test_quantize_clamps_f64_endpoint() { - let q = hilbert_quantize(&[1.0], (0.0, 1.0), bit_depth(31)).unwrap(); + let q = try_hilbert_quantize(&[1.0], (0.0, 1.0), bit_depth(31)).unwrap(); assert_eq!(q, [(1_u32 << 31) - 1]); } @@ -1170,10 +2170,11 @@ mod tests { let n: u32 = 1_u32 << bits; let mut points: Vec<([u32; 2], u128)> = Vec::with_capacity((n * n) as usize); + let params = positive_index_params::<2>(bit_depth(bits)); for x in 0..n { for y in 0..n { let q = [x, y]; - let idx = index_from_quantized(&q, bit_depth(bits)); + let idx = index_from_quantized(&q, params); points.push((q, idx)); } } @@ -1204,12 +2205,13 @@ mod tests { let n: u32 = 1_u32 << bits; let mut points: Vec<([u32; 4], u128)> = Vec::with_capacity((n * n * n * n) as usize); + let params = positive_index_params::<4>(bit_depth(bits)); for x in 0..n { for y in 0..n { for z in 0..n { for w in 0..n { let q = [x, y, z, w]; - let idx = index_from_quantized(&q, bit_depth(bits)); + let idx = index_from_quantized(&q, params); points.push((q, idx)); } } @@ -1239,22 +2241,22 @@ mod tests { #[test] fn test_point_coords_work_with_hilbert() { let p: Point = Point::new([0.25, 0.75]); - let idx = hilbert_index(p.coords(), (0.0, 1.0), bit_depth(16)).unwrap(); + let idx = try_hilbert_index(p.coords(), (0.0, 1.0), bit_depth(16)).unwrap(); assert!(idx > 0); } #[test] fn test_hilbert_bits_boundaries() { let coarsest_bits = bit_depth(1); - let origin = hilbert_index(&[0.0_f64, 0.0], (0.0, 1.0), coarsest_bits).unwrap(); - let corner = hilbert_index(&[1.0_f64, 1.0], (0.0, 1.0), coarsest_bits).unwrap(); + let origin = try_hilbert_index(&[0.0_f64, 0.0], (0.0, 1.0), coarsest_bits).unwrap(); + let corner = try_hilbert_index(&[1.0_f64, 1.0], (0.0, 1.0), coarsest_bits).unwrap(); tracing::debug!(origin, corner, "bits=1 boundaries"); assert_eq!(origin, 0, "bits=1 origin should map to 0"); assert_ne!(origin, corner, "bits=1 should distinguish corners"); let finest_bits = bit_depth(31); - let origin_31 = hilbert_index(&[0.0_f64, 0.0], (0.0, 1.0), finest_bits).unwrap(); - let corner_31 = hilbert_index(&[1.0_f64, 1.0], (0.0, 1.0), finest_bits).unwrap(); + let origin_31 = try_hilbert_index(&[0.0_f64, 0.0], (0.0, 1.0), finest_bits).unwrap(); + let corner_31 = try_hilbert_index(&[1.0_f64, 1.0], (0.0, 1.0), finest_bits).unwrap(); tracing::debug!(origin_31, corner_31, "bits=31 boundaries"); assert_eq!(origin_31, 0, "bits=31 origin should map to 0"); assert_ne!(origin_31, corner_31, "bits=31 should distinguish corners"); @@ -1264,10 +2266,10 @@ mod tests { fn test_hilbert_index_1d_monotonic() { let bounds = (0.0_f64, 1.0_f64); let bits = bit_depth(8); - let a = hilbert_index(&[0.0_f64], bounds, bits).unwrap(); - let b = hilbert_index(&[0.25_f64], bounds, bits).unwrap(); - let c = hilbert_index(&[0.5_f64], bounds, bits).unwrap(); - let d = hilbert_index(&[1.0_f64], bounds, bits).unwrap(); + let a = try_hilbert_index(&[0.0_f64], bounds, bits).unwrap(); + let b = try_hilbert_index(&[0.25_f64], bounds, bits).unwrap(); + let c = try_hilbert_index(&[0.5_f64], bounds, bits).unwrap(); + let d = try_hilbert_index(&[1.0_f64], bounds, bits).unwrap(); tracing::debug!(a, b, c, d, "1d indices"); assert!( a < b && b < c && c < d, @@ -1282,13 +2284,13 @@ mod tests { let bits = bit_depth(8); assert_matches!( - hilbert_quantize(&coords, bounds, bits), + try_hilbert_quantize(&coords, bounds, bits), Err(HilbertError::NonIncreasingBounds { ordering: CoordinateRangeOrdering::Equal }) ); assert_matches!( - hilbert_index(&coords, bounds, bits), + try_hilbert_index(&coords, bounds, bits), Err(HilbertError::NonIncreasingBounds { ordering: CoordinateRangeOrdering::Equal }) @@ -1302,7 +2304,7 @@ mod tests { let bits = bit_depth(8); assert_matches!( - hilbert_quantize(&coords, bounds, bits), + try_hilbert_quantize(&coords, bounds, bits), Err(HilbertError::NonIncreasingBounds { ordering: CoordinateRangeOrdering::Decreasing }) @@ -1314,7 +2316,7 @@ mod tests { let bounds = (0.0_f64, 1.0_f64); let bits = bit_depth(4); let coords = [-1.0_f64, 2.0_f64]; - let q = hilbert_quantize(&coords, bounds, bits).unwrap(); + let q = try_hilbert_quantize(&coords, bounds, bits).unwrap(); let max_val = (1_u32 << bits.get()) - 1; tracing::debug!(?q, max_val, "clamp quantize"); assert_eq!( @@ -1323,8 +2325,8 @@ mod tests { "out-of-range coords should clamp to bounds" ); - let idx = hilbert_index(&coords, bounds, bits).unwrap(); - let idx_clamped = hilbert_index(&[0.0_f64, 1.0_f64], bounds, bits).unwrap(); + let idx = try_hilbert_index(&coords, bounds, bits).unwrap(); + let idx_clamped = try_hilbert_index(&[0.0_f64, 1.0_f64], bounds, bits).unwrap(); tracing::debug!(idx, idx_clamped, "clamp index"); assert_eq!( idx, idx_clamped, @@ -1347,7 +2349,7 @@ mod tests { // Quantize all coordinates let quantized: Vec<[u32; 3]> = coords .iter() - .map(|c| hilbert_quantize(c, bounds, bits).unwrap()) + .map(|c| try_hilbert_quantize(c, bounds, bits).unwrap()) .collect(); // Compute indices via bulk API @@ -1357,7 +2359,7 @@ mod tests { // Compute indices individually let indices_individual: Vec = coords .iter() - .map(|c| hilbert_index(c, bounds, bits).unwrap()) + .map(|c| try_hilbert_index(c, bounds, bits).unwrap()) .collect(); assert_eq!(indices_bulk.len(), coords.len()); @@ -1444,7 +2446,7 @@ mod tests { ]; for (coord, expected_cell) in test_cases { - let q = hilbert_quantize(&[coord], bounds, bits).unwrap(); + let q = try_hilbert_quantize(&[coord], bounds, bits).unwrap(); assert_eq!( q[0], expected_cell, "coordinate {coord} should quantize to cell {expected_cell}, got {}", @@ -1463,7 +2465,7 @@ mod tests { let mut cell_counts = [0_usize; 4]; for i in 0..samples { let coord = f64::from(i) / f64::from(samples); - let q = hilbert_quantize(&[coord], bounds, bits).unwrap(); + let q = try_hilbert_quantize(&[coord], bounds, bits).unwrap(); cell_counts[q[0] as usize] += 1; } diff --git a/src/core/validation.rs b/src/core/validation.rs index 788be638..b200b403 100644 --- a/src/core/validation.rs +++ b/src/core/validation.rs @@ -1561,8 +1561,9 @@ mod tests { use crate::core::tds::{GeometricError, NeighborValidationError, Tds}; use crate::core::vertex::Vertex; use crate::core::vertex::VertexBuilder; + use crate::geometry::coordinate_range::CoordinateRange; use crate::geometry::kernel::FastKernel; - use crate::geometry::util::generate_random_points_seeded; + use crate::geometry::util::generate_random_points_in_range_seeded; use crate::repair::DelaunayRepairPolicy; use crate::triangulation::DelaunayTriangulation; use crate::validation::DelaunayTriangulationValidationError; @@ -3118,7 +3119,8 @@ mod tests { #[test] fn pl_manifold_insertion_keeps_valid_topology_with_explicit_only_validation() { - let points = generate_random_points_seeded::(25, (-100.0, 100.0), 123).unwrap(); + let bounds = CoordinateRange::try_new(-100.0_f64, 100.0).unwrap(); + let points = generate_random_points_in_range_seeded::(25, bounds, 123); let mut dt: DelaunayTriangulation<_, (), (), 3> = DelaunayTriangulation::empty_with_topology_guarantee(TopologyGuarantee::PLManifold); diff --git a/src/delaunay/construction.rs b/src/delaunay/construction.rs index f4e8bc52..3fb0f44f 100644 --- a/src/delaunay/construction.rs +++ b/src/delaunay/construction.rs @@ -70,12 +70,13 @@ use crate::core::tds::{TdsConstructionError, TdsError}; use crate::core::traits::data_type::DataType; use crate::core::triangulation::Triangulation; use crate::core::util::{ - HilbertBitDepth, coords_equal_exact, coords_within_epsilon, hilbert_indices_prequantized, - hilbert_quantize, stable_hash_u64_slice, + HilbertBitDepth, coords_equal_exact, coords_within_epsilon, hilbert_quantize_batch_in_range, + stable_hash_u64_slice, }; use crate::core::validation::{TopologyGuarantee, TriangulationValidationError, ValidationPolicy}; use crate::core::vertex::Vertex; use crate::diagnostics::{BatchLocalRepairTrigger, ConstructionTelemetry, LocalRepairSample}; +use crate::geometry::coordinate_range::CoordinateRange; use crate::geometry::kernel::{AdaptiveKernel, Kernel}; use crate::geometry::point::Point; use crate::geometry::traits::coordinate::{Coordinate, CoordinateScalar, CoordinateValues}; @@ -1838,7 +1839,7 @@ fn hilbert_bits_per_coord() -> Option { return None; }; - // `hilbert_index` encodes D coordinates with `bits` bits each into a `u128`. + // Hilbert indexing encodes D coordinates with `bits` bits each into a `u128`. // Use as many bits as possible (up to the `hilbert` module's `bits <= 31` bound). let bits_per_coord = (128_u32 / d_u32).min(31); HilbertBitDepth::try_new(bits_per_coord).ok() @@ -1861,7 +1862,7 @@ where .collect() } -/// Sort key for Hilbert ordering: `(hilbert_index, quantized_coords, vertex, input_index)`. +/// Sort key for Hilbert ordering: `(Hilbert index, quantized coords, vertex, input index)`. type HilbertSortKey = (u128, [u32; D], Vertex, usize); /// Orders vertices along a Hilbert curve to improve insertion locality while @@ -1911,39 +1912,34 @@ where return order_vertices_lexicographic(vertices); }; - let bounds = (min_t, max_t); - - // Quantize all coordinates. - let quantized: Result, ()> = vertices - .iter() - .map(|vertex| { - hilbert_quantize(vertex.point().coords(), bounds, bits_per_coord).map_err(|_| ()) - }) - .collect(); - - let Ok(quantized) = quantized else { - // On quantization error, fall back to true lexicographic ordering of original coordinates + let Ok(bounds) = CoordinateRange::try_new(min_t, max_t) else { return order_vertices_lexicographic(vertices); }; - // Compute all indices in bulk. - let Ok(indices) = hilbert_indices_prequantized(&quantized, bits_per_coord) else { - // On bulk index computation error, fall back to true lexicographic ordering + // Quantize all coordinates and compute Hilbert indices in one + // proof-carrying pass. `hilbert_quantize_batch_in_range` validates the + // index width and quantization scale once, clamps every coordinate into + // the bit-depth grid, and returns a batch whose indices are infallible, so + // no per-coordinate range rescan is repeated here. + let Ok(batch) = hilbert_quantize_batch_in_range(&vertices, bounds, bits_per_coord, |vertex| { + *vertex.point().coords() + }) else { + // On quantization error, fall back to true lexicographic ordering. return order_vertices_lexicographic(vertices); }; + + let (indices, quantized) = batch.into_indices_and_coordinates(); + if vertices.len() != quantized.len() || vertices.len() != indices.len() { + return order_vertices_lexicographic(vertices); + } + // Pair indices with vertices, quantized coords, and input indices. let mut keyed: Vec> = vertices .into_iter() + .zip(quantized) + .zip(indices) .enumerate() - .map(|(input_index, vertex)| { - let idx = indices - .get(input_index) - .copied() - // Fallback to input index directly as u128 (no u32 truncation) - .unwrap_or(input_index as u128); - let q = quantized[input_index]; - (idx, q, vertex, input_index) - }) + .map(|(input_index, ((vertex, q), idx))| (idx, q, vertex, input_index)) .collect(); keyed.sort_by( diff --git a/src/geometry/util/point_generation.rs b/src/geometry/util/point_generation.rs index 83dfe29c..6ae9ea87 100644 --- a/src/geometry/util/point_generation.rs +++ b/src/geometry/util/point_generation.rs @@ -123,6 +123,23 @@ pub enum RandomPointGenerationError { value: InvalidCoordinateValue, }, + /// Ball rejection sampling could not produce the requested point count. + #[error( + "Could not generate {requested_points} ball points in dimension {dimension} with radius {radius:?} after {attempts} attempts; generated {generated_points}" + )] + BallSamplingFailed { + /// Number of points requested. + requested_points: usize, + /// Number of points generated before attempts were exhausted. + generated_points: usize, + /// Ball dimension. + dimension: usize, + /// Ball radius used by the sampler. + radius: T, + /// Number of candidate samples attempted. + attempts: usize, + }, + /// Failed to convert a discrete count, index, or scalar into a numeric target type. #[error("Failed to convert {value} to {target_type}: {source}")] CoordinateConversionFailed { @@ -209,6 +226,9 @@ const MAX_GRID_BYTES_SAFETY_CAP_DEFAULT: usize = 4_294_967_296; // 4 GiB /// Base number of Poisson disk sampling attempts per requested point. const POISSON_ATTEMPTS_PER_POINT: usize = 30; +/// Base number of ball rejection-sampling attempts per requested point. +const BALL_ATTEMPTS_PER_POINT: usize = 1_024; + /// Get the maximum bytes allowed for grid allocation. /// /// Reads the `MAX_GRID_BYTES_SAFETY_CAP` environment variable if set, @@ -391,6 +411,13 @@ const fn poisson_max_attempts(n_points: usize, dimension: usize) -> usize { .saturating_mul(poisson_dimension_attempt_scaling(dimension)) } +/// Computes the ball rejection-sampling attempt budget without panicking on large inputs. +const fn ball_max_attempts(n_points: usize, dimension: usize) -> usize { + n_points + .saturating_mul(BALL_ATTEMPTS_PER_POINT) + .saturating_mul(poisson_dimension_attempt_scaling(dimension)) +} + /// Generate random points in D-dimensional space with uniform distribution. /// /// This function provides a flexible way to generate random points for testing, @@ -416,29 +443,29 @@ const fn poisson_max_attempts(n_points: usize, dimension: usize) -> usize { /// /// ``` /// use delaunay::prelude::generators::{ -/// RandomPointGenerationError, generate_random_points, +/// RandomPointGenerationError, try_generate_random_points, /// }; /// /// # fn main() -> Result<(), RandomPointGenerationError> { /// // Generate 100 random 2D points with coordinates in [-10.0, 10.0] -/// let points_2d = generate_random_points::(100, (-10.0, 10.0))?; +/// let points_2d = try_generate_random_points::(100, (-10.0, 10.0))?; /// assert_eq!(points_2d.len(), 100); /// /// // Generate 3D points with coordinates in [0.0, 1.0] (unit cube) -/// let points_3d = generate_random_points::(50, (0.0, 1.0))?; +/// let points_3d = try_generate_random_points::(50, (0.0, 1.0))?; /// assert_eq!(points_3d.len(), 50); /// /// // Generate 4D points centered around origin -/// let points_4d = generate_random_points::(25, (-1.0, 1.0))?; +/// let points_4d = try_generate_random_points::(25, (-1.0, 1.0))?; /// assert_eq!(points_4d.len(), 25); /// /// // Error handling -/// let result = generate_random_points::(100, (10.0, -10.0)); +/// let result = try_generate_random_points::(100, (10.0, -10.0)); /// assert!(result.is_err()); // Invalid range /// # Ok(()) /// # } /// ``` -pub fn generate_random_points( +pub fn try_generate_random_points( n_points: usize, range: (T, T), ) -> Result>, RandomPointGenerationError> { @@ -447,7 +474,7 @@ pub fn generate_random_points Result<(), RandomPointGenerationError> { /// // Generate reproducible random points -/// let points1 = generate_random_points_seeded::(100, (-5.0, 5.0), 42)?; -/// let points2 = generate_random_points_seeded::(100, (-5.0, 5.0), 42)?; +/// let points1 = try_generate_random_points_seeded::(100, (-5.0, 5.0), 42)?; +/// let points2 = try_generate_random_points_seeded::(100, (-5.0, 5.0), 42)?; /// assert_eq!(points1, points2); // Same seed produces identical results /// /// // Different seeds produce different results -/// let points3 = generate_random_points_seeded::(100, (-5.0, 5.0), 123)?; +/// let points3 = try_generate_random_points_seeded::(100, (-5.0, 5.0), 123)?; /// assert_ne!(points1, points3); /// /// // Common ranges - unit cube [0,1] -/// let unit_points = generate_random_points_seeded::(50, (0.0, 1.0), 42)?; +/// let unit_points = try_generate_random_points_seeded::(50, (0.0, 1.0), 42)?; /// /// // Centered around origin [-1,1] -/// let centered_points = generate_random_points_seeded::(50, (-1.0, 1.0), 42)?; +/// let centered_points = try_generate_random_points_seeded::(50, (-1.0, 1.0), 42)?; /// # Ok(()) /// # } /// ``` -pub fn generate_random_points_seeded( +pub fn try_generate_random_points_seeded( n_points: usize, range: (T, T), seed: u64, @@ -551,7 +578,7 @@ pub fn generate_random_points_seeded( radius: T, rng: &mut R, ) -> Result>, RandomPointGenerationError> +where + T: CoordinateScalar + SampleUniform, + R: rand::Rng + ?Sized, +{ + generate_random_points_in_ball_with_rng_and_budget( + n_points, + radius, + rng, + ball_max_attempts(n_points, D), + ) +} + +/// Generates ball samples through an injected RNG with an explicit attempt budget. +fn generate_random_points_in_ball_with_rng_and_budget( + n_points: usize, + radius: T, + rng: &mut R, + max_attempts: usize, +) -> Result>, RandomPointGenerationError> where T: CoordinateScalar + SampleUniform, R: rand::Rng + ?Sized, @@ -695,8 +741,10 @@ where } let mut points = Vec::with_capacity(n_points); + let mut attempts = 0; - while points.len() < n_points { + while points.len() < n_points && attempts < max_attempts { + attempts += 1; let coords = [T::zero(); D].map(|_| rng.random_range(bounds.min()..bounds.max())); let norm_sq = coords.iter().fold(T::zero(), |acc, &c| acc + c * c); if norm_sq <= radius_sq { @@ -704,6 +752,16 @@ where } } + if points.len() < n_points { + return Err(RandomPointGenerationError::BallSamplingFailed { + requested_points: n_points, + generated_points: points.len(), + dimension: D, + radius, + attempts, + }); + } + Ok(points) } @@ -726,9 +784,11 @@ where /// # Errors /// /// Returns [`RandomPointGenerationError::InvalidBallRadius`] if `radius` is -/// non-finite or non-positive, or +/// non-finite or non-positive, /// [`RandomPointGenerationError::InvalidBallRadiusSquared`] if squaring a -/// finite radius overflows to a non-finite value. +/// finite radius overflows to a non-finite value, or +/// [`RandomPointGenerationError::BallSamplingFailed`] if rejection sampling +/// exhausts its attempt budget before producing `n_points`. /// /// # Examples /// @@ -765,9 +825,11 @@ pub fn generate_random_points_in_ball( /// /// ``` /// use delaunay::prelude::generators::{ -/// RandomPointGenerationError, generate_poisson_points, +/// RandomPointGenerationError, try_generate_poisson_points, /// }; /// /// # fn main() -> Result<(), RandomPointGenerationError> { /// // Generate ~100 2D points with minimum distance 0.1 in unit square -/// let poisson_2d = generate_poisson_points::(100, (0.0, 1.0), 0.1, 42)?; +/// let poisson_2d = try_generate_poisson_points::(100, (0.0, 1.0), 0.1, 42)?; /// // Actual count may be less than 100 due to spacing constraints /// /// // Generate 3D points in a cube -/// let poisson_3d = generate_poisson_points::(50, (-1.0, 1.0), 0.2, 123)?; +/// let poisson_3d = try_generate_poisson_points::(50, (-1.0, 1.0), 0.2, 123)?; /// # Ok(()) /// # } /// ``` -pub fn generate_poisson_points( +pub fn try_generate_poisson_points( n_points: usize, bounds: (T, T), min_distance: T, @@ -1219,6 +1281,18 @@ mod tests { assert!(display.contains("Invalid squared ball radius")); assert!(display.contains("inf")); + let ball_error = RandomPointGenerationError::BallSamplingFailed { + requested_points: 4, + generated_points: 1, + dimension: 12, + radius: 2.0, + attempts: 8_192, + }; + let display = format!("{ball_error}"); + assert!(display.contains("Could not generate 4 ball points")); + assert!(display.contains("dimension 12")); + assert!(display.contains("generated 1")); + let generated_grid_coordinate_error: RandomPointGenerationError = RandomPointGenerationError::InvalidGeneratedGridCoordinate { axis: 2, @@ -1255,20 +1329,20 @@ mod tests { #[test] fn test_generate_random_points_rejects_nonfinite_tuple_bounds() { assert_matches!( - generate_random_points::(4, (f64::NAN, 1.0)), + try_generate_random_points::(4, (f64::NAN, 1.0)), Err(RandomPointGenerationError::InvalidCoordinateRange { source: CoordinateRangeError::NonFiniteBound { bound, value } }) if bound == CoordinateRangeBound::Minimum && value == InvalidCoordinateValue::Nan ); assert_matches!( - generate_random_points_seeded::(4, (0.0, f64::INFINITY), 42), + try_generate_random_points_seeded::(4, (0.0, f64::INFINITY), 42), Err(RandomPointGenerationError::InvalidCoordinateRange { source: CoordinateRangeError::NonFiniteBound { bound, value } }) if bound == CoordinateRangeBound::Maximum && value == InvalidCoordinateValue::PositiveInfinity ); assert_matches!( - generate_poisson_points::(4, (f64::NEG_INFINITY, 1.0), 0.1, 42), + try_generate_poisson_points::(4, (f64::NEG_INFINITY, 1.0), 0.1, 42), Err(RandomPointGenerationError::InvalidCoordinateRange { source: CoordinateRangeError::NonFiniteBound { bound, value } }) if bound == CoordinateRangeBound::Minimum @@ -1399,10 +1473,19 @@ mod tests { assert_eq!(poisson_max_attempts(usize::MAX, 7), usize::MAX); } + #[test] + fn test_ball_attempt_budget_saturates_for_large_point_counts() { + assert_eq!(ball_max_attempts(10, 2), 10_240); + assert_eq!(ball_max_attempts(10, 4), 20_480); + assert_eq!(ball_max_attempts(10, 6), 40_960); + assert_eq!(ball_max_attempts(10, 7), 81_920); + assert_eq!(ball_max_attempts(usize::MAX, 7), usize::MAX); + } + #[test] fn test_generate_random_points_2d() { // Test 2D random point generation - let points = generate_random_points::(100, (-10.0, 10.0)).unwrap(); + let points = try_generate_random_points::(100, (-10.0, 10.0)).unwrap(); assert_eq!(points.len(), 100); @@ -1417,7 +1500,7 @@ mod tests { #[test] fn test_generate_random_points_3d() { // Test 3D random point generation - let points = generate_random_points::(75, (0.0, 5.0)).unwrap(); + let points = try_generate_random_points::(75, (0.0, 5.0)).unwrap(); assert_eq!(points.len(), 75); @@ -1432,7 +1515,7 @@ mod tests { #[test] fn test_generate_random_points_4d() { // Test 4D random point generation - let points = generate_random_points::(50, (-2.0, 2.0)).unwrap(); + let points = try_generate_random_points::(50, (-2.0, 2.0)).unwrap(); assert_eq!(points.len(), 50); @@ -1447,7 +1530,7 @@ mod tests { #[test] fn test_generate_random_points_5d() { // Test 5D random point generation - let points = generate_random_points::(25, (-1.0, 1.0)).unwrap(); + let points = try_generate_random_points::(25, (-1.0, 1.0)).unwrap(); assert_eq!(points.len(), 25); @@ -1464,23 +1547,23 @@ mod tests { // Test invalid range (non-increasing bounds) across all dimensions // 2D - let result = generate_random_points::(100, (10.0, -10.0)); + let result = try_generate_random_points::(100, (10.0, -10.0)); assert_invalid_coordinate_range(&result, CoordinateRangeOrdering::Decreasing, 10.0, -10.0); // 3D - let result = generate_random_points::(50, (5.0, 5.0)); + let result = try_generate_random_points::(50, (5.0, 5.0)); assert_invalid_coordinate_range(&result, CoordinateRangeOrdering::Equal, 5.0, 5.0); // 4D - let result = generate_random_points::(25, (1.0, 0.5)); + let result = try_generate_random_points::(25, (1.0, 0.5)); assert_invalid_coordinate_range(&result, CoordinateRangeOrdering::Decreasing, 1.0, 0.5); // 5D - let result = generate_random_points::(10, (2.0, 2.0)); + let result = try_generate_random_points::(10, (2.0, 2.0)); assert_invalid_coordinate_range(&result, CoordinateRangeOrdering::Equal, 2.0, 2.0); // Test valid edge case - very small range - let points = generate_random_points::(10, (0.0, 0.001)).unwrap(); + let points = try_generate_random_points::(10, (0.0, 0.001)).unwrap(); assert_eq!(points.len(), 10); for point in points { for &coord in point.coords() { @@ -1492,16 +1575,16 @@ mod tests { #[test] fn test_generate_random_points_zero_points() { // Test generating zero points across all dimensions - let points_2d = generate_random_points::(0, (-1.0, 1.0)).unwrap(); + let points_2d = try_generate_random_points::(0, (-1.0, 1.0)).unwrap(); assert_eq!(points_2d.len(), 0); - let points_3d = generate_random_points::(0, (-1.0, 1.0)).unwrap(); + let points_3d = try_generate_random_points::(0, (-1.0, 1.0)).unwrap(); assert_eq!(points_3d.len(), 0); - let points_4d = generate_random_points::(0, (-1.0, 1.0)).unwrap(); + let points_4d = try_generate_random_points::(0, (-1.0, 1.0)).unwrap(); assert_eq!(points_4d.len(), 0); - let points_5d = generate_random_points::(0, (-1.0, 1.0)).unwrap(); + let points_5d = try_generate_random_points::(0, (-1.0, 1.0)).unwrap(); assert_eq!(points_5d.len(), 0); } @@ -1509,8 +1592,8 @@ mod tests { fn test_generate_random_points_seeded_2d() { // Test seeded 2D generation reproducibility let seed = 42_u64; - let points1 = generate_random_points_seeded::(50, (-5.0, 5.0), seed).unwrap(); - let points2 = generate_random_points_seeded::(50, (-5.0, 5.0), seed).unwrap(); + let points1 = try_generate_random_points_seeded::(50, (-5.0, 5.0), seed).unwrap(); + let points2 = try_generate_random_points_seeded::(50, (-5.0, 5.0), seed).unwrap(); assert_eq!(points1.len(), points2.len()); @@ -1529,8 +1612,8 @@ mod tests { fn test_generate_random_points_seeded_3d() { // Test seeded 3D generation reproducibility let seed = 123_u64; - let points1 = generate_random_points_seeded::(40, (0.0, 10.0), seed).unwrap(); - let points2 = generate_random_points_seeded::(40, (0.0, 10.0), seed).unwrap(); + let points1 = try_generate_random_points_seeded::(40, (0.0, 10.0), seed).unwrap(); + let points2 = try_generate_random_points_seeded::(40, (0.0, 10.0), seed).unwrap(); assert_eq!(points1.len(), points2.len()); @@ -1548,8 +1631,8 @@ mod tests { fn test_generate_random_points_seeded_4d() { // Test seeded 4D generation reproducibility let seed = 789_u64; - let points1 = generate_random_points_seeded::(30, (-2.5, 2.5), seed).unwrap(); - let points2 = generate_random_points_seeded::(30, (-2.5, 2.5), seed).unwrap(); + let points1 = try_generate_random_points_seeded::(30, (-2.5, 2.5), seed).unwrap(); + let points2 = try_generate_random_points_seeded::(30, (-2.5, 2.5), seed).unwrap(); assert_eq!(points1.len(), points2.len()); @@ -1567,8 +1650,8 @@ mod tests { fn test_generate_random_points_seeded_5d() { // Test seeded 5D generation reproducibility let seed = 456_u64; - let points1 = generate_random_points_seeded::(20, (-1.0, 3.0), seed).unwrap(); - let points2 = generate_random_points_seeded::(20, (-1.0, 3.0), seed).unwrap(); + let points1 = try_generate_random_points_seeded::(20, (-1.0, 3.0), seed).unwrap(); + let points2 = try_generate_random_points_seeded::(20, (-1.0, 3.0), seed).unwrap(); assert_eq!(points1.len(), points2.len()); @@ -1587,23 +1670,27 @@ mod tests { // Test that different seeds produce different results across all dimensions // 2D - let points1_2d = generate_random_points_seeded::(50, (0.0, 1.0), 42).unwrap(); - let points2_2d = generate_random_points_seeded::(50, (0.0, 1.0), 123).unwrap(); + let points1_2d = try_generate_random_points_seeded::(50, (0.0, 1.0), 42).unwrap(); + let points2_2d = try_generate_random_points_seeded::(50, (0.0, 1.0), 123).unwrap(); assert_ne!(points1_2d, points2_2d); // 3D - let points1_3d = generate_random_points_seeded::(30, (-5.0, 5.0), 42).unwrap(); - let points2_3d = generate_random_points_seeded::(30, (-5.0, 5.0), 999).unwrap(); + let points1_3d = try_generate_random_points_seeded::(30, (-5.0, 5.0), 42).unwrap(); + let points2_3d = try_generate_random_points_seeded::(30, (-5.0, 5.0), 999).unwrap(); assert_ne!(points1_3d, points2_3d); // 4D - let points1_4d = generate_random_points_seeded::(25, (-1.0, 1.0), 1337).unwrap(); - let points2_4d = generate_random_points_seeded::(25, (-1.0, 1.0), 7331).unwrap(); + let points1_4d = + try_generate_random_points_seeded::(25, (-1.0, 1.0), 1337).unwrap(); + let points2_4d = + try_generate_random_points_seeded::(25, (-1.0, 1.0), 7331).unwrap(); assert_ne!(points1_4d, points2_4d); // 5D - let points1_5d = generate_random_points_seeded::(15, (0.0, 10.0), 2021).unwrap(); - let points2_5d = generate_random_points_seeded::(15, (0.0, 10.0), 2024).unwrap(); + let points1_5d = + try_generate_random_points_seeded::(15, (0.0, 10.0), 2021).unwrap(); + let points2_5d = + try_generate_random_points_seeded::(15, (0.0, 10.0), 2024).unwrap(); assert_ne!(points1_5d, points2_5d); } #[test] @@ -1797,6 +1884,29 @@ mod tests { } } + #[test] + fn test_generate_random_points_in_ball_returns_typed_error_when_budget_exhausts() { + let mut rng = StdRng::seed_from_u64(42); + let result = + generate_random_points_in_ball_with_rng_and_budget::(1, 1.0, &mut rng, 0); + + let Err(RandomPointGenerationError::BallSamplingFailed { + requested_points, + generated_points, + dimension, + radius, + attempts, + }) = result + else { + panic!("expected ball sampling budget exhaustion"); + }; + assert_eq!(requested_points, 1); + assert_eq!(generated_points, 0); + assert_eq!(dimension, 2); + assert_relative_eq!(radius, 1.0, epsilon = f64::EPSILON); + assert_eq!(attempts, 0); + } + #[test] fn test_generate_random_points_in_ball_seeded_same_seed_is_deterministic_4d() { // This is a small smoke test that ensures deterministic output for fixed seed. @@ -1810,7 +1920,7 @@ mod tests { // Test that points cover the range reasonably well across all dimensions // 2D coverage test - let points_2d = generate_random_points::(500, (0.0, 10.0)).unwrap(); + let points_2d = try_generate_random_points::(500, (0.0, 10.0)).unwrap(); let mut min_2d = [f64::INFINITY; 2]; let mut max_2d = [f64::NEG_INFINITY; 2]; @@ -1835,7 +1945,7 @@ mod tests { } // 5D coverage test (smaller sample) - let points_5d = generate_random_points::(200, (-5.0, 5.0)).unwrap(); + let points_5d = try_generate_random_points::(200, (-5.0, 5.0)).unwrap(); let mut min_5d = [f64::INFINITY; 5]; let mut max_5d = [f64::NEG_INFINITY; 5]; @@ -1865,10 +1975,10 @@ mod tests { // Test common useful ranges across dimensions // Unit cube [0,1] for all dimensions - let unit_2d = generate_random_points::(50, (0.0, 1.0)).unwrap(); - let unit_3d = generate_random_points::(50, (0.0, 1.0)).unwrap(); - let unit_4d = generate_random_points::(50, (0.0, 1.0)).unwrap(); - let unit_5d = generate_random_points::(50, (0.0, 1.0)).unwrap(); + let unit_2d = try_generate_random_points::(50, (0.0, 1.0)).unwrap(); + let unit_3d = try_generate_random_points::(50, (0.0, 1.0)).unwrap(); + let unit_4d = try_generate_random_points::(50, (0.0, 1.0)).unwrap(); + let unit_5d = try_generate_random_points::(50, (0.0, 1.0)).unwrap(); assert_eq!(unit_2d.len(), 50); assert_eq!(unit_3d.len(), 50); @@ -1876,10 +1986,10 @@ mod tests { assert_eq!(unit_5d.len(), 50); // Centered cube [-1,1] for all dimensions - let centered_2d = generate_random_points::(30, (-1.0, 1.0)).unwrap(); - let centered_3d = generate_random_points::(30, (-1.0, 1.0)).unwrap(); - let centered_4d = generate_random_points::(30, (-1.0, 1.0)).unwrap(); - let centered_5d = generate_random_points::(30, (-1.0, 1.0)).unwrap(); + let centered_2d = try_generate_random_points::(30, (-1.0, 1.0)).unwrap(); + let centered_3d = try_generate_random_points::(30, (-1.0, 1.0)).unwrap(); + let centered_4d = try_generate_random_points::(30, (-1.0, 1.0)).unwrap(); + let centered_5d = try_generate_random_points::(30, (-1.0, 1.0)).unwrap(); assert_eq!(centered_2d.len(), 30); assert_eq!(centered_3d.len(), 30); @@ -2108,7 +2218,7 @@ mod tests { #[test] fn test_generate_poisson_points_2d() { // Test 2D Poisson disk sampling - let points = generate_poisson_points::(50, (0.0, 10.0), 0.5, 42).unwrap(); + let points = try_generate_poisson_points::(50, (0.0, 10.0), 0.5, 42).unwrap(); // Should generate some points (exact count depends on spacing constraints) assert!(!points.is_empty()); @@ -2141,7 +2251,7 @@ mod tests { #[test] fn test_generate_poisson_points_3d() { // Test 3D Poisson disk sampling - let points = generate_poisson_points::(30, (-1.0, 1.0), 0.2, 123).unwrap(); + let points = try_generate_poisson_points::(30, (-1.0, 1.0), 0.2, 123).unwrap(); assert!(!points.is_empty()); @@ -2174,7 +2284,7 @@ mod tests { #[test] fn test_generate_poisson_points_4d() { // Test 4D Poisson disk sampling - let points = generate_poisson_points::(15, (0.0, 5.0), 0.5, 333).unwrap(); + let points = try_generate_poisson_points::(15, (0.0, 5.0), 0.5, 333).unwrap(); assert!(!points.is_empty()); @@ -2207,7 +2317,7 @@ mod tests { #[test] fn test_generate_poisson_points_5d() { // Test 5D Poisson disk sampling - let points = generate_poisson_points::(10, (-2.0, 2.0), 0.4, 777).unwrap(); + let points = try_generate_poisson_points::(10, (-2.0, 2.0), 0.4, 777).unwrap(); assert!(!points.is_empty()); @@ -2241,8 +2351,8 @@ mod tests { #[test] fn test_generate_poisson_points_reproducible() { // Test that same seed produces same results - let points1 = generate_poisson_points::(25, (0.0, 5.0), 0.3, 456).unwrap(); - let points2 = generate_poisson_points::(25, (0.0, 5.0), 0.3, 456).unwrap(); + let points1 = try_generate_poisson_points::(25, (0.0, 5.0), 0.3, 456).unwrap(); + let points2 = try_generate_poisson_points::(25, (0.0, 5.0), 0.3, 456).unwrap(); assert_eq!(points1.len(), points2.len()); @@ -2256,18 +2366,18 @@ mod tests { } // Different seeds should produce different results - let points3 = generate_poisson_points::(25, (0.0, 5.0), 0.3, 789).unwrap(); + let points3 = try_generate_poisson_points::(25, (0.0, 5.0), 0.3, 789).unwrap(); assert_ne!(points1, points3); } #[test] fn test_generate_poisson_points_error_handling() { // Test invalid range - let result = generate_poisson_points::(50, (10.0, 5.0), 0.1, 42); + let result = try_generate_poisson_points::(50, (10.0, 5.0), 0.1, 42); assert_invalid_coordinate_range(&result, CoordinateRangeOrdering::Decreasing, 10.0, 5.0); // Test non-finite minimum distance rejects at the public boundary. - let nan_distance = generate_poisson_points::(0, (0.0, 1.0), f64::NAN, 42); + let nan_distance = try_generate_poisson_points::(0, (0.0, 1.0), f64::NAN, 42); assert_matches!( nan_distance, Err(RandomPointGenerationError::InvalidMinimumDistance { distance }) @@ -2275,7 +2385,7 @@ mod tests { ); let infinite_distance = - generate_poisson_points::(50, (0.0, 1.0), f64::INFINITY, 42); + try_generate_poisson_points::(50, (0.0, 1.0), f64::INFINITY, 42); assert_matches!( infinite_distance, Err(RandomPointGenerationError::InvalidMinimumDistance { distance }) @@ -2283,7 +2393,7 @@ mod tests { ); // Test minimum distance too large for bounds (should produce few/no points) - let result = generate_poisson_points::(100, (0.0, 1.0), 10.0, 42); + let result = try_generate_poisson_points::(100, (0.0, 1.0), 10.0, 42); match result { Ok(points) => { // Should produce very few points or fail @@ -2296,11 +2406,11 @@ mod tests { } // Test zero distance optimization (should return exact count without spacing checks) - let points = generate_poisson_points::(100, (0.0, 10.0), 0.0, 42).unwrap(); + let points = try_generate_poisson_points::(100, (0.0, 10.0), 0.0, 42).unwrap(); assert_eq!(points.len(), 100); // Should get exactly the requested number // Test negative distance optimization (should return exact count without spacing checks) - let points = generate_poisson_points::(50, (0.0, 10.0), -1.0, 42).unwrap(); + let points = try_generate_poisson_points::(50, (0.0, 10.0), -1.0, 42).unwrap(); assert_eq!(points.len(), 50); // Should get exactly the requested number } @@ -2311,40 +2421,40 @@ mod tests { #[test] fn test_generate_random_points_invalid_range() { // Test invalid range (non-increasing bounds) - let result = generate_random_points::(100, (10.0, 5.0)); + let result = try_generate_random_points::(100, (10.0, 5.0)); assert_matches!( result, Err(RandomPointGenerationError::InvalidCoordinateRange { .. }) ); // Test equal min and max - let result = generate_random_points::(100, (5.0, 5.0)); + let result = try_generate_random_points::(100, (5.0, 5.0)); assert_matches!( result, Err(RandomPointGenerationError::InvalidCoordinateRange { .. }) ); // Test valid range - let points = generate_random_points::(10, (0.0, 1.0)).unwrap(); + let points = try_generate_random_points::(10, (0.0, 1.0)).unwrap(); assert_eq!(points.len(), 10); } #[test] fn test_generate_random_points_seeded_invalid_range() { // Test invalid range with seed - let result = generate_random_points_seeded::(50, (100.0, 10.0), 42); + let result = try_generate_random_points_seeded::(50, (100.0, 10.0), 42); assert_matches!( result, Err(RandomPointGenerationError::InvalidCoordinateRange { .. }) ); // Test valid range produces consistent results - let points1 = generate_random_points_seeded::(5, (0.0, 1.0), 42).unwrap(); - let points2 = generate_random_points_seeded::(5, (0.0, 1.0), 42).unwrap(); + let points1 = try_generate_random_points_seeded::(5, (0.0, 1.0), 42).unwrap(); + let points2 = try_generate_random_points_seeded::(5, (0.0, 1.0), 42).unwrap(); assert_eq!(points1, points2); // Different seeds produce different results - let points3 = generate_random_points_seeded::(5, (0.0, 1.0), 123).unwrap(); + let points3 = try_generate_random_points_seeded::(5, (0.0, 1.0), 123).unwrap(); assert_ne!(points1, points3); } @@ -2367,21 +2477,21 @@ mod tests { #[test] fn test_generate_poisson_points_edge_cases() { // Test very small spacing with valid number of points - let result = generate_poisson_points::(100, (0.0, 1.0), 0.001, 42); + let result = try_generate_poisson_points::(100, (0.0, 1.0), 0.001, 42); if let Ok(points) = result { assert!(!points.is_empty()); } // May fail due to too many points // Test with zero points (should succeed with empty result) - let result = generate_poisson_points::(0, (0.0, 1.0), 0.1, 42); + let result = try_generate_poisson_points::(0, (0.0, 1.0), 0.1, 42); if let Ok(points) = result { assert!(points.is_empty()); } // Also acceptable if Err // Test very large spacing (should work but produce fewer points) - let result = generate_poisson_points::(10, (0.0, 1.0), 2.0, 42); + let result = try_generate_poisson_points::(10, (0.0, 1.0), 2.0, 42); if let Ok(points) = result { assert!(points.len() <= 10); } diff --git a/src/geometry/util/triangulation_generation.rs b/src/geometry/util/triangulation_generation.rs index 18bcfc8a..b2f29c21 100644 --- a/src/geometry/util/triangulation_generation.rs +++ b/src/geometry/util/triangulation_generation.rs @@ -291,7 +291,7 @@ where /// /// ```no_run /// use delaunay::prelude::construction::DelaunayTriangulationConstructionError; -/// use delaunay::prelude::generators::generate_random_triangulation; +/// use delaunay::prelude::generators::try_generate_random_triangulation; /// use std::num::NonZeroUsize; /// /// # fn main() -> Result<(), DelaunayTriangulationConstructionError> { @@ -305,7 +305,7 @@ where /// # return Ok(()); /// # }; /// // Generate a 2D triangulation with 50 points, no seed (random each time) -/// let triangulation_2d = generate_random_triangulation::<(), (), 2>( +/// let triangulation_2d = try_generate_random_triangulation::<(), (), 2>( /// fifty, /// (-10.0, 10.0), /// None, @@ -313,7 +313,7 @@ where /// ); /// /// // Generate a 3D triangulation with 30 points, seeded for reproducibility -/// let triangulation_3d = generate_random_triangulation::<(), (), 3>( +/// let triangulation_3d = try_generate_random_triangulation::<(), (), 3>( /// thirty, /// (-5.0, 5.0), /// None, @@ -321,7 +321,7 @@ where /// ); /// /// // Generate a 4D triangulation with custom vertex data -/// let triangulation_4d = generate_random_triangulation::( +/// let triangulation_4d = try_generate_random_triangulation::( /// twenty, /// (0.0, 1.0), /// Some(123), @@ -329,7 +329,7 @@ where /// ); /// /// // For string-like data, use fixed-size character arrays (Copy types) -/// let triangulation_with_strings = generate_random_triangulation::<[char; 8], (), 2>( +/// let triangulation_with_strings = try_generate_random_triangulation::<[char; 8], (), 2>( /// twenty, /// (0.0, 1.0), /// Some(['v', 'e', 'r', 't', 'e', 'x', '_', 'A']), @@ -361,11 +361,11 @@ where /// /// # See Also /// -/// - [`generate_random_points`](crate::geometry::util::generate_random_points) - For generating points without triangulation -/// - [`generate_random_points_seeded`](crate::geometry::util::generate_random_points_seeded) - For seeded random point generation only +/// - [`try_generate_random_points`](crate::geometry::util::try_generate_random_points) - For generating points without triangulation from raw bounds +/// - [`try_generate_random_points_seeded`](crate::geometry::util::try_generate_random_points_seeded) - For seeded random point generation from raw bounds /// - [`DelaunayTriangulationBuilder`](crate::DelaunayTriangulationBuilder) - For creating triangulations from existing vertices /// - [`RandomTriangulationBuilder`] - For more control over construction options -pub fn generate_random_triangulation( +pub fn try_generate_random_triangulation( n_points: NonZeroUsize, bounds: (f64, f64), vertex_data: Option, @@ -384,7 +384,7 @@ where n_points = n_points.get(), dimension = D, seed = ?seed, - "triangulation_generation::generate_random_triangulation called" + "triangulation_generation::try_generate_random_triangulation called" ); } let bounds = CoordinateRange::try_from(bounds) @@ -419,7 +419,7 @@ where /// /// ```no_run /// use delaunay::prelude::construction::DelaunayTriangulationConstructionError; -/// use delaunay::prelude::generators::generate_random_triangulation_with_topology_guarantee; +/// use delaunay::prelude::generators::try_generate_random_triangulation_with_topology_guarantee; /// use delaunay::prelude::TopologyGuarantee; /// use std::num::NonZeroUsize; /// @@ -427,7 +427,7 @@ where /// # let Some(twenty) = NonZeroUsize::new(20) else { /// # return Ok(()); /// # }; -/// let dt = generate_random_triangulation_with_topology_guarantee::<(), (), 3>( +/// let dt = try_generate_random_triangulation_with_topology_guarantee::<(), (), 3>( /// twenty, /// (-1.0, 1.0), /// None, @@ -438,7 +438,7 @@ where /// # Ok(()) /// # } /// ``` -pub fn generate_random_triangulation_with_topology_guarantee( +pub fn try_generate_random_triangulation_with_topology_guarantee( n_points: NonZeroUsize, bounds: (f64, f64), vertex_data: Option, @@ -691,7 +691,7 @@ where /// /// This builder provides a fluent API for constructing random triangulations with control over: /// - Insertion order strategy (`Input`, `Hilbert`) -/// - Topology guarantee (`None`, `PLManifold`) +/// - Topology guarantee (`Pseudomanifold`, `PLManifold`, `PLManifoldStrict`) /// - Construction options (deduplication, retry policy) /// - Topology/Euler validation (the final triangulation must pass Level-3 checks) /// @@ -1150,9 +1150,13 @@ mod tests { #[test] fn test_generate_random_triangulation_basic() { // Test 2D triangulation creation - let triangulation_2d = - generate_random_triangulation::<(), (), 2>(nonzero(10), (-5.0, 5.0), None, Some(42)) - .unwrap(); + let triangulation_2d = try_generate_random_triangulation::<(), (), 2>( + nonzero(10), + (-5.0, 5.0), + None, + Some(42), + ) + .unwrap(); assert!( triangulation_2d.number_of_vertices() >= 3, @@ -1163,7 +1167,7 @@ mod tests { triangulation_2d.is_valid().unwrap(); // Test 3D triangulation creation with data - let triangulation_3d = generate_random_triangulation::( + let triangulation_3d = try_generate_random_triangulation::( nonzero(8), (0.0, 1.0), Some(123), @@ -1180,13 +1184,21 @@ mod tests { triangulation_3d.is_valid().unwrap(); // Exercise repeatable construction with two deterministic seeds. - let triangulation_seeded = - generate_random_triangulation::<(), (), 2>(nonzero(5), (-1.0, 1.0), None, Some(789)) - .unwrap(); + let triangulation_seeded = try_generate_random_triangulation::<(), (), 2>( + nonzero(5), + (-1.0, 1.0), + None, + Some(789), + ) + .unwrap(); - let triangulation_different_seed = - generate_random_triangulation::<(), (), 2>(nonzero(5), (-1.0, 1.0), None, Some(790)) - .unwrap(); + let triangulation_different_seed = try_generate_random_triangulation::<(), (), 2>( + nonzero(5), + (-1.0, 1.0), + None, + Some(790), + ) + .unwrap(); triangulation_seeded.is_valid().unwrap(); triangulation_different_seed.is_valid().unwrap(); @@ -1204,7 +1216,7 @@ mod tests { #[test] fn test_generate_random_triangulation_error_cases() { - let result = generate_random_triangulation::<(), (), 2>( + let result = try_generate_random_triangulation::<(), (), 2>( nonzero(10), (5.0, 1.0), // min > max None, @@ -1233,7 +1245,7 @@ mod tests { #[test] fn test_generate_random_triangulation_rejects_nonfinite_bounds_as_generation_error() { - let nan_bounds = generate_random_triangulation::<(), (), 2>( + let nan_bounds = try_generate_random_triangulation::<(), (), 2>( nonzero(10), (f64::NAN, 1.0), None, @@ -1353,27 +1365,33 @@ mod tests { assert_relative_eq!(min, 5.0, epsilon = f64::EPSILON); assert_relative_eq!(max, 1.0, epsilon = f64::EPSILON); - let seeded_invalid_bounds = - RandomTriangulationBuilder::::try_new(nonzero(10), (5.0, 1.0)); - let Err(CoordinateRangeError::NonIncreasing { ordering, min, max }) = seeded_invalid_bounds - else { - panic!("expected seeded invalid bounds to fail"); + let equal_bounds = RandomTriangulationBuilder::::try_new(nonzero(10), (2.0, 2.0)); + let Err(CoordinateRangeError::NonIncreasing { ordering, min, max }) = equal_bounds else { + panic!("expected equal bounds to fail"); }; - assert_eq!(ordering, CoordinateRangeOrdering::Decreasing); - assert_relative_eq!(min, 5.0, epsilon = f64::EPSILON); - assert_relative_eq!(max, 1.0, epsilon = f64::EPSILON); + assert_eq!(ordering, CoordinateRangeOrdering::Equal); + assert_relative_eq!(min, 2.0, epsilon = f64::EPSILON); + assert_relative_eq!(max, 2.0, epsilon = f64::EPSILON); } #[test] fn test_generate_random_triangulation_reproducibility() { // Same seed should produce identical triangulations - let triangulation1 = - generate_random_triangulation::<(), (), 3>(nonzero(6), (-2.0, 2.0), None, Some(12345)) - .unwrap(); + let triangulation1 = try_generate_random_triangulation::<(), (), 3>( + nonzero(6), + (-2.0, 2.0), + None, + Some(12345), + ) + .unwrap(); - let triangulation2 = - generate_random_triangulation::<(), (), 3>(nonzero(6), (-2.0, 2.0), None, Some(12345)) - .unwrap(); + let triangulation2 = try_generate_random_triangulation::<(), (), 3>( + nonzero(6), + (-2.0, 2.0), + None, + Some(12345), + ) + .unwrap(); // Should have same structural properties assert_eq!( @@ -1415,30 +1433,46 @@ mod tests { // point sets in each dimension. // 2D with sufficient points for full triangulation - let tri_2d = - generate_random_triangulation::<(), (), 2>(nonzero(15), (0.0, 10.0), None, Some(555)) - .unwrap(); + let tri_2d = try_generate_random_triangulation::<(), (), 2>( + nonzero(15), + (0.0, 10.0), + None, + Some(555), + ) + .unwrap(); assert_eq!(tri_2d.dim(), 2); assert!(tri_2d.number_of_simplices() > 0); // 3D with sufficient points for full triangulation - let tri_3d = - generate_random_triangulation::<(), (), 3>(nonzero(20), (-3.0, 3.0), None, Some(666)) - .unwrap(); + let tri_3d = try_generate_random_triangulation::<(), (), 3>( + nonzero(20), + (-3.0, 3.0), + None, + Some(666), + ) + .unwrap(); assert_eq!(tri_3d.dim(), 3); assert!(tri_3d.number_of_simplices() > 0); // 4D with sufficient points for full triangulation - let tri_4d = - generate_random_triangulation::<(), (), 4>(nonzero(12), (-1.0, 1.0), None, Some(777)) - .unwrap(); + let tri_4d = try_generate_random_triangulation::<(), (), 4>( + nonzero(12), + (-1.0, 1.0), + None, + Some(777), + ) + .unwrap(); assert_eq!(tri_4d.dim(), 4); assert!(tri_4d.number_of_simplices() > 0); // 5D with sufficient points for full triangulation - let tri_5d = - generate_random_triangulation::<(), (), 5>(nonzero(10), (0.0, 5.0), None, Some(888)) - .unwrap(); + let tri_5d = try_generate_random_triangulation::<(), (), 5>( + nonzero(10), + (0.0, 5.0), + None, + Some(888), + ) + .unwrap(); assert_eq!(tri_5d.dim(), 5); assert!(tri_5d.number_of_simplices() > 0); } @@ -1450,7 +1484,7 @@ mod tests { // Test with fixed-size character array (Copy type that can represent strings) // NOTE: This is a workaround for the DataType trait requiring Copy, which // prevents using String or &str directly due to lifetime/ownership constraints - let tri_with_char_array = generate_random_triangulation::<[char; 8], (), 2>( + let tri_with_char_array = try_generate_random_triangulation::<[char; 8], (), 2>( nonzero(6), (-2.0, 2.0), Some(['v', 'e', 'r', 't', 'e', 'x', '_', 'd']), @@ -1476,7 +1510,7 @@ mod tests { None; let mut last_err: Option = None; for seed in seeds { - match generate_random_triangulation::( + match try_generate_random_triangulation::( nonzero(8), (0.0, 5.0), Some(42u32), @@ -1503,9 +1537,13 @@ mod tests { tri_with_int_data.tds().is_valid().unwrap(); // Test without data (None) - let tri_no_data = - generate_random_triangulation::<(), (), 2>(nonzero(5), (-1.0, 1.0), None, Some(111)) - .unwrap(); + let tri_no_data = try_generate_random_triangulation::<(), (), 2>( + nonzero(5), + (-1.0, 1.0), + None, + Some(111), + ) + .unwrap(); assert!( tri_no_data.number_of_vertices() >= 3, diff --git a/src/lib.rs b/src/lib.rs index 8b7a3c58..b35e36e4 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1063,9 +1063,12 @@ pub mod prelude { // In particular, exporting a local `uuid` module conflicts with the external `uuid` // crate name, making `use uuid::Uuid;` ambiguous for downstream users. pub use self::ordering::{ - HilbertBitDepth, HilbertError, MAX_HILBERT_BITS, hilbert_index, - hilbert_indices_prequantized, hilbert_quantize, hilbert_sort_by_stable, - hilbert_sort_by_unstable, hilbert_sorted_indices, + HilbertBitDepth, HilbertError, HilbertQuantizedBatch, HilbertQuantizedVec, + MAX_HILBERT_BITS, hilbert_index_in_range, hilbert_indices_for_quantized_batch, + hilbert_indices_prequantized, hilbert_quantize_batch_in_range, hilbert_quantize_in_range, + hilbert_sort_by_stable_in_range, hilbert_sort_by_unstable_in_range, + hilbert_sorted_indices_in_range, try_hilbert_index, try_hilbert_quantize, + try_hilbert_sort_by_stable, try_hilbert_sort_by_unstable, try_hilbert_sorted_indices, }; pub use crate::core::util::{ DeduplicationError, DelaunayValidationError, dedup_vertices_epsilon, dedup_vertices_exact, @@ -1615,14 +1618,14 @@ pub mod prelude { }; pub use crate::geometry::util::{ InvalidPositiveScalar, RandomPointGenerationError, RandomTriangulationBuilder, - generate_grid_points, generate_poisson_points, generate_poisson_points_in_range, - generate_random_points, generate_random_points_in_ball, + generate_grid_points, generate_poisson_points_in_range, generate_random_points_in_ball, generate_random_points_in_ball_seeded, generate_random_points_in_range, generate_random_points_in_range_seeded, generate_random_points_periodic, - generate_random_points_seeded, generate_random_triangulation, generate_random_triangulation_in_range, generate_random_triangulation_in_range_with_topology_guarantee, - generate_random_triangulation_with_topology_guarantee, scaled_bounds_by_point_count, + scaled_bounds_by_point_count, try_generate_poisson_points, try_generate_random_points, + try_generate_random_points_seeded, try_generate_random_triangulation, + try_generate_random_triangulation_with_topology_guarantee, }; } @@ -1635,20 +1638,24 @@ pub mod prelude { /// # Examples /// /// ```rust - /// use delaunay::prelude::ordering::{HilbertBitDepth, HilbertError, hilbert_sorted_indices}; + /// use delaunay::prelude::ordering::{HilbertBitDepth, HilbertError, try_hilbert_sorted_indices}; /// /// let coords = [[0.9_f64, 0.9], [0.1, 0.1], [0.5, 0.5]]; /// let bits = HilbertBitDepth::try_new(8)?; - /// let order = hilbert_sorted_indices(&coords, (0.0, 1.0), bits)?; + /// let order = try_hilbert_sorted_indices(&coords, (0.0, 1.0), bits)?; /// /// assert_eq!(order.len(), coords.len()); /// # Ok::<(), HilbertError>(()) /// ``` pub mod ordering { pub use crate::core::util::{ - HilbertBitDepth, HilbertError, MAX_HILBERT_BITS, hilbert_index, - hilbert_indices_prequantized, hilbert_quantize, hilbert_sort_by_stable, - hilbert_sort_by_unstable, hilbert_sorted_indices, + HilbertBitDepth, HilbertError, HilbertQuantizedBatch, HilbertQuantizedVec, + MAX_HILBERT_BITS, hilbert_index_in_range, hilbert_indices_for_quantized_batch, + hilbert_indices_prequantized, hilbert_quantize_batch_in_range, + hilbert_quantize_in_range, hilbert_sort_by_stable_in_range, + hilbert_sort_by_unstable_in_range, hilbert_sorted_indices_in_range, try_hilbert_index, + try_hilbert_quantize, try_hilbert_sort_by_stable, try_hilbert_sort_by_unstable, + try_hilbert_sorted_indices, }; } diff --git a/src/topology/traits/topological_space.rs b/src/topology/traits/topological_space.rs index ccdb0aab..d1a2e449 100644 --- a/src/topology/traits/topological_space.rs +++ b/src/topology/traits/topological_space.rs @@ -666,6 +666,7 @@ pub trait TopologicalSpace { #[cfg(test)] mod tests { use super::*; + use approx::assert_relative_eq; use std::assert_matches; #[test] @@ -920,6 +921,26 @@ mod tests { ); } + #[test] + fn test_toroidal_domain_try_from_and_into_periods_preserve_validation() { + let domain = ToroidalDomain::<3>::try_from([1.0, 2.0, 4.0]).unwrap(); + assert_relative_eq!(domain.periods()[0], 1.0); + assert_relative_eq!(domain.periods()[1], 2.0); + assert_relative_eq!(domain.periods()[2], 4.0); + + let periods = domain.into_periods(); + assert_relative_eq!(periods[0], 1.0); + assert_relative_eq!(periods[1], 2.0); + assert_relative_eq!(periods[2], 4.0); + + let invalid = ToroidalDomain::<3>::try_from([1.0, f64::NEG_INFINITY, 4.0]).unwrap_err(); + assert_matches!( + invalid, + ToroidalDomainError::InvalidPeriod { axis: 1, period } + if period.is_infinite() && period.is_sign_negative() + ); + } + #[test] fn test_global_topology_try_toroidal_parses_domain() { let topology = diff --git a/tests/delaunay_edge_cases.rs b/tests/delaunay_edge_cases.rs index 7b770ccc..3dc94f01 100644 --- a/tests/delaunay_edge_cases.rs +++ b/tests/delaunay_edge_cases.rs @@ -14,7 +14,10 @@ use delaunay::prelude::construction::{ }; #[cfg(feature = "diagnostics")] use delaunay::prelude::diagnostics::debug_print_first_delaunay_violation; -use delaunay::prelude::generators::generate_random_points_in_ball_seeded; +use delaunay::prelude::generators::{ + generate_random_points_in_ball_seeded, + try_generate_random_triangulation_with_topology_guarantee, +}; use delaunay::prelude::geometry::RobustKernel; use rand::SeedableRng; use rand::seq::SliceRandom; @@ -473,11 +476,7 @@ fn test_regression_non_manifold_3d_seed123_50pts() { // Exact configuration from CI failure (matches ci_performance_suite.rs) let n_points = nonzero(50); let raw_n_points = n_points.get(); - let result = delaunay::geometry::util::generate_random_triangulation_with_topology_guarantee::< - (), - (), - 3, - >( + let result = try_generate_random_triangulation_with_topology_guarantee::<(), (), 3>( n_points, // Point count from CI benchmark (-100.0, 100.0), // Bounds from benchmark None, // No vertex data @@ -529,18 +528,13 @@ fn test_regression_non_manifold_nearby_seeds() { let min_vertices = (raw_n_points / 6).max(4); for seed in test_seeds { - let result = - delaunay::geometry::util::generate_random_triangulation_with_topology_guarantee::< - (), - (), - 3, - >( - n_points, - (-100.0, 100.0), - None, - Some(seed), - TopologyGuarantee::PLManifold, - ); + let result = try_generate_random_triangulation_with_topology_guarantee::<(), (), 3>( + n_points, + (-100.0, 100.0), + None, + Some(seed), + TopologyGuarantee::PLManifold, + ); assert!( result.is_ok(), diff --git a/tests/large_scale_debug.rs b/tests/large_scale_debug.rs index f59dd9a1..8ad19adb 100644 --- a/tests/large_scale_debug.rs +++ b/tests/large_scale_debug.rs @@ -100,15 +100,17 @@ #![forbid(unsafe_code)] use delaunay::geometry::kernel::{ExactPredicates, Kernel, RobustKernel}; -use delaunay::geometry::util::{ - generate_random_points_in_ball_seeded, generate_random_points_seeded, safe_usize_to_scalar, -}; +use delaunay::geometry::util::safe_usize_to_scalar; use delaunay::prelude::construction::{ ConstructionOptions, ConstructionStatistics, DelaunayRepairPolicy, DelaunayTriangulation, DelaunayTriangulationConstructionErrorWithStatistics, InitialSimplexStrategy, TopologyGuarantee, Vertex, vertex, }; use delaunay::prelude::diagnostics::ConstructionTelemetry; +use delaunay::prelude::generators::{ + generate_random_points_in_ball_seeded, generate_random_points_in_range_seeded, +}; +use delaunay::prelude::geometry::CoordinateRange; #[cfg(feature = "diagnostics")] use delaunay::prelude::insertion::InsertionResult; use delaunay::prelude::insertion::{InsertionOutcome, InsertionStatistics}; @@ -1264,10 +1266,9 @@ where }) } PointDistribution::Box => { - let range = (-box_half_width, box_half_width); - generate_random_points_seeded::(n_points, range, seed).unwrap_or_else(|e| { - panic!("failed to generate deterministic box points (range={range:?}): {e}") - }) + let range = CoordinateRange::try_new(-box_half_width, box_half_width) + .expect("box half-width should define a finite non-empty coordinate range"); + generate_random_points_in_range_seeded::(n_points, range, seed) } }; println!("Generated {} points in {:?}", points.len(), t_gen.elapsed()); diff --git a/tests/prelude_exports.rs b/tests/prelude_exports.rs index 09edaed9..e6b8cbc6 100644 --- a/tests/prelude_exports.rs +++ b/tests/prelude_exports.rs @@ -35,7 +35,7 @@ use delaunay::prelude::construction::{ InvalidPositiveScalar as ConstructionInvalidPositiveScalar, RandomPointGenerationError, SimplexValidationError, SpatialIndexConstructionFailure as ConstructionSpatialIndexConstructionFailure, - TopologyGuarantee, Vertex, vertex, + TopologyGuarantee, ToroidalDomain as ConstructionToroidalDomain, Vertex, vertex, }; use delaunay::prelude::delaunayize::{ DelaunayTriangulationBuilder as DelaunayizeDelaunayTriangulationBuilder, DelaunayizeConfig, @@ -51,7 +51,8 @@ use delaunay::prelude::diagnostics::{ use delaunay::prelude::flips::BistellarFlips; use delaunay::prelude::generators::{ CoordinateRange, CoordinateRangeError, InvalidPositiveScalar, RandomTriangulationBuilder, - generate_grid_points, generate_random_points_in_range_seeded, generate_random_points_seeded, + generate_grid_points, generate_random_points_in_range_seeded, + try_generate_random_points_seeded, }; #[cfg(feature = "diagnostics")] use delaunay::prelude::geometry::{AdaptiveKernel, Coordinate}; @@ -67,8 +68,11 @@ use delaunay::prelude::insertion::{ NeighborRebuildError, Tds as InsertionTds, TdsMutationError, repair_neighbor_pointers_local, }; use delaunay::prelude::ordering::{ - HilbertBitDepth, HilbertError, MAX_HILBERT_BITS, hilbert_index, hilbert_indices_prequantized, - hilbert_quantize, hilbert_sort_by_stable, hilbert_sort_by_unstable, hilbert_sorted_indices, + HilbertBitDepth, HilbertError, HilbertQuantizedBatch, MAX_HILBERT_BITS, hilbert_index_in_range, + hilbert_indices_for_quantized_batch, hilbert_indices_prequantized, hilbert_quantize_in_range, + hilbert_sort_by_stable_in_range, hilbert_sort_by_unstable_in_range, + hilbert_sorted_indices_in_range, try_hilbert_index, try_hilbert_quantize, + try_hilbert_sort_by_stable, try_hilbert_sort_by_unstable, try_hilbert_sorted_indices, }; use delaunay::prelude::query::{ConvexHull, QueryError}; use delaunay::prelude::repair::{ @@ -87,7 +91,8 @@ use delaunay::prelude::topology::spaces::{ GlobalTopology, ToroidalConstructionMode, ToroidalDomain, ToroidalDomainError, }; use delaunay::prelude::topology::validation::{ - ManifoldError, RidgeVertices, RidgeVerticesError, ridge_star_simplices, + GlobalTopology as TopologyValidationGlobalTopology, ManifoldError, RidgeVertices, + RidgeVerticesError, ridge_star_simplices, }; use delaunay::prelude::triangulation::{ FacetIssuesMap as TriangulationFacetIssuesMap, FastKernel as TriangulationFastKernel, @@ -222,7 +227,8 @@ fn root_exports_cover_flattened_public_api() -> Result<(), RootApiExportTestErro #[test] fn preludes_cover_bench_apis() -> Result<(), PreludeExportTestError> { - let _generated_points: Vec> = generate_random_points_seeded(3, (0.0, 1.0), 42)?; + let _generated_points: Vec> = + try_generate_random_points_seeded(3, (0.0, 1.0), 42)?; let vertices: Vec> = vec![ vertex!([0.0, 0.0, 0.0]), @@ -669,6 +675,12 @@ fn topology_spaces_prelude_covers_toroidal_domain_api() -> Result<(), PreludeExp assert_relative_eq!(domain.periods()[2], 3.0); assert_eq!(domain.period(1), Some(2.0)); + let construction_domain = ConstructionToroidalDomain::<3>::try_new([1.0, 2.0, 3.0])?; + assert_relative_eq!(construction_domain.periods()[2], 3.0); + + let validation_topology = TopologyValidationGlobalTopology::<3>::default(); + assert!(validation_topology.is_euclidean()); + let topology = GlobalTopology::try_toroidal( [1.0, 2.0, 3.0], ToroidalConstructionMode::PeriodicImagePoint, @@ -862,26 +874,47 @@ fn ordering_prelude_covers_hilbert_apis() -> Result<(), HilbertError> { let coords = [[0.9_f64, 0.9], [0.1, 0.1], [0.5, 0.5]]; assert_eq!(MAX_HILBERT_BITS, 31); let bits = HilbertBitDepth::try_new(8)?; - let order = hilbert_sorted_indices(&coords, (0.0, 1.0), bits)?; + let order = try_hilbert_sorted_indices(&coords, (0.0, 1.0), bits)?; assert_eq!(order.len(), coords.len()); + let bounds = CoordinateRange::try_new(0.0_f64, 1.0).expect("test bounds should be valid"); + let range_order = hilbert_sorted_indices_in_range(&coords, bounds, bits)?; + assert_eq!(range_order, order); let quantized: Vec<[u32; 2]> = coords .iter() - .map(|coord| hilbert_quantize(coord, (0.0, 1.0), bits)) + .map(|coord| try_hilbert_quantize(coord, (0.0, 1.0), bits)) .collect::>()?; + let range_quantized = hilbert_quantize_in_range(&coords[0], bounds, bits)?; + assert_eq!(range_quantized, quantized[0]); let indices = hilbert_indices_prequantized(&quantized, bits)?; assert_eq!(indices.len(), coords.len()); + let quantized_batch = HilbertQuantizedBatch::try_new(&quantized, bits)?; + assert_eq!(quantized_batch.coordinates(), quantized.as_slice()); + assert_eq!(quantized_batch.bits(), bits); + assert_eq!(quantized_batch.indices(), indices); + assert_eq!( + hilbert_indices_for_quantized_batch(quantized_batch), + indices + ); - let index = hilbert_index(&coords[0], (0.0, 1.0), bits)?; + let index = try_hilbert_index(&coords[0], (0.0, 1.0), bits)?; assert_eq!(index, indices[0]); + let range_index = hilbert_index_in_range(&coords[0], bounds, bits)?; + assert_eq!(range_index, index); let mut stable_payload = vec![0_usize, 1, 2]; - hilbert_sort_by_stable(&mut stable_payload, (0.0, 1.0), bits, |&i| coords[i])?; + try_hilbert_sort_by_stable(&mut stable_payload, (0.0, 1.0), bits, |&i| coords[i])?; assert_eq!(stable_payload, order); + let mut range_stable_payload = vec![0_usize, 1, 2]; + hilbert_sort_by_stable_in_range(&mut range_stable_payload, bounds, bits, |&i| coords[i])?; + assert_eq!(range_stable_payload, order); let mut unstable_payload = vec![0_usize, 1, 2]; - hilbert_sort_by_unstable(&mut unstable_payload, (0.0, 1.0), bits, |&i| coords[i])?; + try_hilbert_sort_by_unstable(&mut unstable_payload, (0.0, 1.0), bits, |&i| coords[i])?; assert_eq!(unstable_payload, order); + let mut range_unstable_payload = vec![0_usize, 1, 2]; + hilbert_sort_by_unstable_in_range(&mut range_unstable_payload, bounds, bits, |&i| coords[i])?; + assert_eq!(range_unstable_payload, order); Ok(()) } diff --git a/tests/proptest_euler_characteristic.rs b/tests/proptest_euler_characteristic.rs index 906a2f9c..cf9d592d 100644 --- a/tests/proptest_euler_characteristic.rs +++ b/tests/proptest_euler_characteristic.rs @@ -16,8 +16,8 @@ //! //! For deterministic tests with known configurations, see `euler_characteristic.rs`. -use delaunay::geometry::util::generate_random_triangulation_with_topology_guarantee; use delaunay::prelude::construction::{DelaunayTriangulation, TopologyGuarantee, vertex}; +use delaunay::prelude::generators::try_generate_random_triangulation_with_topology_guarantee; use delaunay::topology::characteristics::{euler, validation}; use proptest::prelude::*; use std::num::NonZeroUsize; @@ -171,7 +171,7 @@ macro_rules! test_euler_properties { #[test] fn test_seeded_random_generator_euler_consistent() { - let dt_2d = generate_random_triangulation_with_topology_guarantee::<(), (), 2>( + let dt_2d = try_generate_random_triangulation_with_topology_guarantee::<(), (), 2>( nonzero(15), (0.0, 10.0), None, @@ -190,7 +190,7 @@ fn test_seeded_random_generator_euler_consistent() { result_2d.counts.count(2), ); - let dt_3d = generate_random_triangulation_with_topology_guarantee::<(), (), 3>( + let dt_3d = try_generate_random_triangulation_with_topology_guarantee::<(), (), 3>( nonzero(20), (-3.0, 3.0), None, diff --git a/tests/regressions.rs b/tests/regressions.rs index 3857258c..2b1a99af 100644 --- a/tests/regressions.rs +++ b/tests/regressions.rs @@ -11,12 +11,13 @@ use delaunay::prelude::construction::{ #[cfg(feature = "diagnostics")] use delaunay::prelude::diagnostics::debug_print_first_delaunay_violation; use delaunay::prelude::generators::generate_random_points_in_ball_seeded; -use delaunay::prelude::geometry::{Point, RobustKernel}; +use delaunay::prelude::geometry::{Coordinate, CoordinateRange, Point, RobustKernel}; use delaunay::prelude::insertion::{ HullExtensionReason, InsertionError, InsertionErrorKind, InsertionErrorSummary, }; use delaunay::prelude::ordering::{ - HilbertBitDepth, hilbert_indices_prequantized, hilbert_quantize, + HilbertBitDepth, hilbert_indices_prequantized, hilbert_quantize_batch_in_range, + hilbert_quantize_in_range, }; /// Replays a full Hilbert ordering while keeping only the prefix that first @@ -25,12 +26,12 @@ fn hilbert_ordered_prefix( points: Vec>, prefix_len: usize, ) -> Vec> { - let (min, max) = coordinate_bounds(&points); + let bounds = coordinate_bounds(&points); let bits_per_coord = HilbertBitDepth::try_new(31).expect("test bit depth must be valid"); let quantized: Vec<[u32; D]> = points .iter() .map(|point| { - hilbert_quantize(point.coords(), (min, max), bits_per_coord) + hilbert_quantize_in_range(point.coords(), bounds, bits_per_coord) .expect("finite generated points should quantize") }) .collect(); @@ -73,13 +74,82 @@ fn hilbert_ordered_prefix( /// Computes the scalar range used by batch Hilbert ordering so regression /// prefixes match the original full construction order. -fn coordinate_bounds(points: &[Point]) -> (f64, f64) { - points +fn coordinate_bounds(points: &[Point]) -> CoordinateRange { + let (min, max) = points .iter() .flat_map(Point::coords) .fold((f64::INFINITY, f64::NEG_INFINITY), |(min, max), &coord| { (min.min(coord), max.max(coord)) - }) + }); + CoordinateRange::try_new(min, max) + .expect("generated regression points should span a finite non-empty range") +} + +/// Locks the equivalence between the single-pass proof-carrying batch quantizer +/// used by Hilbert construction ordering and the original two-step +/// `quantize` + `hilbert_indices_prequantized` path. +/// +/// `order_vertices_hilbert` switched to `hilbert_quantize_batch_in_range` to +/// drop a redundant per-coordinate range rescan (and a per-point quantization +/// scale recompute). This regression guards that the change does not alter the +/// quantized cells or Hilbert indices — and therefore the deterministic +/// insertion order — across representative dimensions and adversarial inputs. +#[test] +fn regression_hilbert_batch_quantize_matches_two_step_path() { + fn assert_paths_match(points: &[Point]) { + let bounds = coordinate_bounds(points); + // Mirror `order_vertices_hilbert`'s per-dimension precision so the + // `D * bits <= 128` index-width invariant holds for every D. + let bits_per_coord = (128_u32 / u32::try_from(D).expect("dimension fits in u32")).min(31); + let bits = HilbertBitDepth::try_new(bits_per_coord).expect("test bit depth must be valid"); + + // Original two-step path: per-point quantize, then bulk index. + let two_step_quantized: Vec<[u32; D]> = points + .iter() + .map(|point| { + hilbert_quantize_in_range(point.coords(), bounds, bits) + .expect("finite points should quantize") + }) + .collect(); + let two_step_indices = hilbert_indices_prequantized(&two_step_quantized, bits) + .expect("indices should fit in u128"); + + // New single-pass proof-carrying batch path. + let batch = hilbert_quantize_batch_in_range(points, bounds, bits, |point| *point.coords()) + .expect("finite points should quantize"); + let (batch_indices, batch_quantized) = batch.into_indices_and_coordinates(); + + assert_eq!( + batch_quantized, two_step_quantized, + "batch quantizer must produce identical quantized cells in {D}D" + ); + assert_eq!( + batch_indices, two_step_indices, + "batch quantizer must produce identical Hilbert indices in {D}D" + ); + } + + // Adversarial mixes: negative/asymmetric ranges, clamping at both ends, + // duplicate cells, and exact endpoints. + assert_paths_match::<2>(&[ + Point::new([-2.0, -1.0]), + Point::new([-1.5, 0.25]), + Point::new([0.1, -0.7]), + Point::new([3.0, 3.0]), + Point::new([3.0, 3.0]), + ]); + assert_paths_match::<3>(&[ + Point::new([-2.0, -1.0, 0.0]), + Point::new([-1.5, 0.25, 1.75]), + Point::new([0.1, -0.7, 2.2]), + Point::new([3.0, 3.0, -2.0]), + ]); + assert_paths_match::<5>(&[ + Point::new([-2.0, -1.0, 0.0, 1.0, 2.0]), + Point::new([-1.5, 0.25, 1.75, 2.5, -0.5]), + Point::new([0.1, -0.7, 2.2, -1.8, 1.4]), + Point::new([3.0, 3.0, -2.0, -2.0, 0.5]), + ]); } #[test]