Recover from numerical instability.

RELEASES.md

@@ -1,3 +1,8 @@
+# Version 0.4.1 (2024-03-26)
+
+ * Leverage move-to-front heuristic to prevent panic due to numerical instability.
+ * Let `Enclosing::contains` panic for infinite points.
+
 # Version 0.4.0 (2024-03-24)
 
  * Lower trait bounds on `Ball` and `Enclosing`.

src/ball.rs

@@ -34,7 +34,9 @@ where
 {
  #[inline]
  fn contains(&self, point: &OPoint<T, D>) -> bool {
- (point - &self.center).norm_squared() <= self.radius_squared
+ let norm_squared = (point - &self.center).norm_squared();
+ assert!(norm_squared.is_finite());
+ norm_squared <= self.radius_squared
  }
  fn with_bounds(bounds: &[OPoint<T, D>]) -> Option<Self>
  where
@@ -65,16 +67,17 @@ where
  }
  });
  let vector = vector.view((0, 0), (length, 1));
- matrix.try_inverse().map(|matrix| {
+ matrix.try_inverse().and_then(|matrix| {
  let vector = matrix * vector;
  let mut center = OVector::<T, D>::zeros();
  for point in 0..length {
  center += points.column(point) * vector[point].clone();
  }
- Self {
+ let radius_squared = center.norm_squared();
+ radius_squared.is_finite().then(|| Self {
  center: &bounds[0] + &center,
- radius_squared: center.norm_squared(),
- }
+ radius_squared,
+ })
  })
  }
 }

src/enclosing.rs

@@ -147,13 +147,18 @@ where
  DefaultAllocator: Allocator<T, D, D> + Allocator<OPoint<T, D>, DimNameSum<D, U1>>,
  <DefaultAllocator as Allocator<OPoint<T, D>, DimNameSum<D, U1>>>::Buffer: Default,
  {
- maybe_grow(Self::RED_ZONE, Self::STACK_SIZE, || {
- Self::enclosing_points_with_bounds(
- points,
- &mut OVec::<OPoint<T, D>, DimNameSum<D, U1>>::new(),
- )
- .expect("Empty point set")
- })
+ 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::<OPoint<T, D>, DimNameSum<D, U1>>::new(),
+ )
+ }) {
+ return ball;
+ }
+ }
+ unreachable!("numerical instability");
  }
  /// Returns minimum ball enclosing `points` with `bounds`.
  ///