From 7e6df5fa51392f9d6ac28fa84eb1dda0ccabe5c8 Mon Sep 17 00:00:00 2001 From: janbridley Date: Thu, 18 Jun 2026 13:11:16 -0400 Subject: [PATCH 01/10] clean f32 support --- hoomd-geometry/src/shape/convex_polytope.rs | 31 +++++++++++++++------ 1 file changed, 22 insertions(+), 9 deletions(-) diff --git a/hoomd-geometry/src/shape/convex_polytope.rs b/hoomd-geometry/src/shape/convex_polytope.rs index 7db4a0e38..58b0a2e03 100644 --- a/hoomd-geometry/src/shape/convex_polytope.rs +++ b/hoomd-geometry/src/shape/convex_polytope.rs @@ -207,23 +207,36 @@ impl ConvexPolytope } } +/// Compute the dot product of two `[f32; N]` arrays. +#[inline(always)] +fn dot_arrays(l: &[f32; N], r: &[f32; N]) -> f32 { + (0..N).map(|i| l[i] * r[i]).sum() +} + impl SupportMapping> for ConvexPolytope { #[inline] + #[expect(clippy::cast_possible_truncation, reason = "Truncation is ok.")] fn support_mapping(&self, n: &Cartesian) -> Cartesian { match N { 0 => Cartesian::::default(), 1 => self.vertices[0], - _ => *self - .vertices - .iter() - .max_by(|a, b| { - a.dot(n) - .partial_cmp(&b.dot(n)) - .unwrap_or(std::cmp::Ordering::Equal) - }) - .expect("the 0 match statement should handle empty vectors"), + _ => { + let n_f32 = std::array::from_fn(|i| n[i] as f32); + let vertices_f32 = self + .vertices + .iter() + .map(|v| v.coordinates.map(|x| x as f32)); + let support = vertices_f32 + .max_by(|a, b| { + dot_arrays(a, &n_f32) + .partial_cmp(&dot_arrays(b, &n_f32)) + .unwrap_or(std::cmp::Ordering::Equal) + }) + .expect("the 0 match statement should handle empty vectors"); + support.map(f64::from).into() + } } } } From 8255470ab5ccc3f3351c5be5a0cb18a0d31db520 Mon Sep 17 00:00:00 2001 From: janbridley Date: Thu, 18 Jun 2026 20:45:16 -0400 Subject: [PATCH 02/10] cherry-pick the hyperplatonic commits --- hoomd-geometry/src/shape/convex_polytope.rs | 170 ++++++++++++++++++++ 1 file changed, 170 insertions(+) diff --git a/hoomd-geometry/src/shape/convex_polytope.rs b/hoomd-geometry/src/shape/convex_polytope.rs index 58b0a2e03..462543f8c 100644 --- a/hoomd-geometry/src/shape/convex_polytope.rs +++ b/hoomd-geometry/src/shape/convex_polytope.rs @@ -3,6 +3,8 @@ //! N-Dimensional generalization of a convex polyhedron. +use std::f64::consts::SQRT_2; + use serde::{Deserialize, Serialize}; use crate::{BoundingSphereRadius, Error, SupportMapping}; @@ -205,6 +207,174 @@ impl ConvexPolytope .try_into() .expect("convex polytope should have a positive bounding radius") } + + /// Build the N-dimensional generalization of an octahedron with unit edge length. + /// + /// This shape, also referred to as the hyperoctahedron or cross polytope, is + /// defined by the set of coordinates `{0...±√2/2...0}` for all combinations of + /// `±√2/2` and `N-1` zeros. + /// + /// # Example + /// ``` + /// use approxim::assert_relative_eq; + /// use hoomd_geometry::{ + /// BoundingSphereRadius, Volume, + /// shape::{ConvexPolyhedron, ConvexPolytope}, + /// }; + /// + /// # fn main() -> Result<(), hoomd_geometry::Error> { + /// let diamond = ConvexPolytope::<2>::orthoplex(); + /// assert_relative_eq!(diamond.volume(), 1.0); + /// + /// let octahedron = ConvexPolyhedron::orthoplex(); + /// assert_relative_eq!( + /// octahedron.bounding_sphere_radius().get(), + /// f64::sqrt(2.0) + /// ); + /// + /// // Cross polytopes always have 2*N vertices + /// let octachoron = ConvexPolytope::<4, 8>::orthoplex(); + /// assert_eq!(octachoron.vertices().len(), 8); + /// # Ok(()) + /// # } + /// ``` + /// + /// # Panics + /// If `N=0` or `N > MAX_VERTICES/2`. + #[inline] + #[must_use] + pub fn orthoplex() -> Self { + assert!( + N != 0, + "An orthoplex is not well-defined in zero dimensions!" + ); + let sqrt_2_halves = SQRT_2 / 2.0; + let bounding_radius = SQRT_2.try_into().expect("Hard-coded value"); + let mut vertices = ArrayVec::<_, MAX_VERTICES>::new(); + + for nonzero_index in 0..N { + let coord = Cartesian::::from(std::array::from_fn(|i| { + f64::from(i == nonzero_index) * sqrt_2_halves + })); + vertices.push(coord); + vertices.push(-coord); + } + Self { + vertices, + bounding_radius, + } + } + + /// Build the N-dimensional generalization of a tetrahedron with unit edge length. + /// + /// # Example + /// ``` + /// use approxim::assert_relative_eq; + /// use hoomd_geometry::{ + /// BoundingSphereRadius, Volume, + /// shape::{ConvexPolygon, ConvexPolyhedron, ConvexPolytope}, + /// }; + /// + /// # fn main() -> Result<(), hoomd_geometry::Error> { + /// let triangle = ConvexPolygon::simplex(); + /// assert_relative_eq!(triangle.volume(), f64::sqrt(3.0) / 4.0); + /// + /// let tetrahedron = ConvexPolyhedron::simplex(); + /// assert_relative_eq!( + /// tetrahedron.bounding_sphere_radius().get(), + /// f64::sqrt(3.0 / 8.0) + /// ); + /// + /// // Simplices always have N+1 vertices + /// let pentachoron = ConvexPolytope::<4, 5>::simplex(); + /// assert_eq!(pentachoron.vertices().len(), 5); + /// # Ok(()) + /// # } + /// ``` + /// + /// # Panics + /// If `N+1 > MAX_VERTICES`. + #[inline] + #[must_use] + pub fn simplex() -> Self { + // https://en.wikipedia.org/wiki/Simplex#Cartesian_coordinates_for_a_regular_n-dimensional_simplex_in_Rn + let inv_sqrt2 = SQRT_2.recip(); + + let mut vertices = ArrayVec::<_, MAX_VERTICES>::new(); + + // Un-centered: N vertices at scaled standard basis positions + 1 shared vertex + for k in 0..N { + vertices.push(std::array::from_fn(|i| f64::from(i == k) * inv_sqrt2).into()); + } + + let c = (1.0 - f64::sqrt(N as f64 + 1.0)) / (f64::sqrt(2.0) * N as f64); + vertices.push([c; N].into()); + + // Center by subtracting centroid + let center = Cartesian::from([(inv_sqrt2 + c) / (N as f64 + 1.0); N]); + vertices.iter_mut().for_each(|vertex| *vertex -= center); + + Self { + vertices, + bounding_radius: f64::sqrt(N as f64 / (2.0 * (N as f64 + 1.0))) + .try_into() + .expect("sqrt of positive is positive"), + } + } + + /// Build the N-dimensional generalization of a cube with unit edge length. + /// + /// This shape, also referred to as the hypercube or orthotope, is defined by + /// the signed, length-N combinations of `{±0.5}`. + /// + /// # Example + /// ``` + /// use approxim::assert_relative_eq; + /// use hoomd_geometry::{ + /// BoundingSphereRadius, Volume, + /// shape::{ConvexPolygon, ConvexPolyhedron, ConvexPolytope}, + /// }; + /// + /// # fn main() -> Result<(), hoomd_geometry::Error> { + /// let square = ConvexPolygon::hypercube(); + /// assert_relative_eq!(square.volume(), 1.0); + /// + /// let cube = ConvexPolyhedron::hypercube(); + /// assert_relative_eq!(cube.bounding_sphere_radius().get(), f64::sqrt(3.0)); + /// + /// // Hypercubes have 2^N vertices + /// let hypercube = ConvexPolytope::<4, 16>::hypercube(); + /// assert_eq!(hypercube.vertices().len(), 16); + /// # Ok(()) + /// # } + /// ``` + /// + /// # Panics + /// If `N=0` or `2^N > MAX_VERTICES`. + #[inline] + #[must_use] + pub fn hypercube() -> Self { + #[inline] + fn signed_combinations_of_one_half() -> Vec> { + (0..(1usize << N)) + .map(|bits| { + std::array::from_fn(|i| if bits & (1 << i) == 0 { 0.5 } else { -0.5 }).into() + }) + .collect() + } + assert!( + N != 0, + "Cross polytope is not well-defined in zero dimensions!" + ); + let bounding_radius = f64::sqrt(N as f64) + .try_into() + .expect("sqrt(positive) is positive."); + + Self { + vertices: ArrayVec::<_, MAX_VERTICES>::from_iter(signed_combinations_of_one_half::()), + bounding_radius, + } + } } /// Compute the dot product of two `[f32; N]` arrays. From aefb99194b76382a41fa788d727061967223b33b Mon Sep 17 00:00:00 2001 From: janbridley Date: Thu, 18 Jun 2026 20:45:40 -0400 Subject: [PATCH 03/10] Restore 4d xenocollide, force push is scary --- hoomd-geometry/src/xenocollide.rs | 480 +++++++++++++++++++++++++++++- 1 file changed, 479 insertions(+), 1 deletion(-) diff --git a/hoomd-geometry/src/xenocollide.rs b/hoomd-geometry/src/xenocollide.rs index 0b5d06c6f..05d7d3c57 100644 --- a/hoomd-geometry/src/xenocollide.rs +++ b/hoomd-geometry/src/xenocollide.rs @@ -291,6 +291,214 @@ impl MinkowskiPortalRefinement<3> for Cartesian<3> { } } +/// Find a vector perpendicular to the plane spanned by `a` and `b` in R^4. +/// +/// `counary_cross` requires 3 inputs but during early portal discovery we only have +/// 2 vectors. This function uses a standard basis vector as the 3rd input — at least +/// one of the 4 basis vectors is guaranteed to be non-coplanar with any 2 given vectors. +#[inline] +fn perp_to_2d_subspace(a: Cartesian<4>, b: Cartesian<4>, tol: f64) -> Option> { + for i in 0..4 { + let mut e = [0.0f64; 4]; + e[i] = 1.0; + let n = Cartesian::<4>::counary_cross(&[a, b, e.into()]); + if n.into_iter().any(|x| x.abs() > tol) { + return Some(n); + } + } + None +} + +impl MinkowskiPortalRefinement<4> for Cartesian<4> { + const TOLERANCE: f64 = 2e-12; + + #[inline] + fn outward_normal(portal: &[Cartesian<4>; 4], interior: &Cartesian<4>) -> Cartesian<4> { + let e1 = portal[1] - portal[0]; + let e2 = portal[2] - portal[0]; + let e3 = portal[3] - portal[0]; + let mut n = Self::counary_cross(&[e1, e2, e3]); + if (portal[0] - *interior).dot(&n) < 0.0 { + n = -n; + } + n + } + + fn discover_portal>, B: SupportMapping>>( + s: &MinkowskiDifference<4, A, B>, + v0: &Cartesian<4>, + ) -> Discovery<4> { + // Interior point at origin implies overlap + if v0.into_iter().all(|x| x.abs() < Self::TOLERANCE) { + return Discovery::Known(true); + } + + // v1: first support point + let mut v1 = s.composite_support_mapping(-*v0); + if v1.dot(v0) > 0.0 { + return Discovery::Known(false); + } + + // v2: support in direction perpendicular to the v0-v1 plane + let Some(n) = perp_to_2d_subspace(*v0, v1, Self::TOLERANCE) else { + return Discovery::Known(true); // v0, v1, origin collinear + }; + + let mut v2 = s.composite_support_mapping(n); + if v2.dot(&n) < 0.0 { + return Discovery::Known(false); + } + + // v3: support in direction perpendicular to the v0-v1-v2 subspace + let Some(mut n) = perp_to_2d_subspace(v1 - *v0, v2 - *v0, Self::TOLERANCE) else { + return Discovery::Known(true); + }; + // Orient normal away from v0 (maintain handedness) + if n.dot(v0) > 0.0 { + (v1, v2) = (v2, v1); + n = -n; + } + let v3 = s.composite_support_mapping(n); + if v3.dot(&n) <= 0.0 { + return Discovery::Known(false); + } + + // Recompute portal normal from 3 edge vectors + let mut n = Self::counary_cross(&[v1 - *v0, v2 - *v0, v3 - *v0]); + // Coplanar v1, v2, v3 with v0 produce a zero normal — origin lies on a flat + // boundary region which we treat as overlap + if n.into_iter().all(|x| x.abs() < Self::TOLERANCE) { + return Discovery::Known(true); + } + if n.dot(v0) > 0.0 { + n = -n; + } + + // v4-validation loop + // The 3 face checks use the 4D quadruple product: + // counary_cross(&[v_a, v_b, v4]).dot(&v0) = det([v_a, v_b, v4, v0]) + // Negative means origin is on the wrong side of that face. + let mut v3 = v3; + let mut count = 0_usize; + let v4 = loop { + count += 1; + if count >= XENOCOLLIDE_MAX_ITER { + return Discovery::Known(true); + } + + let v4 = s.composite_support_mapping(n); + if v4.dot(&n) <= 0.0 { + return Discovery::Known(false); + } + + // Face (v1, v2, v4) — opposite v3 + if Self::counary_cross(&[v1, v2, v4]).dot(v0) < 0.0 { + v3 = v4; + n = Self::counary_cross(&[v1 - *v0, v2 - *v0, v3 - *v0]); + if n.dot(v0) > 0.0 { + n = -n; + } + continue; + } + // Face (v2, v3, v4) — opposite v1 + if Self::counary_cross(&[v2, v3, v4]).dot(v0) < 0.0 { + v1 = v4; + n = Self::counary_cross(&[v1 - *v0, v2 - *v0, v3 - *v0]); + if n.dot(v0) > 0.0 { + n = -n; + } + continue; + } + // Face (v1, v4, v3) — opposite v2 + if Self::counary_cross(&[v1, v4, v3]).dot(v0) < 0.0 { + v2 = v4; + n = Self::counary_cross(&[v1 - *v0, v2 - *v0, v3 - *v0]); + if n.dot(v0) > 0.0 { + n = -n; + } + continue; + } + break v4; + }; + + Discovery::Found([v1, v2, v3, v4]) + } + + #[inline] + fn tolerance_check( + portal: &[Cartesian<4>; 4], + v_new: &Cartesian<4>, + normal: &Cartesian<4>, + ) -> Option { + let tolerance = Self::TOLERANCE * normal.norm(); + + let d = (*v_new - portal[0]).dot(normal); + if d.abs() < tolerance { + return Some(false); + } + let d = portal[0].dot(normal); + if d.abs() < tolerance { + return Some(true); + } + None + } + #[inline] + fn narrow_portal(interior: &Cartesian<4>, portal: &mut [Cartesian<4>; 4], v_new: Cartesian<4>) { + // The 4-simplex [portal[0..4] , v_new] has 5 tetrahedral faces. + // The entry face is the current portal (portal[0..4]). + // The 4 exit faces each contain v_new and all portal vertices except portal[i]. + // + // The origin ray enters through the portal and must exit through one of these + // 4 faces. To determine which, we compute the outward-facing normal to each exit + // face and check whether the origin lies on the outward side: + // + // n_i = inward normal of the face opposite portal[i] + // + // The signed distance from a point P to the hyperplane through v_new with normal n + // is (P − v_new) . n. The origin (P = 0) has distance −v_new . n, and interior + // has distance (interior − v_new) . n. We orient n toward portal[i] so that the + // origin is on the opposite (outward) side when v_new . n > 0. + // Among faces the ray crosses (d > 0 and num > 0), we pick the one with smallest + // t = num / (num + d) and fall back to the face with largest d if none qualify. + // + // NOTE: This is fundementally different from 3d, where the cross product + // uniquely identifies the correct facet (plane, in that case) to exit through. + // In 4D, we need to find the *best* facet + let mut exit_face: Option<(usize, f64)> = None; + let mut fallback_i = 0; + let mut fallback_d = f64::NEG_INFINITY; + + for i in 0..4 { + let edges = std::array::from_fn(|k| portal[(i + k + 1) % 4] - v_new); + let n = Self::counary_cross(&edges); + // Skip near-degenerate faces. + if n.norm_squared() < Self::TOLERANCE * Self::TOLERANCE * interior.norm_squared() { + continue; + } + let n = if (portal[i] - v_new).dot(&n) < 0.0 { + -n + } else { + n + }; + let d = v_new.dot(&n); + let num = (*interior - v_new).dot(&n); + + if d > fallback_d { + fallback_d = d; + fallback_i = i; + } + if d > 0.0 && num > 0.0 { + let t = num / (num + d); + match exit_face { + Some((_, best_t)) if t >= best_t => {} + _ => exit_face = Some((i, t)), + } + } + } + portal[exit_face.map_or(fallback_i, |(i, _)| i)] = v_new; + } +} + /// Stateful type that efficiently computes repeated Minkowski differences. pub(crate) struct MinkowskiDifference< 'a, @@ -391,6 +599,14 @@ where let normal = Cartesian::::outward_normal(&portal, &v0); + // If the portal's (N-1)-simplex is degenerate, the normal is near zero and the + // rest of the algorithm does not make sense (so we fail safe -> overlap) + if normal.norm_squared() + < Cartesian::::TOLERANCE * Cartesian::::TOLERANCE * v0.norm_squared() + { + return true; + } + // Hit test: is the origin enclosed by the portal? if portal[0].dot(&normal) >= 0.0 { return true; @@ -444,15 +660,28 @@ where collide::<3, R, A, B>(sa, sb, v_ij, q_ij) } +/// Detect collision between two convex 4D objects via Minkowski Portal Refinement. +#[inline(never)] +pub fn collide4d(sa: &A, sb: &B, v_ij: &Cartesian<4>, q_ij: &R) -> bool +where + A: SupportMapping>, + B: SupportMapping>, + R: Copy, + RotationMatrix<4>: From, +{ + collide::<4, R, A, B>(sa, sb, v_ij, q_ij) +} + #[cfg(test)] mod tests { use super::*; use crate::IntersectsAt; use rstest::*; - use crate::shape::{Circle, Hypercuboid, Hypersphere}; + use crate::shape::{Circle, ConvexPolytope, Hypercuboid, Hypersphere}; use hoomd_utility::valid::PositiveReal; use hoomd_vector::{Angle, Rotation, Versor}; + use rand::{RngExt, SeedableRng, rngs::StdRng}; #[rstest( v => [[0.1, 0.1], [999.9, 0.0], [0.0, 5.123_f64.next_down()], [0.0, 5.123_000_001]], @@ -535,4 +764,253 @@ mod tests { let overlaps = collide3d(&c0, &c1, &v.into(), &theta); assert_eq!(overlaps, c0.intersects_aligned(&c1, &v.into())); } + + #[rstest( + v => [ + [0.1, 0.1, 0.1, 0.1], + [999.9, 0.0, 0.0, -10.9], + [0.0, 5.123, 0.0, 0.0], + [0.0, 0.0, 0.0, 5.123_000_001], + ], + radius => [0.001, 1.0, 4.123], + )] + fn test_4d_spheres_collide(v: [f64; 4], radius: f64) { + let (s0, s1) = ( + Hypersphere { + radius: 1.0.try_into().expect("test value is a positive real"), + }, + Hypersphere::<4> { + radius: radius.try_into().expect("test value is a positive real"), + }, + ); + let q_ij = RotationMatrix::<4>::default(); + + let overlaps = collide4d(&s0, &s1, &v.into(), &q_ij); + + assert_eq!( + overlaps, + s0.intersects_at(&s1, &Cartesian::from(v), &q_ij), + "4D Xenocollide result did not match standard implementation!" + ); + } + + #[rstest( + v => [ + [0.1, 0.1, 0.1, 0.1], + [999.9, 0.0, 0.0, 0.05], + [0.0, 5.123, 0.0, 0.0], + [0.0, 0.0, 0.0, 5.123_000_000_001], + ], + tesseract => [ + [1.0.try_into().expect("test value is a positive real"); 4], + [999.0.try_into().expect("test value is a positive real"), 0.1.try_into().expect("test value is a positive real"), 0.5.try_into().expect("test value is a positive real"), 1.0.try_into().expect("test value is a positive real")], + ], + )] + fn test_tesseracts_collide(v: [f64; 4], tesseract: [PositiveReal; 4]) { + let c0 = Hypercuboid { + edge_lengths: tesseract, + }; + let c1 = Hypercuboid { + edge_lengths: [1.0.try_into().expect("test value is a positive real"); 4], + }; + let q_ij = RotationMatrix::<4>::default(); + + let overlaps = collide4d(&c0, &c1, &v.into(), &q_ij); + assert_eq!(overlaps, c0.intersects_aligned(&c1, &v.into())); + } + + /// Sweep two unit hypercubes from overlapping to separated along a diagonal direction, + /// verifying collide4d matches the analytical result at every step. + #[rstest( + direction => [ + [1.0_f64, 1.0, 1.0, 1.0], + [2.0_f64, 1.0, 1.0, 1.0], + [1.0_f64, 1.0, 1.0, 3.0], + [1.0_f64, 2.0, 3.0, 1.0], + ], + )] + fn test_4d_hypercuboid_diagonal_separation_sweep(direction: [f64; 4]) { + let one: PositiveReal = 1.0.try_into().unwrap(); + let c0 = Hypercuboid { + edge_lengths: [one; 4], + }; + let c1 = Hypercuboid { + edge_lengths: [one; 4], + }; + + // Half-edge-lengths: 0.5 each, sum = 1.0 each. + // Critical t = min_i(1.0 / |dir_i|) + let critical_t = (0..4) + .map(|i| 1.0 / direction[i].abs()) + .reduce(f64::min) + .unwrap(); + + let t_start = critical_t - 0.001; + let t_end = critical_t + 0.001; + let steps = 10_000; + let dt = (t_end - t_start) / f64::from(steps); + + let q_ij = RotationMatrix::<4>::default(); + + for step in 0..=steps { + let t = t_start + dt * f64::from(step); + let d: Cartesian<4> = direction.map(|d_i| d_i * t).into(); + let expected = c0.intersects_aligned(&c1, &d); + let result = collide4d(&c0, &c1, &d, &q_ij); + assert_eq!( + result, expected, + "Mismatch at step {step}, t = {t:.12}, critical_t = {critical_t:.12}" + ); + } + } + + /// Sweep two tesseracts (vertex-based `ConvexPolytope::hypercube()`) along + /// axis-aligned directions, crossing the Minkowski-difference boundary at a cubical + /// facet *center*. + #[rstest( + direction => [ + [1.0_f64, 0.0, 0.0, 0.0], + [0.0_f64, 1.0, 0.0, 0.0], + [0.0_f64, 0.0, 1.0, 0.0], + [0.0_f64, 0.0, 0.0, 1.0], + [1.0_f64, 1.0, 0.0, 0.0], + [0.0_f64, 1.0, 1.0, 0.0], + [0.0_f64, 0.0, 1.0, 1.0], + [2.0_f64, 1.0, 0.0, 0.0], + ], + )] + fn test_4d_hypercube_polytope_axis_separation_sweep(direction: [f64; 4]) { + let c0 = ConvexPolytope::<4, 16>::hypercube(); + let c1 = ConvexPolytope::<4, 16>::hypercube(); + + // Half-edge 0.5 each; collide4d reports overlap (incl. touching) iff |d[i]| <= 1.0 + // for all i. Critical t = min_i(1.0 / |dir_i|); zero components contribute infinity. + let critical_t = (0..4) + .map(|i| 1.0 / direction[i].abs()) + .reduce(f64::min) + .unwrap(); + + let t_start = critical_t - 0.001; + let t_end = critical_t + 0.001; + let steps = 10_000; + let dt = (t_end - t_start) / f64::from(steps); + + let q_ij = RotationMatrix::<4>::default(); + + for step in 0..=steps { + let t = t_start + dt * f64::from(step); + let d: Cartesian<4> = direction.map(|d_i| d_i * t).into(); + let expected = (0..4).all(|i| d[i].abs() <= 1.0); + let result = collide4d(&c0, &c1, &d, &q_ij); + assert_eq!( + result, expected, + "Mismatch at step {step}, t = {t:.12}, critical_t = {critical_t:.12}" + ); + } + } + + /// Stress-test collide4d with two hyperspheres near their overlap boundary. + /// + /// Generates random displacement directions uniformly on S^3 using random + /// unit quaternions ([`Versor`], Muller/Marsaglia Method 19) and random + /// radii in a thin shell of width 1e-3 centered on the overlap radius. + #[rstest( + r0 => [1.0, 0.5, 3.7], + r1 => [1.0, 2.0, 0.8], + seed => [0, 1, 42], + )] + fn test_4d_hypersphere_shell_overlap(r0: f64, r1: f64, seed: u64) { + let boundary = r0 + r1; + let shell_half_width = 5e-10; + let n_samples = 10_000; + + let s0 = Hypersphere::<4> { + radius: r0.try_into().unwrap(), + }; + let s1 = Hypersphere::<4> { + radius: r1.try_into().unwrap(), + }; + let q_ij = RotationMatrix::<4>::default(); + + let mut rng = StdRng::seed_from_u64(seed); + + for i in 0..n_samples { + // Random unit quaternion → uniform direction on S^3 + let q: Versor = rng.random(); + let q = q.get(); + let dir: Cartesian<4> = [q.scalar, q.vector[0], q.vector[1], q.vector[2]].into(); + + let r = boundary - shell_half_width + rng.random::() * 2.0 * shell_half_width; + let d = dir * r; + + let analytical = r < boundary; + let result = collide4d(&s0, &s1, &d, &q_ij); + assert_eq!( + result, analytical, + "Mismatch at sample {i}, |d| = {r:.12}, boundary = {boundary:.12}" + ); + } + } + + /// Construct a 4-simplex from 5 vertices drawn from c·{±1}⁴. + /// + /// 9 of the 10 edges have length 2c√2; the v₃–v₄ edge has length 4c. + fn pentachoron(c: f64) -> ConvexPolytope<4, 8> { + ConvexPolytope::<4, 8>::with_vertices([ + Cartesian::from([c, c, c, c]), + Cartesian::from([c, -c, -c, c]), + Cartesian::from([-c, c, -c, c]), + Cartesian::from([-c, -c, c, c]), + Cartesian::from([c, c, -c, -c]), + ]) + .unwrap() + } + + /// Sweep two simplices along vertex-to-vertex (edge) directions. + /// + /// For two identical simplices with vertices {`v_i`}, a sweep along the + /// edge direction `d0 = v_a − v_b` has an exact analytical boundary: + /// the t* at intersection equals the edge length `|v_a − v_b|` (9 edges have + /// length 2c*sqrt(2), 1 edge has length 4c). + #[rstest( + edge => [ + [0.0_f64, 2.0, 2.0, 0.0], // v0 - v1 + [2.0_f64, 0.0, 2.0, 0.0], // v0 - v2 + [2.0_f64, 2.0, 0.0, 0.0], // v0 - v3 + [0.0_f64, 0.0, 2.0, 2.0], // v0 - v4 + [2.0_f64, -2.0, 0.0, 0.0], // v1 - v2 + [2.0_f64, 0.0, -2.0, 0.0], // v1 - v3 + [0.0_f64, -2.0, 0.0, 2.0], // v1 - v4 + [0.0_f64, 2.0, -2.0, 0.0], // v2 - v3 + [-2.0_f64, 0.0, 0.0, 2.0], // v2 - v4 + [-2.0_f64, -2.0, 2.0, 2.0], // v3 - v4 + ], + c => [1.0, 0.5, 2.0], + )] + fn test_4d_pentachoron_edge_sweep(edge: [f64; 4], c: f64) { + let s0 = pentachoron(c); + let s1 = pentachoron(c); + + let edge: Cartesian<4> = edge.map(|e| e * c).into(); + let edge_length = edge.norm(); + let direction = edge / edge_length; + + let t_start = edge_length - 1.0 / 3.0; + let t_end = edge_length + 0.01; + let steps = 10_000; + let dt = (t_end - t_start) / f64::from(steps); + + let q_ij = RotationMatrix::<4>::default(); + + for step in 0..=steps { + let t = t_start + dt * f64::from(step); + let d = direction * t; + let analytical = t < edge_length; + let result = collide4d(&s0, &s1, &d, &q_ij); + assert_eq!( + result, analytical, + "Mismatch at step {step}, t = {t:.12}, edge_length = {edge_length:.12}" + ); + } + } } From 3ef1f33bf1edab74d574a3181b4ef69ebdb01e11 Mon Sep 17 00:00:00 2001 From: janbridley Date: Thu, 18 Jun 2026 20:45:59 -0400 Subject: [PATCH 04/10] Minor changes to f32 support --- hoomd-geometry/src/shape/convex_polytope.rs | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/hoomd-geometry/src/shape/convex_polytope.rs b/hoomd-geometry/src/shape/convex_polytope.rs index 462543f8c..2767517dc 100644 --- a/hoomd-geometry/src/shape/convex_polytope.rs +++ b/hoomd-geometry/src/shape/convex_polytope.rs @@ -394,18 +394,19 @@ impl SupportMapping> 1 => self.vertices[0], _ => { let n_f32 = std::array::from_fn(|i| n[i] as f32); - let vertices_f32 = self + + let support = self .vertices .iter() - .map(|v| v.coordinates.map(|x| x as f32)); - let support = vertices_f32 - .max_by(|a, b| { - dot_arrays(a, &n_f32) - .partial_cmp(&dot_arrays(b, &n_f32)) - .unwrap_or(std::cmp::Ordering::Equal) + .map(|v| { + let v_f32 = v.coordinates.map(|x| x as f32); + let dot = dot_arrays(&v_f32, &n_f32); + (v_f32, dot) }) - .expect("the 0 match statement should handle empty vectors"); - support.map(f64::from).into() + .max_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal)) + .expect("N > 1 is guaranteed by match statement"); + + support.0.map(f64::from).into() } } } From 662ad78d484dda29c4cfec50357e9424196d23ea Mon Sep 17 00:00:00 2001 From: janbridley Date: Thu, 18 Jun 2026 20:47:09 -0400 Subject: [PATCH 05/10] Revert "Minor changes to f32 support" This reverts commit 3ef1f33bf1edab74d574a3181b4ef69ebdb01e11. --- hoomd-geometry/src/shape/convex_polytope.rs | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/hoomd-geometry/src/shape/convex_polytope.rs b/hoomd-geometry/src/shape/convex_polytope.rs index 2767517dc..462543f8c 100644 --- a/hoomd-geometry/src/shape/convex_polytope.rs +++ b/hoomd-geometry/src/shape/convex_polytope.rs @@ -394,19 +394,18 @@ impl SupportMapping> 1 => self.vertices[0], _ => { let n_f32 = std::array::from_fn(|i| n[i] as f32); - - let support = self + let vertices_f32 = self .vertices .iter() - .map(|v| { - let v_f32 = v.coordinates.map(|x| x as f32); - let dot = dot_arrays(&v_f32, &n_f32); - (v_f32, dot) + .map(|v| v.coordinates.map(|x| x as f32)); + let support = vertices_f32 + .max_by(|a, b| { + dot_arrays(a, &n_f32) + .partial_cmp(&dot_arrays(b, &n_f32)) + .unwrap_or(std::cmp::Ordering::Equal) }) - .max_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal)) - .expect("N > 1 is guaranteed by match statement"); - - support.0.map(f64::from).into() + .expect("the 0 match statement should handle empty vectors"); + support.map(f64::from).into() } } } From f736ee6261ce844c8594dbb064034f409de3f8ec Mon Sep 17 00:00:00 2001 From: janbridley Date: Thu, 18 Jun 2026 20:47:15 -0400 Subject: [PATCH 06/10] Revert "clean f32 support" This reverts commit 7e6df5fa51392f9d6ac28fa84eb1dda0ccabe5c8. --- hoomd-geometry/src/shape/convex_polytope.rs | 31 ++++++--------------- 1 file changed, 9 insertions(+), 22 deletions(-) diff --git a/hoomd-geometry/src/shape/convex_polytope.rs b/hoomd-geometry/src/shape/convex_polytope.rs index 462543f8c..cc7f9de06 100644 --- a/hoomd-geometry/src/shape/convex_polytope.rs +++ b/hoomd-geometry/src/shape/convex_polytope.rs @@ -377,36 +377,23 @@ impl ConvexPolytope } } -/// Compute the dot product of two `[f32; N]` arrays. -#[inline(always)] -fn dot_arrays(l: &[f32; N], r: &[f32; N]) -> f32 { - (0..N).map(|i| l[i] * r[i]).sum() -} - impl SupportMapping> for ConvexPolytope { #[inline] - #[expect(clippy::cast_possible_truncation, reason = "Truncation is ok.")] fn support_mapping(&self, n: &Cartesian) -> Cartesian { match N { 0 => Cartesian::::default(), 1 => self.vertices[0], - _ => { - let n_f32 = std::array::from_fn(|i| n[i] as f32); - let vertices_f32 = self - .vertices - .iter() - .map(|v| v.coordinates.map(|x| x as f32)); - let support = vertices_f32 - .max_by(|a, b| { - dot_arrays(a, &n_f32) - .partial_cmp(&dot_arrays(b, &n_f32)) - .unwrap_or(std::cmp::Ordering::Equal) - }) - .expect("the 0 match statement should handle empty vectors"); - support.map(f64::from).into() - } + _ => *self + .vertices + .iter() + .max_by(|a, b| { + a.dot(n) + .partial_cmp(&b.dot(n)) + .unwrap_or(std::cmp::Ordering::Equal) + }) + .expect("the 0 match statement should handle empty vectors"), } } } From 4e68c2185f320a30f15ed4e696ccac31ea78d6ed Mon Sep 17 00:00:00 2001 From: janbridley Date: Thu, 18 Jun 2026 21:15:10 -0400 Subject: [PATCH 07/10] try restoring old counary code, something is off with rounding --- hoomd-vector/src/cartesian.rs | 70 ++++++++++++++++++++--------------- 1 file changed, 41 insertions(+), 29 deletions(-) diff --git a/hoomd-vector/src/cartesian.rs b/hoomd-vector/src/cartesian.rs index 09d8464b9..f2a1b1d74 100644 --- a/hoomd-vector/src/cartesian.rs +++ b/hoomd-vector/src/cartesian.rs @@ -468,36 +468,48 @@ impl Cartesian<4> { #[inline] #[must_use] pub fn counary_cross(vectors: &[Self; 3]) -> Self { - /// Calculate a 2x2 matrix determinant while compensating for errors. - /// - #[inline(always)] - fn diff_of_products(a: f64, b: f64, c: f64, d: f64) -> f64 { - let cd = c * d; - let err = (-c).mul_add(d, cd); - let dop = a.mul_add(b, -cd); - dop + err - } - let u = &vectors[0]; - let v = &vectors[1]; - let w = &vectors[2]; - - // 2x2 Minors via Kahan's Exact FMA - let m01 = diff_of_products(v[0], w[1], v[1], w[0]); - let m02 = diff_of_products(v[0], w[2], v[2], w[0]); - let m03 = diff_of_products(v[0], w[3], v[3], w[0]); - let m12 = diff_of_products(v[1], w[2], v[2], w[1]); - let m13 = diff_of_products(v[1], w[3], v[3], w[1]); - let m23 = diff_of_products(v[2], w[3], v[3], w[2]); - - // Association is important here, as we can accumulate small sign errors if we - // do not correctly group terms! - [ - (u[1] * m23 - u[2] * m13 + u[3] * m12), - -(u[0] * m23 - u[2] * m03 + u[3] * m02), - (u[0] * m13 - u[1] * m03 + u[3] * m01), - -(u[0] * m12 - u[1] * m02 + u[2] * m01), - ] + std::array::from_fn(|skip| { + Matrix33 { + rows: std::array::from_fn(|r| { + // Skip the column 'skip' by adding 1 to the index when c >= skip + std::array::from_fn(|c| vectors[r][c + usize::from(c >= skip)]) + }), + } + .determinant() + * if skip.is_multiple_of(2) { 1.0 } else { -1.0 } + }) .into() + + // /// Calculate a 2x2 matrix determinant while compensating for errors. + // /// + // #[inline(always)] + // fn diff_of_products(a: f64, b: f64, c: f64, d: f64) -> f64 { + // let cd = c * d; + // let err = (-c).mul_add(d, cd); + // let dop = a.mul_add(b, -cd); + // dop + err + // } + // let u = &vectors[0]; + // let v = &vectors[1]; + // let w = &vectors[2]; + + // // 2x2 Minors via Kahan's Exact FMA + // let m01 = diff_of_products(v[0], w[1], v[1], w[0]); + // let m02 = diff_of_products(v[0], w[2], v[2], w[0]); + // let m03 = diff_of_products(v[0], w[3], v[3], w[0]); + // let m12 = diff_of_products(v[1], w[2], v[2], w[1]); + // let m13 = diff_of_products(v[1], w[3], v[3], w[1]); + // let m23 = diff_of_products(v[2], w[3], v[3], w[2]); + + // // Association is important here, as we can accumulate small sign errors if we + // // do not correctly group terms! + // [ + // (u[1] * m23 - u[2] * m13 + u[3] * m12), + // -(u[0] * m23 - u[2] * m03 + u[3] * m02), + // (u[0] * m13 - u[1] * m03 + u[3] * m01), + // -(u[0] * m12 - u[1] * m02 + u[2] * m01), + // ] + // .into() } } From e29b90816ed0a64b0d34267a99c3041acef437c2 Mon Sep 17 00:00:00 2001 From: janbridley Date: Thu, 18 Jun 2026 21:35:20 -0400 Subject: [PATCH 08/10] restore, wip --- hoomd-vector/src/cartesian.rs | 70 +++++++++++++++-------------------- 1 file changed, 29 insertions(+), 41 deletions(-) diff --git a/hoomd-vector/src/cartesian.rs b/hoomd-vector/src/cartesian.rs index f2a1b1d74..09d8464b9 100644 --- a/hoomd-vector/src/cartesian.rs +++ b/hoomd-vector/src/cartesian.rs @@ -468,48 +468,36 @@ impl Cartesian<4> { #[inline] #[must_use] pub fn counary_cross(vectors: &[Self; 3]) -> Self { - std::array::from_fn(|skip| { - Matrix33 { - rows: std::array::from_fn(|r| { - // Skip the column 'skip' by adding 1 to the index when c >= skip - std::array::from_fn(|c| vectors[r][c + usize::from(c >= skip)]) - }), - } - .determinant() - * if skip.is_multiple_of(2) { 1.0 } else { -1.0 } - }) + /// Calculate a 2x2 matrix determinant while compensating for errors. + /// + #[inline(always)] + fn diff_of_products(a: f64, b: f64, c: f64, d: f64) -> f64 { + let cd = c * d; + let err = (-c).mul_add(d, cd); + let dop = a.mul_add(b, -cd); + dop + err + } + let u = &vectors[0]; + let v = &vectors[1]; + let w = &vectors[2]; + + // 2x2 Minors via Kahan's Exact FMA + let m01 = diff_of_products(v[0], w[1], v[1], w[0]); + let m02 = diff_of_products(v[0], w[2], v[2], w[0]); + let m03 = diff_of_products(v[0], w[3], v[3], w[0]); + let m12 = diff_of_products(v[1], w[2], v[2], w[1]); + let m13 = diff_of_products(v[1], w[3], v[3], w[1]); + let m23 = diff_of_products(v[2], w[3], v[3], w[2]); + + // Association is important here, as we can accumulate small sign errors if we + // do not correctly group terms! + [ + (u[1] * m23 - u[2] * m13 + u[3] * m12), + -(u[0] * m23 - u[2] * m03 + u[3] * m02), + (u[0] * m13 - u[1] * m03 + u[3] * m01), + -(u[0] * m12 - u[1] * m02 + u[2] * m01), + ] .into() - - // /// Calculate a 2x2 matrix determinant while compensating for errors. - // /// - // #[inline(always)] - // fn diff_of_products(a: f64, b: f64, c: f64, d: f64) -> f64 { - // let cd = c * d; - // let err = (-c).mul_add(d, cd); - // let dop = a.mul_add(b, -cd); - // dop + err - // } - // let u = &vectors[0]; - // let v = &vectors[1]; - // let w = &vectors[2]; - - // // 2x2 Minors via Kahan's Exact FMA - // let m01 = diff_of_products(v[0], w[1], v[1], w[0]); - // let m02 = diff_of_products(v[0], w[2], v[2], w[0]); - // let m03 = diff_of_products(v[0], w[3], v[3], w[0]); - // let m12 = diff_of_products(v[1], w[2], v[2], w[1]); - // let m13 = diff_of_products(v[1], w[3], v[3], w[1]); - // let m23 = diff_of_products(v[2], w[3], v[3], w[2]); - - // // Association is important here, as we can accumulate small sign errors if we - // // do not correctly group terms! - // [ - // (u[1] * m23 - u[2] * m13 + u[3] * m12), - // -(u[0] * m23 - u[2] * m03 + u[3] * m02), - // (u[0] * m13 - u[1] * m03 + u[3] * m01), - // -(u[0] * m12 - u[1] * m02 + u[2] * m01), - // ] - // .into() } } From c9dbffebd69cb45e18c870276f4d84af232a2671 Mon Sep 17 00:00:00 2001 From: janbridley Date: Thu, 18 Jun 2026 22:14:20 -0400 Subject: [PATCH 09/10] Tolerance for narrow portal --- hoomd-geometry/src/xenocollide.rs | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/hoomd-geometry/src/xenocollide.rs b/hoomd-geometry/src/xenocollide.rs index 05d7d3c57..faa8f2ca8 100644 --- a/hoomd-geometry/src/xenocollide.rs +++ b/hoomd-geometry/src/xenocollide.rs @@ -391,8 +391,19 @@ impl MinkowskiPortalRefinement<4> for Cartesian<4> { return Discovery::Known(false); } + // Each face test is a 4×4 determinant. For axis-aligned shapes, the + // support vertices are frequently *coplanar* with `v0`, so the determinant + // lands within a few ULP of zero. `< 0.0` then fires on floating-point + // noise, replacing a vertex and bouncing the portal between states until + // `MAX_ITER` bails out with a false-positive overlap. + // + // A determinant within `det_band` of zero means the origin lies ON the + // face (within precision), which we treat as enclosed. Only a clearly- + // negative determinant triggers a vertex replacement. + let det_tol = Self::TOLERANCE * v1.norm() * v2.norm() * v3.norm() * v0.norm(); + // Face (v1, v2, v4) — opposite v3 - if Self::counary_cross(&[v1, v2, v4]).dot(v0) < 0.0 { + if Self::counary_cross(&[v1, v2, v4]).dot(v0) < -det_tol { v3 = v4; n = Self::counary_cross(&[v1 - *v0, v2 - *v0, v3 - *v0]); if n.dot(v0) > 0.0 { @@ -401,7 +412,7 @@ impl MinkowskiPortalRefinement<4> for Cartesian<4> { continue; } // Face (v2, v3, v4) — opposite v1 - if Self::counary_cross(&[v2, v3, v4]).dot(v0) < 0.0 { + if Self::counary_cross(&[v2, v3, v4]).dot(v0) < -det_tol { v1 = v4; n = Self::counary_cross(&[v1 - *v0, v2 - *v0, v3 - *v0]); if n.dot(v0) > 0.0 { @@ -410,7 +421,7 @@ impl MinkowskiPortalRefinement<4> for Cartesian<4> { continue; } // Face (v1, v4, v3) — opposite v2 - if Self::counary_cross(&[v1, v4, v3]).dot(v0) < 0.0 { + if Self::counary_cross(&[v1, v4, v3]).dot(v0) < -det_tol { v2 = v4; n = Self::counary_cross(&[v1 - *v0, v2 - *v0, v3 - *v0]); if n.dot(v0) > 0.0 { From cd9c85c8d577172051e44e973fa123af7f8d87f9 Mon Sep 17 00:00:00 2001 From: janbridley Date: Thu, 18 Jun 2026 22:18:48 -0400 Subject: [PATCH 10/10] more tests and regression --- hoomd-geometry/src/xenocollide.rs | 63 +++++++++++++++++++++++++++++++ 1 file changed, 63 insertions(+) diff --git a/hoomd-geometry/src/xenocollide.rs b/hoomd-geometry/src/xenocollide.rs index faa8f2ca8..dd0d9b85d 100644 --- a/hoomd-geometry/src/xenocollide.rs +++ b/hoomd-geometry/src/xenocollide.rs @@ -875,6 +875,69 @@ mod tests { } } + #[rstest(seed => [0_usize, 1, 2, 7, 42, 2024])] + fn test_4d_hypercuboid_random_near_boundary(seed: usize) { + let one: PositiveReal = 1.0.try_into().unwrap(); + let c0 = Hypercuboid { + edge_lengths: [one; 4], + }; + let c1 = Hypercuboid { + edge_lengths: [one; 4], + }; + let q_ij = RotationMatrix::<4>::default(); + + let mut rng = StdRng::seed_from_u64(seed as u64); + let shell = 1e-9_f64; + let samples = 100_000_usize; + for _ in 0..samples { + // Uniform direction in [-1, 1]^4 (scaling is absorbed by `t0`). + let dir: Cartesian<4> = rng.random(); + let max_abs = dir.into_iter().map(f64::abs).fold(0.0_f64, f64::max); + if max_abs < f64::EPSILON { + continue; + } + let t0 = 1.0 / max_abs; + let t = t0 + (rng.random::() * 2.0 - 1.0) * shell; + let d = dir * t; + let expected = c0.intersects_aligned(&c1, &d); + let result = collide4d(&c0, &c1, &d, &q_ij); + assert_eq!( + result, expected, + "near-boundary mismatch: dir = {dir:?}, t = {t}, d = {d:?}" + ); + } + } + + /// Two unit hypercubes are separated by `d = [2t, t, t, t]`; the analytical boundary + /// is `t = 0.5` (`|d_0| = 1`). Just past it (e.g. `t = 0.5000158`) the + /// shapes do not overlap, but a coplanar portal previously drove MPR into a cycle + /// that falesely reported overlap. + #[test] + fn test_4d_hypercuboid_diagonal_2111_near_boundary() { + let one: PositiveReal = 1.0.try_into().unwrap(); + let c0 = Hypercuboid { + edge_lengths: [one; 4], + }; + let c1 = Hypercuboid { + edge_lengths: [one; 4], + }; + let q_ij = RotationMatrix::<4>::default(); + + for t in [ + 0.5_f64, + 0.5 + 5e-6, + 0.5 + 1.58e-5, // Exact failing case from previous code + 0.5 + 3e-5, + 0.5 + 1e-4, + 0.5 - 1e-4, + ] { + let d: Cartesian<4> = [2.0 * t, t, t, t].into(); + let expected = c0.intersects_aligned(&c1, &d); + let result = collide4d(&c0, &c1, &d, &q_ij); + assert_eq!(result, expected, "t = {t}, d = {d:?}"); + } + } + /// Sweep two tesseracts (vertex-based `ConvexPolytope::hypercube()`) along /// axis-aligned directions, crossing the Minkowski-difference boundary at a cubical /// facet *center*.