diff --git a/Cargo.toml b/Cargo.toml index 7689221..8e3e60d 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -2,7 +2,7 @@ name = "miniball" description = "Minimum enclosing ball" authors = ["Rouven Spreckels "] -version = "0.4.0" +version = "0.4.1" rust-version = "1.61.0" edition = "2021" documentation = "https://docs.rs/miniball" @@ -37,7 +37,7 @@ all-features = true rustdoc-args = ["--cfg", "docsrs"] [dependencies] -nalgebra = { version = "0.32.4", default-features = false, features = ["alloc"] } +nalgebra = { version = "0.32.5", default-features = false, features = ["alloc"] } stacker = { version = "0.1.15", optional = true } [features] @@ -45,7 +45,9 @@ default = ["std"] std = ["dep:stacker"] [dev-dependencies] -nalgebra = { version = "0.32.4", features = ["rand"] } +nalgebra = { version = "0.32.5", features = ["alloc", "rand"] } +rand_distr = { version = "0.4.3", default-features = false } +rand = { version = "0.8.5", default-features = false } [profile.test] opt-level = 2 diff --git a/README.md b/README.md index 3356077..4ea86ee 100644 --- a/README.md +++ b/README.md @@ -25,6 +25,10 @@ Minimum enclosing ball. * Find minimum-volume enclosing *n*-ellipsoid. * Improve numerical stability and performance. +# Features + + * `std` for spilling recursion stack over to the heap if necessary. Enabled by `default`. + See the [release history] to keep track of the development. [release history]: RELEASES.md diff --git a/RELEASES.md b/RELEASES.md index 15de415..08bcd79 100644 --- a/RELEASES.md +++ b/RELEASES.md @@ -1,5 +1,6 @@ -# Version 0.4.1 (2024-03-26) +# Version 0.4.1 (2024-04-12) + * Attempt to improve numerical stability by enclosing approximately co-spherical points. * Leverage move-to-front heuristic to prevent panic due to numerical instability. * Let `Enclosing::contains` panic for infinite points. diff --git a/examples/co-spherical.rs b/examples/co-spherical.rs new file mode 100644 index 0000000..91c2269 --- /dev/null +++ b/examples/co-spherical.rs @@ -0,0 +1,34 @@ +use miniball::{ + nalgebra::{Point3, Vector3}, + {Ball, Enclosing}, +}; +use rand::thread_rng; +use rand_distr::{Distribution, UnitSphere}; +use std::collections::VecDeque; + +fn main() { + let m = 100_000; + println!("n = 3, m = {m}"); + println!(); + let center = Vector3::new(-3.0, 7.0, 4.8); + let radius = 3.0; + let radius_squared = radius * radius; + let mut points = UnitSphere + .sample_iter(&mut thread_rng()) + .take(m) + .map(Point3::::from) + .map(|point| point * radius + center) + .collect::>(); + let ball = (0..8) + .map(|_| { + let ball = Ball::enclosing_points(&mut points); + let epsilon = ball.radius_squared / radius_squared - 1.0; + println!("sample with accuracy: 1{epsilon:+.1e}"); + ball + }) + .min() + .unwrap(); + println!(); + let epsilon = ball.radius_squared / radius_squared - 1.0; + println!("result with accuracy: 1{epsilon:+.1e}"); +} diff --git a/src/ball.rs b/src/ball.rs index b18fcf0..f15447f 100644 --- a/src/ball.rs +++ b/src/ball.rs @@ -5,12 +5,13 @@ // file, You can obtain one at https://mozilla.org/MPL/2.0/. use super::Enclosing; +use core::cmp::Ordering; use nalgebra::{ base::allocator::Allocator, DefaultAllocator, DimName, OMatrix, OPoint, OVector, RealField, }; /// Ball over real field `T` of dimension `D` with center and radius squared. -#[derive(Debug, Clone, PartialEq)] +#[derive(Debug, Clone)] pub struct Ball where DefaultAllocator: Allocator, @@ -28,6 +29,41 @@ where { } +impl PartialEq for Ball +where + DefaultAllocator: Allocator, +{ + fn eq(&self, other: &Self) -> bool { + assert!( + self.radius_squared.is_finite() && other.radius_squared.is_finite(), + "infinite ball" + ); + self.radius_squared == other.radius_squared + } +} + +impl Eq for Ball where DefaultAllocator: Allocator {} + +impl PartialOrd for Ball +where + DefaultAllocator: Allocator, +{ + fn partial_cmp(&self, other: &Self) -> Option { + Some(self.cmp(other)) + } +} + +impl Ord for Ball +where + DefaultAllocator: Allocator, +{ + fn cmp(&self, other: &Self) -> Ordering { + self.radius_squared + .partial_cmp(&other.radius_squared) + .expect("infinite ball") + } +} + impl Enclosing for Ball where DefaultAllocator: Allocator, @@ -35,8 +71,8 @@ where #[inline] fn contains(&self, point: &OPoint) -> bool { let norm_squared = (point - &self.center).norm_squared(); - assert!(norm_squared.is_finite()); - norm_squared <= self.radius_squared + assert!(norm_squared.is_finite(), "infinite point"); + self.radius_squared.clone() / norm_squared >= T::one() - T::default_epsilon().sqrt() } fn with_bounds(bounds: &[OPoint]) -> Option where @@ -74,8 +110,9 @@ where center += points.column(point) * vector[point].clone(); } let radius_squared = center.norm_squared(); + let center = &bounds[0] + ¢er; radius_squared.is_finite().then(|| Self { - center: &bounds[0] + ¢er, + center, radius_squared, }) }) diff --git a/src/enclosing.rs b/src/enclosing.rs index 9133eb9..9d743d7 100644 --- a/src/enclosing.rs +++ b/src/enclosing.rs @@ -91,6 +91,11 @@ where /// *n*-dimensional points. The complexity constant in *m* is significantly reduced by reusing /// permuted points of previous invocations. /// + /// # Stability + /// + /// Due to floating-point inaccuracies, the returned ball might not exactly be the minimum for + /// degenerate (e.g., co-spherical) `points`. + /// /// # Example /// /// Finds minimum 4-ball enclosing 4-cube (tesseract): @@ -148,17 +153,14 @@ where , DimNameSum>>::Buffer: Default, { assert!(!points.is_empty(), "empty point set"); - for _ in 0..points.len() { - if let Some(ball) = maybe_grow(Self::RED_ZONE, Self::STACK_SIZE, || { - Self::enclosing_points_with_bounds( - points, - &mut OVec::, DimNameSum>::new(), - ) - }) { - return ball; - } - } - unreachable!("numerical instability"); + let mut bounds = OVec::, DimNameSum>::new(); + (0..bounds.capacity()) + .find_map(|_| { + maybe_grow(Self::RED_ZONE, Self::STACK_SIZE, || { + Self::enclosing_points_with_bounds(points, &mut bounds) + }) + }) + .expect("numerical instability") } /// Returns minimum ball enclosing `points` with `bounds`. /// diff --git a/src/ovec.rs b/src/ovec.rs index 1815db9..54f651d 100644 --- a/src/ovec.rs +++ b/src/ovec.rs @@ -30,6 +30,12 @@ where pub fn new() -> Self { Self::default() } + /// Maximum number of items. + #[must_use] + #[inline] + pub fn capacity(&self) -> usize { + self.data.len() + } /// Number of items. #[must_use] #[inline] diff --git a/tests/ball_enclosing_points.rs b/tests/ball_enclosing_points.rs index 3840100..3df7ff9 100644 --- a/tests/ball_enclosing_points.rs +++ b/tests/ball_enclosing_points.rs @@ -106,7 +106,7 @@ fn minimum_3_ball_enclosing_3_line() { #[test] fn minimum_6_ball_enclosing_6_cube() { for _randomize in 0..100 { - // Uniform distribution in 4-cube centered around `offset` with room `diagonal_halved`. + // Uniform distribution in 6-cube centered around `offset` with room `diagonal_halved`. let offset = Vector6::new(-3.0, 7.0, 4.8, 1.2, 5.3, 7.4); let diagonal_halved = 3.0; let mut points = (0..10_000) @@ -115,26 +115,26 @@ fn minimum_6_ball_enclosing_6_cube() { .map(|point| point + offset) .collect::>(); for _reuse in 0..10 { - // Computes 4-ball enclosing 4-cube. + // Computes 6-ball enclosing 6-cube. let Ball { center, radius_squared, } = Ball::enclosing_points(&mut points); let radius = radius_squared.sqrt(); - // Ensures enclosing 4-ball is roughly centered around uniform distribution in 4-cube + // Ensures enclosing 6-ball is roughly centered around uniform distribution in 6-cube // and radius roughly matches room diagonal halved, guaranteeing certain uniformity of // randomly distributed points. assert!((center - offset).map(f64::abs) < Vector6::from_element(1.0).into()); assert!((radius - diagonal_halved).abs() < 1.0); - // Epsilon of numerical stability for computing circumscribed 4-ball. This is related to + // Epsilon of numerical stability for computing circumscribed 6-ball. This is related to // robustness of `Enclosing::with_bounds()` regarding floating-point inaccuracies. let epsilon = f64::EPSILON.sqrt(); - // Ensures all points are enclosed by 4-ball. + // Ensures all points are enclosed by 6-ball. let all_enclosed = points .iter() .all(|point| distance(point, ¢er) <= radius + epsilon); assert!(all_enclosed); - // Ensures at least 2 points are on surface of 4-ball, mandatory to be minimum. + // Ensures at least 2 points are on surface of 6-ball, mandatory to be minimum. let bounds_count = points .iter() .map(|point| distance(point, ¢er))