From 23adf4af4b81724ea4edfb4994cf37b2c432a5cc Mon Sep 17 00:00:00 2001 From: janbridley Date: Fri, 26 Jun 2026 16:40:44 -0400 Subject: [PATCH 1/8] Generalize GenerateGhosts for Periodic> on N --- .../src/boundary/periodic/cuboid.rs | 785 +++++++++++------- 1 file changed, 473 insertions(+), 312 deletions(-) diff --git a/hoomd-microstate/src/boundary/periodic/cuboid.rs b/hoomd-microstate/src/boundary/periodic/cuboid.rs index 2758ad25e..5cc8e0ffb 100644 --- a/hoomd-microstate/src/boundary/periodic/cuboid.rs +++ b/hoomd-microstate/src/boundary/periodic/cuboid.rs @@ -93,19 +93,46 @@ where } } -impl GenerateGhosts for Periodic> +impl GenerateGhosts for Periodic> where - S: Position> + Copy + Default, + S: Position> + Copy + Default, { #[inline] fn maximum_interaction_range(&self) -> f64 { self.maximum_interaction_range } - /// Place periodic images of sites near the edge of the periodic boundary. + /// Place periodic images of sites near the periodic boundary. /// - /// For 2D cuboids, `generate_ghosts` places ghosts near the 4 edges and 4 - /// vertices. + /// `generate_ghosts` emits one periodic image of `site_properties` for each + /// non-empty combination of the boundary directions that `site_properties` lies + /// within `maximum_interaction_range` of. A site in the bulk produces no ghosts; a + /// site near the middle of an `(N - 1)`-face produces one; a + /// site near an edge (where two boundaries meet) produces three; and so on, up to + /// `2.pow(N) - 1` ghosts for a site near all `2 * N` boundaries at once. + /// The images are emitted in descending order of the set of active directions. + /// Images are currently emitted in descending order of the subset bitmask (with + /// higher bit positions corresponding to folds across a higher-dimensional facet) + /// + /// # Example + /// + /// ``` + /// use hoomd_geometry::shape::Rectangle; + /// use hoomd_microstate::{ + /// boundary::{Periodic, GenerateGhosts}, + /// property::Point, + /// }; + /// use hoomd_vector::Cartesian; + /// + /// # fn main() -> Result<(), Box> { + /// let periodic = + /// Periodic::new(1.0, Rectangle::with_equal_edges(10.0.try_into()?))?; + /// // A site near the right edge produces an image shifted across it. + /// let ghosts = periodic.generate_ghosts(&Point::new(Cartesian::from([4.6, 0.0]))); + /// assert_eq!(ghosts.len(), 1); + /// # Ok(()) + /// # } + /// ``` #[inline] fn generate_ghosts(&self, site_properties: &S) -> ArrayVec { let mut result = ArrayVec::new(); @@ -118,182 +145,142 @@ where return result; } - let new_site = |x, y| { - let mut new_site = *site_properties; - new_site.position_mut()[0] += x * self.shape.edge_lengths[0].get(); - new_site.position_mut()[1] += y * self.shape.edge_lengths[1].get(); - new_site - }; - - let near_left = r[0] < min[0] + self.maximum_interaction_range; - let near_right = r[0] > max[0] - self.maximum_interaction_range; - let near_top = r[1] > max[1] - self.maximum_interaction_range; - let near_bottom = r[1] < min[1] + self.maximum_interaction_range; - - if near_right { - result.push(new_site(-1.0, 0.0)); - } - if near_left { - result.push(new_site(1.0, 0.0)); - } - if near_top { - result.push(new_site(0.0, -1.0)); - } - if near_bottom { - result.push(new_site(0.0, 1.0)); - } - if near_right && near_top { - result.push(new_site(-1.0, -1.0)); - } - if near_right && near_bottom { - result.push(new_site(-1.0, 1.0)); - } - if near_left && near_top { - result.push(new_site(1.0, -1.0)); + // Record each direction the site is near, and the translation that folds a + // ghost back across that boundary. `near_mask` is the set of active directions + // and `dim_offset` is the per-direction periodic translation. + let mut near_mask = 0u32; + let mut dim_offset = [0.0_f64; N]; + for i in 0..N { + if r[i] > max[i] - self.maximum_interaction_range { + near_mask |= 1 << i; + dim_offset[i] = -self.shape.edge_lengths[i].get(); + } else if r[i] < min[i] + self.maximum_interaction_range { + near_mask |= 1 << i; + dim_offset[i] = self.shape.edge_lengths[i].get(); + } } - if near_left && near_bottom { - result.push(new_site(1.0, 1.0)); + + // Emit one image for every non-empty subset of the active directions by + // walking the subsets of `near_mask` in descending order. + let mut subset = near_mask; + while subset != 0 { + let mut ghost = *site_properties; + let pos = ghost.position_mut(); + for i in 0..N { + if subset & (1 << i) != 0 { + pos[i] += dim_offset[i]; + } + } + result.push(ghost); + subset = (subset - 1) & near_mask; } result } } -impl GenerateGhosts for Periodic> -where - S: Position> + Copy + Default, -{ - #[inline] - fn maximum_interaction_range(&self) -> f64 { - self.maximum_interaction_range - } +#[cfg(test)] +mod tests { + use super::*; + use crate::property::Point; - /// Place periodic images of sites near the edge of the periodic boundary. - /// - /// For 3D cuboids, `generate_ghosts` places ghosts near the 6 faces, 12 edges, - /// and 8 vertices. - #[inline] - fn generate_ghosts(&self, site_properties: &S) -> ArrayVec { - let mut result = ArrayVec::new(); + use approxim::assert_relative_eq; + use rand::{SeedableRng, distr::Distribution, rngs::StdRng}; - let r = site_properties.position(); - let max = self.shape.maximal_extents(); - let min = self.shape.minimal_extents(); + const N_SAMPLES: usize = 1024; - if !self.shape.is_point_inside(r) { - return result; + /// Assert the ghosts equal the expected positions, compared as sets. + /// + /// `generate_ghosts` is permitted to emit its periodic images in any order, + /// so tests compare the set of ghost positions rather than a specific + /// ordering. + fn assert_ghost_positions( + ghosts: &ArrayVec>, MAX_GHOSTS>, + expected: &[[f64; N]], + ) { + assert_eq!(ghosts.len(), expected.len()); + let mut got: Vec> = ghosts.iter().map(|ghost| ghost.position).collect(); + let mut want: Vec> = expected + .iter() + .map(|coords| Cartesian::from(*coords)) + .collect(); + got.sort_by(|a, b| a.coordinates.partial_cmp(&b.coordinates).unwrap()); + want.sort_by(|a, b| a.coordinates.partial_cmp(&b.coordinates).unwrap()); + for (got, want) in got.iter().zip(want.iter()) { + assert_relative_eq!(got, want); } + } - let new_site = |x, y, z| { - let mut new_site = *site_properties; - new_site.position_mut()[0] += x * self.shape.edge_lengths[0].get(); - new_site.position_mut()[1] += y * self.shape.edge_lengths[1].get(); - new_site.position_mut()[2] += z * self.shape.edge_lengths[2].get(); - new_site - }; - - let near_left = r[0] < min[0] + self.maximum_interaction_range; - let near_right = r[0] > max[0] - self.maximum_interaction_range; - let near_top = r[1] > max[1] - self.maximum_interaction_range; - let near_bottom = r[1] < min[1] + self.maximum_interaction_range; - let near_front = r[2] > max[2] - self.maximum_interaction_range; - let near_back = r[2] < min[2] + self.maximum_interaction_range; - - if near_right { - result.push(new_site(-1.0, 0.0, 0.0)); - } - if near_left { - result.push(new_site(1.0, 0.0, 0.0)); - } - if near_top { - result.push(new_site(0.0, -1.0, 0.0)); - } - if near_bottom { - result.push(new_site(0.0, 1.0, 0.0)); - } - if near_front { - result.push(new_site(0.0, 0.0, -1.0)); - } - if near_back { - result.push(new_site(0.0, 0.0, 1.0)); - } + mod cuboid_1 { + use super::*; - if near_right && near_top { - result.push(new_site(-1.0, -1.0, 0.0)); - } - if near_right && near_bottom { - result.push(new_site(-1.0, 1.0, 0.0)); - } - if near_right && near_front { - result.push(new_site(-1.0, 0.0, -1.0)); - } - if near_right && near_back { - result.push(new_site(-1.0, 0.0, 1.0)); - } - if near_left && near_top { - result.push(new_site(1.0, -1.0, 0.0)); - } - if near_left && near_bottom { - result.push(new_site(1.0, 1.0, 0.0)); - } - if near_left && near_front { - result.push(new_site(1.0, 0.0, -1.0)); - } - if near_left && near_back { - result.push(new_site(1.0, 0.0, 1.0)); + fn pos(value: f64) -> PositiveReal { + value + .try_into() + .expect("hard-coded constant should be positive") } - if near_top && near_front { - result.push(new_site(0.0, -1.0, -1.0)); - } - if near_bottom && near_front { - result.push(new_site(0.0, 1.0, -1.0)); - } - if near_top && near_back { - result.push(new_site(0.0, -1.0, 1.0)); - } - if near_bottom && near_back { - result.push(new_site(0.0, 1.0, 1.0)); + #[test] + fn maximum_allowable() { + let cuboid = Hypercuboid { + edge_lengths: [pos(10.0)], + }; + assert_eq!(cuboid.maximum_allowable_interaction_range(), 5.0); } - if near_right && near_top && near_front { - result.push(new_site(-1.0, -1.0, -1.0)); - } - if near_right && near_top && near_back { - result.push(new_site(-1.0, -1.0, 1.0)); - } - if near_right && near_bottom && near_front { - result.push(new_site(-1.0, 1.0, -1.0)); - } - if near_right && near_bottom && near_back { - result.push(new_site(-1.0, 1.0, 1.0)); - } - if near_left && near_top && near_front { - result.push(new_site(1.0, -1.0, -1.0)); - } - if near_left && near_top && near_back { - result.push(new_site(1.0, -1.0, 1.0)); - } - if near_left && near_bottom && near_front { - result.push(new_site(1.0, 1.0, -1.0)); - } - if near_left && near_bottom && near_back { - result.push(new_site(1.0, 1.0, 1.0)); + #[test] + fn wrap() { + let cuboid = Hypercuboid { + edge_lengths: [pos(20.0)], + }; + let periodic = Periodic::new(0.0, cuboid).expect("hard-coded range should be valid"); + + let point = Point::new([5.0].into()); + assert_eq!(periodic.wrap(point), Ok(point)); + + let point = Point::new([10.0].into()); + assert_eq!(periodic.wrap(point), Ok(Point::new([-10.0].into()))); + + let point = Point::new([25.0].into()); + assert_eq!(periodic.wrap(point), Ok(Point::new([5.0].into()))); } - result - } -} + #[test] + fn no_ghosts() { + let cuboid = Hypercuboid { + edge_lengths: [pos(20.0)], + }; + let periodic = Periodic::new(1.0, cuboid).expect("hard-coded range should be valid"); -#[cfg(test)] -mod tests { - use super::*; - use crate::property::Point; + // a site in the bulk produces no ghosts + assert!( + periodic + .generate_ghosts(&Point::new([0.0].into())) + .is_empty() + ); + // a site outside the boundary produces no ghosts + assert!( + periodic + .generate_ghosts(&Point::new([10.5].into())) + .is_empty() + ); + } - use approxim::assert_relative_eq; - use rand::{SeedableRng, distr::Distribution, rngs::StdRng}; + #[test] + fn ghosts() { + let cuboid = Hypercuboid { + edge_lengths: [pos(20.0)], + }; + let periodic = Periodic::new(1.0, cuboid).expect("hard-coded range should be valid"); - const N_SAMPLES: usize = 1024; + // a site near each end produces one image shifted across that boundary + let ghosts = periodic.generate_ghosts(&Point::new([9.5].into())); + assert_ghost_positions(&ghosts, &[[-10.5]]); + + let ghosts = periodic.generate_ghosts(&Point::new([-9.5].into())); + assert_ghost_positions(&ghosts, &[[10.5]]); + } + } mod cuboid_2 { use super::*; @@ -419,47 +406,31 @@ mod tests { let ghosts = periodic.generate_ghosts(&Point::new([0.0, 5.5].into())); assert!(ghosts.is_empty()); - // edges + // faces (one ghost each) let ghosts = periodic.generate_ghosts(&Point::new([9.5, 0.0].into())); - assert_eq!(ghosts.len(), 1); - assert_relative_eq!(ghosts[0].position, [-10.5, 0.0].into()); + assert_ghost_positions(&ghosts, &[[-10.5, 0.0]]); let ghosts = periodic.generate_ghosts(&Point::new([-9.5, 0.0].into())); - assert_eq!(ghosts.len(), 1); - assert_relative_eq!(ghosts[0].position, [10.5, 0.0].into()); + assert_ghost_positions(&ghosts, &[[10.5, 0.0]]); let ghosts = periodic.generate_ghosts(&Point::new([0.0, 4.5].into())); - assert_eq!(ghosts.len(), 1); - assert_relative_eq!(ghosts[0].position, [0.0, -5.5].into()); + assert_ghost_positions(&ghosts, &[[0.0, -5.5]]); let ghosts = periodic.generate_ghosts(&Point::new([0.0, -4.5].into())); - assert_eq!(ghosts.len(), 1); - assert_relative_eq!(ghosts[0].position, [0.0, 5.5].into()); + assert_ghost_positions(&ghosts, &[[0.0, 5.5]]); - // vertices + // vertices (three ghosts each) let ghosts = periodic.generate_ghosts(&Point::new([9.5, 4.5].into())); - assert_eq!(ghosts.len(), 3); - assert_relative_eq!(ghosts[0].position, [-10.5, 4.5].into()); - assert_relative_eq!(ghosts[1].position, [9.5, -5.5].into()); - assert_relative_eq!(ghosts[2].position, [-10.5, -5.5].into()); + assert_ghost_positions(&ghosts, &[[-10.5, 4.5], [9.5, -5.5], [-10.5, -5.5]]); let ghosts = periodic.generate_ghosts(&Point::new([9.5, -4.5].into())); - assert_eq!(ghosts.len(), 3); - assert_relative_eq!(ghosts[0].position, [-10.5, -4.5].into()); - assert_relative_eq!(ghosts[1].position, [9.5, 5.5].into()); - assert_relative_eq!(ghosts[2].position, [-10.5, 5.5].into()); + assert_ghost_positions(&ghosts, &[[-10.5, -4.5], [9.5, 5.5], [-10.5, 5.5]]); let ghosts = periodic.generate_ghosts(&Point::new([-9.5, 4.5].into())); - assert_eq!(ghosts.len(), 3); - assert_relative_eq!(ghosts[0].position, [10.5, 4.5].into()); - assert_relative_eq!(ghosts[1].position, [-9.5, -5.5].into()); - assert_relative_eq!(ghosts[2].position, [10.5, -5.5].into()); + assert_ghost_positions(&ghosts, &[[10.5, 4.5], [-9.5, -5.5], [10.5, -5.5]]); let ghosts = periodic.generate_ghosts(&Point::new([-9.5, -4.5].into())); - assert_eq!(ghosts.len(), 3); - assert_relative_eq!(ghosts[0].position, [10.5, -4.5].into()); - assert_relative_eq!(ghosts[1].position, [-9.5, 5.5].into()); - assert_relative_eq!(ghosts[2].position, [10.5, 5.5].into()); + assert_ghost_positions(&ghosts, &[[10.5, -4.5], [-9.5, 5.5], [10.5, 5.5]]); } } @@ -617,184 +588,374 @@ mod tests { let ghosts = periodic.generate_ghosts(&Point::new([0.0, 5.5, 0.0].into())); assert!(ghosts.is_empty()); - // faces + // faces (one ghost each) let ghosts = periodic.generate_ghosts(&Point::new([9.5, 0.0, 0.0].into())); - assert_eq!(ghosts.len(), 1); - assert_relative_eq!(ghosts[0].position, [-10.5, 0.0, 0.0].into()); + assert_ghost_positions(&ghosts, &[[-10.5, 0.0, 0.0]]); let ghosts = periodic.generate_ghosts(&Point::new([-9.5, 0.0, 0.0].into())); - assert_eq!(ghosts.len(), 1); - assert_relative_eq!(ghosts[0].position, [10.5, 0.0, 0.0].into()); + assert_ghost_positions(&ghosts, &[[10.5, 0.0, 0.0]]); let ghosts = periodic.generate_ghosts(&Point::new([0.0, 4.5, 0.0].into())); - assert_eq!(ghosts.len(), 1); - assert_relative_eq!(ghosts[0].position, [0.0, -5.5, 0.0].into()); + assert_ghost_positions(&ghosts, &[[0.0, -5.5, 0.0]]); let ghosts = periodic.generate_ghosts(&Point::new([0.0, -4.5, 0.0].into())); - assert_eq!(ghosts.len(), 1); - assert_relative_eq!(ghosts[0].position, [0.0, 5.5, 0.0].into()); + assert_ghost_positions(&ghosts, &[[0.0, 5.5, 0.0]]); let ghosts = periodic.generate_ghosts(&Point::new([0.0, 0.0, 19.5].into())); - assert_eq!(ghosts.len(), 1); - assert_relative_eq!(ghosts[0].position, [0.0, 0.0, -20.5].into()); + assert_ghost_positions(&ghosts, &[[0.0, 0.0, -20.5]]); let ghosts = periodic.generate_ghosts(&Point::new([0.0, 0.0, -19.5].into())); - assert_eq!(ghosts.len(), 1); - assert_relative_eq!(ghosts[0].position, [0.0, 0.0, 20.5].into()); + assert_ghost_positions(&ghosts, &[[0.0, 0.0, 20.5]]); - // edges + // edges (three ghosts each) let ghosts = periodic.generate_ghosts(&Point::new([9.5, 4.5, 0.0].into())); - assert_eq!(ghosts.len(), 3); - assert_relative_eq!(ghosts[0].position, [-10.5, 4.5, 0.0].into()); - assert_relative_eq!(ghosts[1].position, [9.5, -5.5, 0.0].into()); - assert_relative_eq!(ghosts[2].position, [-10.5, -5.5, 0.0].into()); + assert_ghost_positions( + &ghosts, + &[[-10.5, 4.5, 0.0], [9.5, -5.5, 0.0], [-10.5, -5.5, 0.0]], + ); let ghosts = periodic.generate_ghosts(&Point::new([9.5, -4.5, 0.0].into())); - assert_eq!(ghosts.len(), 3); - assert_relative_eq!(ghosts[0].position, [-10.5, -4.5, 0.0].into()); - assert_relative_eq!(ghosts[1].position, [9.5, 5.5, 0.0].into()); - assert_relative_eq!(ghosts[2].position, [-10.5, 5.5, 0.0].into()); + assert_ghost_positions( + &ghosts, + &[[-10.5, -4.5, 0.0], [9.5, 5.5, 0.0], [-10.5, 5.5, 0.0]], + ); let ghosts = periodic.generate_ghosts(&Point::new([-9.5, 4.5, 0.0].into())); - assert_eq!(ghosts.len(), 3); - assert_relative_eq!(ghosts[0].position, [10.5, 4.5, 0.0].into()); - assert_relative_eq!(ghosts[1].position, [-9.5, -5.5, 0.0].into()); - assert_relative_eq!(ghosts[2].position, [10.5, -5.5, 0.0].into()); + assert_ghost_positions( + &ghosts, + &[[10.5, 4.5, 0.0], [-9.5, -5.5, 0.0], [10.5, -5.5, 0.0]], + ); let ghosts = periodic.generate_ghosts(&Point::new([-9.5, -4.5, 0.0].into())); - assert_eq!(ghosts.len(), 3); - assert_relative_eq!(ghosts[0].position, [10.5, -4.5, 0.0].into()); - assert_relative_eq!(ghosts[1].position, [-9.5, 5.5, 0.0].into()); - assert_relative_eq!(ghosts[2].position, [10.5, 5.5, 0.0].into()); + assert_ghost_positions( + &ghosts, + &[[10.5, -4.5, 0.0], [-9.5, 5.5, 0.0], [10.5, 5.5, 0.0]], + ); let ghosts = periodic.generate_ghosts(&Point::new([9.5, 0.0, 19.5].into())); - assert_eq!(ghosts.len(), 3); - assert_relative_eq!(ghosts[0].position, [-10.5, 0.0, 19.5].into()); - assert_relative_eq!(ghosts[1].position, [9.5, 0.0, -20.5].into()); - assert_relative_eq!(ghosts[2].position, [-10.5, 0.0, -20.5].into()); + assert_ghost_positions( + &ghosts, + &[[-10.5, 0.0, 19.5], [9.5, 0.0, -20.5], [-10.5, 0.0, -20.5]], + ); let ghosts = periodic.generate_ghosts(&Point::new([-9.5, 0.0, 19.5].into())); - assert_eq!(ghosts.len(), 3); - assert_relative_eq!(ghosts[0].position, [10.5, 0.0, 19.5].into()); - assert_relative_eq!(ghosts[1].position, [-9.5, 0.0, -20.5].into()); - assert_relative_eq!(ghosts[2].position, [10.5, 0.0, -20.5].into()); + assert_ghost_positions( + &ghosts, + &[[10.5, 0.0, 19.5], [-9.5, 0.0, -20.5], [10.5, 0.0, -20.5]], + ); let ghosts = periodic.generate_ghosts(&Point::new([0.0, 4.5, 19.5].into())); - assert_eq!(ghosts.len(), 3); - assert_relative_eq!(ghosts[0].position, [0.0, -5.5, 19.5].into()); - assert_relative_eq!(ghosts[1].position, [0.0, 4.5, -20.5].into()); - assert_relative_eq!(ghosts[2].position, [0.0, -5.5, -20.5].into()); + assert_ghost_positions( + &ghosts, + &[[0.0, -5.5, 19.5], [0.0, 4.5, -20.5], [0.0, -5.5, -20.5]], + ); let ghosts = periodic.generate_ghosts(&Point::new([0.0, -4.5, 19.5].into())); - assert_eq!(ghosts.len(), 3); - assert_relative_eq!(ghosts[0].position, [0.0, 5.5, 19.5].into()); - assert_relative_eq!(ghosts[1].position, [0.0, -4.5, -20.5].into()); - assert_relative_eq!(ghosts[2].position, [0.0, 5.5, -20.5].into()); + assert_ghost_positions( + &ghosts, + &[[0.0, 5.5, 19.5], [0.0, -4.5, -20.5], [0.0, 5.5, -20.5]], + ); let ghosts = periodic.generate_ghosts(&Point::new([9.5, 0.0, -19.5].into())); - assert_eq!(ghosts.len(), 3); - assert_relative_eq!(ghosts[0].position, [-10.5, 0.0, -19.5].into()); - assert_relative_eq!(ghosts[1].position, [9.5, 0.0, 20.5].into()); - assert_relative_eq!(ghosts[2].position, [-10.5, 0.0, 20.5].into()); + assert_ghost_positions( + &ghosts, + &[[-10.5, 0.0, -19.5], [9.5, 0.0, 20.5], [-10.5, 0.0, 20.5]], + ); let ghosts = periodic.generate_ghosts(&Point::new([-9.5, 0.0, -19.5].into())); - assert_eq!(ghosts.len(), 3); - assert_relative_eq!(ghosts[0].position, [10.5, 0.0, -19.5].into()); - assert_relative_eq!(ghosts[1].position, [-9.5, 0.0, 20.5].into()); - assert_relative_eq!(ghosts[2].position, [10.5, 0.0, 20.5].into()); + assert_ghost_positions( + &ghosts, + &[[10.5, 0.0, -19.5], [-9.5, 0.0, 20.5], [10.5, 0.0, 20.5]], + ); let ghosts = periodic.generate_ghosts(&Point::new([0.0, 4.5, -19.5].into())); - assert_eq!(ghosts.len(), 3); - assert_relative_eq!(ghosts[0].position, [0.0, -5.5, -19.5].into()); - assert_relative_eq!(ghosts[1].position, [0.0, 4.5, 20.5].into()); - assert_relative_eq!(ghosts[2].position, [0.0, -5.5, 20.5].into()); + assert_ghost_positions( + &ghosts, + &[[0.0, -5.5, -19.5], [0.0, 4.5, 20.5], [0.0, -5.5, 20.5]], + ); let ghosts = periodic.generate_ghosts(&Point::new([0.0, -4.5, -19.5].into())); - assert_eq!(ghosts.len(), 3); - assert_relative_eq!(ghosts[0].position, [0.0, 5.5, -19.5].into()); - assert_relative_eq!(ghosts[1].position, [0.0, -4.5, 20.5].into()); - assert_relative_eq!(ghosts[2].position, [0.0, 5.5, 20.5].into()); + assert_ghost_positions( + &ghosts, + &[[0.0, 5.5, -19.5], [0.0, -4.5, 20.5], [0.0, 5.5, 20.5]], + ); - // vertices + // vertices (seven ghosts each) let ghosts = periodic.generate_ghosts(&Point::new([9.5, 4.5, 19.5].into())); - assert_eq!(ghosts.len(), 7); - assert_relative_eq!(ghosts[0].position, [-10.5, 4.5, 19.5].into()); - assert_relative_eq!(ghosts[1].position, [9.5, -5.5, 19.5].into()); - assert_relative_eq!(ghosts[2].position, [9.5, 4.5, -20.5].into()); - assert_relative_eq!(ghosts[3].position, [-10.5, -5.5, 19.5].into()); - assert_relative_eq!(ghosts[4].position, [-10.5, 4.5, -20.5].into()); - assert_relative_eq!(ghosts[5].position, [9.5, -5.5, -20.5].into()); - assert_relative_eq!(ghosts[6].position, [-10.5, -5.5, -20.5].into()); + assert_ghost_positions( + &ghosts, + &[ + [-10.5, 4.5, 19.5], + [9.5, -5.5, 19.5], + [9.5, 4.5, -20.5], + [-10.5, -5.5, 19.5], + [-10.5, 4.5, -20.5], + [9.5, -5.5, -20.5], + [-10.5, -5.5, -20.5], + ], + ); let ghosts = periodic.generate_ghosts(&Point::new([9.5, 4.5, -19.5].into())); - assert_eq!(ghosts.len(), 7); - assert_relative_eq!(ghosts[0].position, [-10.5, 4.5, -19.5].into()); - assert_relative_eq!(ghosts[1].position, [9.5, -5.5, -19.5].into()); - assert_relative_eq!(ghosts[2].position, [9.5, 4.5, 20.5].into()); - assert_relative_eq!(ghosts[3].position, [-10.5, -5.5, -19.5].into()); - assert_relative_eq!(ghosts[4].position, [-10.5, 4.5, 20.5].into()); - assert_relative_eq!(ghosts[5].position, [9.5, -5.5, 20.5].into()); - assert_relative_eq!(ghosts[6].position, [-10.5, -5.5, 20.5].into()); + assert_ghost_positions( + &ghosts, + &[ + [-10.5, 4.5, -19.5], + [9.5, -5.5, -19.5], + [9.5, 4.5, 20.5], + [-10.5, -5.5, -19.5], + [-10.5, 4.5, 20.5], + [9.5, -5.5, 20.5], + [-10.5, -5.5, 20.5], + ], + ); let ghosts = periodic.generate_ghosts(&Point::new([9.5, -4.5, 19.5].into())); - assert_eq!(ghosts.len(), 7); - assert_relative_eq!(ghosts[0].position, [-10.5, -4.5, 19.5].into()); - assert_relative_eq!(ghosts[1].position, [9.5, 5.5, 19.5].into()); - assert_relative_eq!(ghosts[2].position, [9.5, -4.5, -20.5].into()); - assert_relative_eq!(ghosts[3].position, [-10.5, 5.5, 19.5].into()); - assert_relative_eq!(ghosts[4].position, [-10.5, -4.5, -20.5].into()); - assert_relative_eq!(ghosts[5].position, [9.5, 5.5, -20.5].into()); - assert_relative_eq!(ghosts[6].position, [-10.5, 5.5, -20.5].into()); + assert_ghost_positions( + &ghosts, + &[ + [-10.5, -4.5, 19.5], + [9.5, 5.5, 19.5], + [9.5, -4.5, -20.5], + [-10.5, 5.5, 19.5], + [-10.5, -4.5, -20.5], + [9.5, 5.5, -20.5], + [-10.5, 5.5, -20.5], + ], + ); let ghosts = periodic.generate_ghosts(&Point::new([9.5, -4.5, -19.5].into())); - assert_eq!(ghosts.len(), 7); - assert_relative_eq!(ghosts[0].position, [-10.5, -4.5, -19.5].into()); - assert_relative_eq!(ghosts[1].position, [9.5, 5.5, -19.5].into()); - assert_relative_eq!(ghosts[2].position, [9.5, -4.5, 20.5].into()); - assert_relative_eq!(ghosts[3].position, [-10.5, 5.5, -19.5].into()); - assert_relative_eq!(ghosts[4].position, [-10.5, -4.5, 20.5].into()); - assert_relative_eq!(ghosts[5].position, [9.5, 5.5, 20.5].into()); - assert_relative_eq!(ghosts[6].position, [-10.5, 5.5, 20.5].into()); + assert_ghost_positions( + &ghosts, + &[ + [-10.5, -4.5, -19.5], + [9.5, 5.5, -19.5], + [9.5, -4.5, 20.5], + [-10.5, 5.5, -19.5], + [-10.5, -4.5, 20.5], + [9.5, 5.5, 20.5], + [-10.5, 5.5, 20.5], + ], + ); let ghosts = periodic.generate_ghosts(&Point::new([-9.5, 4.5, 19.5].into())); - assert_eq!(ghosts.len(), 7); - assert_relative_eq!(ghosts[0].position, [10.5, 4.5, 19.5].into()); - assert_relative_eq!(ghosts[1].position, [-9.5, -5.5, 19.5].into()); - assert_relative_eq!(ghosts[2].position, [-9.5, 4.5, -20.5].into()); - assert_relative_eq!(ghosts[3].position, [10.5, -5.5, 19.5].into()); - assert_relative_eq!(ghosts[4].position, [10.5, 4.5, -20.5].into()); - assert_relative_eq!(ghosts[5].position, [-9.5, -5.5, -20.5].into()); - assert_relative_eq!(ghosts[6].position, [10.5, -5.5, -20.5].into()); + assert_ghost_positions( + &ghosts, + &[ + [10.5, 4.5, 19.5], + [-9.5, -5.5, 19.5], + [-9.5, 4.5, -20.5], + [10.5, -5.5, 19.5], + [10.5, 4.5, -20.5], + [-9.5, -5.5, -20.5], + [10.5, -5.5, -20.5], + ], + ); let ghosts = periodic.generate_ghosts(&Point::new([-9.5, 4.5, -19.5].into())); - assert_eq!(ghosts.len(), 7); - assert_relative_eq!(ghosts[0].position, [10.5, 4.5, -19.5].into()); - assert_relative_eq!(ghosts[1].position, [-9.5, -5.5, -19.5].into()); - assert_relative_eq!(ghosts[2].position, [-9.5, 4.5, 20.5].into()); - assert_relative_eq!(ghosts[3].position, [10.5, -5.5, -19.5].into()); - assert_relative_eq!(ghosts[4].position, [10.5, 4.5, 20.5].into()); - assert_relative_eq!(ghosts[5].position, [-9.5, -5.5, 20.5].into()); - assert_relative_eq!(ghosts[6].position, [10.5, -5.5, 20.5].into()); + assert_ghost_positions( + &ghosts, + &[ + [10.5, 4.5, -19.5], + [-9.5, -5.5, -19.5], + [-9.5, 4.5, 20.5], + [10.5, -5.5, -19.5], + [10.5, 4.5, 20.5], + [-9.5, -5.5, 20.5], + [10.5, -5.5, 20.5], + ], + ); let ghosts = periodic.generate_ghosts(&Point::new([-9.5, -4.5, 19.5].into())); - assert_eq!(ghosts.len(), 7); - assert_relative_eq!(ghosts[0].position, [10.5, -4.5, 19.5].into()); - assert_relative_eq!(ghosts[1].position, [-9.5, 5.5, 19.5].into()); - assert_relative_eq!(ghosts[2].position, [-9.5, -4.5, -20.5].into()); - assert_relative_eq!(ghosts[3].position, [10.5, 5.5, 19.5].into()); - assert_relative_eq!(ghosts[4].position, [10.5, -4.5, -20.5].into()); - assert_relative_eq!(ghosts[5].position, [-9.5, 5.5, -20.5].into()); - assert_relative_eq!(ghosts[6].position, [10.5, 5.5, -20.5].into()); + assert_ghost_positions( + &ghosts, + &[ + [10.5, -4.5, 19.5], + [-9.5, 5.5, 19.5], + [-9.5, -4.5, -20.5], + [10.5, 5.5, 19.5], + [10.5, -4.5, -20.5], + [-9.5, 5.5, -20.5], + [10.5, 5.5, -20.5], + ], + ); let ghosts = periodic.generate_ghosts(&Point::new([-9.5, -4.5, -19.5].into())); - assert_eq!(ghosts.len(), 7); - assert_relative_eq!(ghosts[0].position, [10.5, -4.5, -19.5].into()); - assert_relative_eq!(ghosts[1].position, [-9.5, 5.5, -19.5].into()); - assert_relative_eq!(ghosts[2].position, [-9.5, -4.5, 20.5].into()); - assert_relative_eq!(ghosts[3].position, [10.5, 5.5, -19.5].into()); - assert_relative_eq!(ghosts[4].position, [10.5, -4.5, 20.5].into()); - assert_relative_eq!(ghosts[5].position, [-9.5, 5.5, 20.5].into()); - assert_relative_eq!(ghosts[6].position, [10.5, 5.5, 20.5].into()); + assert_ghost_positions( + &ghosts, + &[ + [10.5, -4.5, -19.5], + [-9.5, 5.5, -19.5], + [-9.5, -4.5, 20.5], + [10.5, 5.5, -19.5], + [10.5, -4.5, 20.5], + [-9.5, 5.5, 20.5], + [10.5, 5.5, 20.5], + ], + ); + } + } + + mod cuboid_4 { + use super::*; + use rstest::rstest; + + fn pos(value: f64) -> PositiveReal { + value + .try_into() + .expect("hard-coded constant should be positive") + } + + #[rstest] + #[case::min_edge_is_last( + [pos(10.0), pos(6.0), pos(4.0), pos(2.0)], + 1.0, + )] + #[case::min_edge_is_third( + [pos(20.0), pos(10.0), pos(6.0), pos(8.0)], + 3.0, + )] + fn maximum_allowable(#[case] edge_lengths: [PositiveReal; 4], #[case] expected: f64) { + let cuboid = Hypercuboid { edge_lengths }; + assert_eq!(cuboid.maximum_allowable_interaction_range(), expected); + } + + #[rstest] + #[case::interior( + [5.0, 3.0, 8.0, -1.0], + [5.0, 3.0, 8.0, -1.0], + )] + #[case::lower_bound( + [-10.0, -10.0, -10.0, -10.0], + [-10.0, -10.0, -10.0, -10.0], + )] + #[case::just_below_upper( + [10.0_f64.next_down(), 10.0_f64.next_down(), 10.0_f64.next_down(), 10.0_f64.next_down()], + [10.0_f64.next_down(), 10.0_f64.next_down(), 10.0_f64.next_down(), 10.0_f64.next_down()], + )] + #[case::at_upper( + [10.0, 10.0, 10.0, 10.0], + [-10.0, -10.0, -10.0, -10.0], + )] + #[case::wrap_around( + [25.0, -35.0, 55.0, -15.0], + [5.0, 5.0, -5.0, 5.0], + )] + fn wrap(#[case] input: [f64; 4], #[case] expected: [f64; 4]) { + let cuboid = Hypercuboid { + edge_lengths: [pos(20.0); 4], + }; + let periodic = Periodic::new(0.0, cuboid).expect("hard-coded range should be valid"); + assert_eq!( + periodic.wrap(Point::new(input.into())), + Ok(Point::new(expected.into())), + ); + } + + #[test] + fn no_ghosts() { + let cuboid = Hypercuboid { + edge_lengths: [pos(40.0), pos(20.0), pos(10.0), pos(8.0)], + }; + let periodic = Periodic::new(1.0, cuboid).expect("hard-coded range should be valid"); + + let inner = Hypercuboid { + edge_lengths: [pos(38.0), pos(18.0), pos(8.0), pos(6.0)], + }; + + let mut rng = StdRng::seed_from_u64(1); + for _ in 0..N_SAMPLES { + let point = inner.sample(&mut rng); + let ghosts = periodic.generate_ghosts(&Point::new(point)); + assert!(ghosts.is_empty()); + } + } + + /// lx=20, ly=10, lz=40, lw=8 with interaction range 1.0. + fn unequal_box() -> Periodic> { + let cuboid = Hypercuboid { + edge_lengths: [pos(20.0), pos(10.0), pos(40.0), pos(8.0)], + }; + Periodic::new(1.0, cuboid).expect("hard-coded range should be valid") + } + + #[rstest] + #[case::outside_x([10.5, 0.0, 0.0, 0.0], 0)] + #[case::outside_w([0.0, 0.0, 0.0, 4.5], 0)] + #[case::boundary_max_x([9.0, 0.0, 0.0, 0.0], 0)] + #[case::boundary_min_x([-9.0, 0.0, 0.0, 0.0], 0)] + #[case::boundary_max_y([0.0, 4.0, 0.0, 0.0], 0)] + #[case::boundary_min_y([0.0, -4.0, 0.0, 0.0], 0)] + #[case::boundary_max_z([0.0, 0.0, 19.0, 0.0], 0)] + #[case::boundary_min_z([0.0, 0.0, -19.0, 0.0], 0)] + #[case::boundary_max_w([0.0, 0.0, 0.0, 3.0], 0)] + #[case::boundary_min_w([0.0, 0.0, 0.0, -3.0], 0)] + #[case::face_pos_x([9.5, 0.0, 0.0, 0.0], 1)] + #[case::face_neg_x([-9.5, 0.0, 0.0, 0.0], 1)] + #[case::face_pos_y([0.0, 4.5, 0.0, 0.0], 1)] + #[case::face_pos_w([0.0, 0.0, 0.0, 3.5], 1)] + #[case::face_neg_w([0.0, 0.0, 0.0, -3.5], 1)] + #[case::edge_pos_x_pos_y([9.5, 4.5, 0.0, 0.0], 3)] + #[case::edge_pos_x_neg_y_pos_z([9.5, -4.5, 19.5, 0.0], 7)] + fn ghosts(#[case] input: [f64; 4], #[case] expected_count: usize) { + let periodic = unequal_box(); + let ghosts = periodic.generate_ghosts(&Point::new(input.into())); + assert_eq!(ghosts.len(), expected_count); + } + + #[test] + fn ghosts_corner_all_pos() { + if MAX_GHOSTS < 15 { + return; + } + let periodic = unequal_box(); + let ghosts = periodic.generate_ghosts(&Point::new([9.5, 4.5, 19.5, 3.5].into())); + assert_eq!(ghosts.len(), 15); + } + + #[rstest] + #[case::face_pos_x( + [9.5, 0.0, 0.0, 0.0], + [[-10.5, 0.0, 0.0, 0.0]].to_vec(), + )] + #[case::face_neg_x( + [-9.5, 0.0, 0.0, 0.0], + [[10.5, 0.0, 0.0, 0.0]].to_vec(), + )] + #[case::face_pos_y( + [0.0, 4.5, 0.0, 0.0], + [[0.0, -5.5, 0.0, 0.0]].to_vec(), + )] + #[case::face_pos_w( + [0.0, 0.0, 0.0, 3.5], + [[0.0, 0.0, 0.0, -4.5]].to_vec(), + )] + #[case::face_neg_w( + [0.0, 0.0, 0.0, -3.5], + [[0.0, 0.0, 0.0, 4.5]].to_vec(), + )] + #[case::edge_pos_x_pos_y( + [9.5, 4.5, 0.0, 0.0], + [[-10.5, -5.5, 0.0, 0.0], [9.5, -5.5, 0.0, 0.0], [-10.5, 4.5, 0.0, 0.0]].to_vec(), + )] + #[case::edge_pos_x_neg_y_pos_z( + [9.5, -4.5, 19.5, 0.0], + [ + [-10.5, 5.5, -20.5, 0.0], + [9.5, 5.5, -20.5, 0.0], + [-10.5, -4.5, -20.5, 0.0], + [9.5, -4.5, -20.5, 0.0], + [-10.5, 5.5, 19.5, 0.0], + [9.5, 5.5, 19.5, 0.0], + [-10.5, -4.5, 19.5, 0.0], + ].to_vec(), + )] + fn ghost_positions(#[case] input: [f64; 4], #[case] expected: Vec<[f64; 4]>) { + let periodic = unequal_box(); + let ghosts = periodic.generate_ghosts(&Point::new(input.into())); + assert_eq!(ghosts.len(), expected.len()); + for (ghost, expected_pos) in ghosts.iter().zip(expected.iter()) { + assert_relative_eq!(ghost.position, (*expected_pos).into()); + } } } } From 5b39a5ee3b02cdc8e464e2f0497fad97f8a9a142 Mon Sep 17 00:00:00 2001 From: janbridley Date: Fri, 26 Jun 2026 16:40:47 -0400 Subject: [PATCH 2/8] bench --- hoomd-microstate/benches/cuboid.rs | 343 +++++++++++++++++++++++++++++ 1 file changed, 343 insertions(+) create mode 100644 hoomd-microstate/benches/cuboid.rs diff --git a/hoomd-microstate/benches/cuboid.rs b/hoomd-microstate/benches/cuboid.rs new file mode 100644 index 000000000..bcb411aaf --- /dev/null +++ b/hoomd-microstate/benches/cuboid.rs @@ -0,0 +1,343 @@ +// Copyright (c) 2024-2026 The Regents of the University of Michigan. +// Part of hoomd-rs, released under the BSD 3-Clause License. + +#![expect( + clippy::missing_docs_in_private_items, + reason = "benches don't need public documentation" +)] +#![expect(clippy::wildcard_imports, reason = "simplifies code")] +#![expect( + clippy::needless_pass_by_value, + reason = "divan takes Bencher by value" +)] + +//! Benchmark periodic ghost generation for [`Hypercuboid`]. +//! +//! The crate ships a single generic `GenerateGhosts` impl for +//! `Periodic>` that enumerates periodic images with a bitmask +//! loop. This benchmark compares it to the hand-unrolled, dimension-specialized +//! 2D/3D implementations it replaced. +//! +//! The buffer capacity is the crate's real `MAX_GHOSTS` (set at build time via +//! `HOOMD_MAX_GHOSTS`), read here from the same environment variable so the +//! self-contained `specific_*` variants allocate the *same* buffer the shipped +//! `generate_ghosts` does. Because `MAX_GHOSTS` is a single global constant, +//! each dimension is measured at its own optimal capacity by running the bench +//! separately with `HOOMD_MAX_GHOSTS = 2^N - 1` and filtering to that +//! dimension's rows (e.g. `HOOMD_MAX_GHOSTS=7 cargo bench --bench cuboid -- 3d`). +//! +//! * `trait_method` — the shipped `generate_ghosts` (the new generic bitmask +//! algorithm). Run for the dimension matching the configured `MAX_GHOSTS`. +//! * `specific` — the original hand-unrolled dimension-specialized algorithm, +//! reimplemented here for direct comparison. +//! +//! Each is measured on a `corner` site (the worst case, `2^N - 1` ghosts) and an +//! `interior` site (the common case, no ghosts). + +use arrayvec::ArrayVec; +use divan::{Bencher, black_box}; +use hoomd_geometry::{IsPointInside, shape::Hypercuboid}; +use hoomd_microstate::{ + boundary::{GenerateGhosts, Periodic}, + property::Point, +}; +use hoomd_utility::valid::PositiveReal; +use hoomd_vector::Cartesian; + +fn main() { + divan::main(); +} + +/// The crate's `MAX_GHOSTS`, read from the same `HOOMD_MAX_GHOSTS` build-time +/// environment variable so the bench and the shipped code agree on the buffer +/// size. +const MAX_GHOSTS: usize = match option_env!("HOOMD_MAX_GHOSTS") { + // `from_str` is not const, so manually parse like the crate does. + Some(val) => match usize::from_str_radix(val, 10) { + Ok(n) => n, + Err(_) => panic!("HOOMD_MAX_GHOSTS must be a non-negative integer"), + }, + None => 12, +}; + +/// Edge length and interaction range used by every benchmark. +const EDGE: f64 = 20.0; +const RANGE: f64 = 1.0; + +/// Build a `PositiveReal` from a hard-coded constant. +fn pos(value: f64) -> PositiveReal { + value + .try_into() + .expect("hard-coded constant should be positive") +} + +/// Build a cubic periodic box. +fn periodic_box(edge: f64, range: f64) -> Periodic> { + let shape = Hypercuboid::with_equal_edges(pos(edge)); + Periodic::new(range, shape).expect("hard-coded range should be valid") +} + +// ------------------------------------------------------------------------------------------------- +// The original hand-unrolled, dimension-specialized implementations. +// +// These are verbatim reimplementations of the pre-generalization `generate_ghosts` +// bodies, so they measure exactly the algorithm the generic impl replaced. +// `min`/`max` are hoisted once, matching the originals. +// ------------------------------------------------------------------------------------------------- + +/// The original hand-unrolled 2D implementation. +fn specific_2d( + shape: &Hypercuboid<2>, + range: f64, + site: &Point>, +) -> ArrayVec>, MAX_GHOSTS> { + let mut result = ArrayVec::new(); + let r = site.position; + if !shape.is_point_inside(&r) { + return result; + } + + let max = shape.maximal_extents(); + let min = shape.minimal_extents(); + + let new_site = |x: f64, y: f64| { + let mut new_site = *site; + new_site.position[0] += x * shape.edge_lengths[0].get(); + new_site.position[1] += y * shape.edge_lengths[1].get(); + new_site + }; + + let near_left = r[0] < min[0] + range; + let near_right = r[0] > max[0] - range; + let near_top = r[1] > max[1] - range; + let near_bottom = r[1] < min[1] + range; + + if near_right { + result.push(new_site(-1.0, 0.0)); + } + if near_left { + result.push(new_site(1.0, 0.0)); + } + if near_top { + result.push(new_site(0.0, -1.0)); + } + if near_bottom { + result.push(new_site(0.0, 1.0)); + } + if near_right && near_top { + result.push(new_site(-1.0, -1.0)); + } + if near_right && near_bottom { + result.push(new_site(-1.0, 1.0)); + } + if near_left && near_top { + result.push(new_site(1.0, -1.0)); + } + if near_left && near_bottom { + result.push(new_site(1.0, 1.0)); + } + + result +} + +/// The original hand-unrolled 3D implementation. +#[allow( + clippy::too_many_lines, + reason = "mirrors the original hand-unrolled code" +)] +fn specific_3d( + shape: &Hypercuboid<3>, + range: f64, + site: &Point>, +) -> ArrayVec>, MAX_GHOSTS> { + let mut result = ArrayVec::new(); + let r = site.position; + if !shape.is_point_inside(&r) { + return result; + } + + let max = shape.maximal_extents(); + let min = shape.minimal_extents(); + + let new_site = |x: f64, y: f64, z: f64| { + let mut new_site = *site; + new_site.position[0] += x * shape.edge_lengths[0].get(); + new_site.position[1] += y * shape.edge_lengths[1].get(); + new_site.position[2] += z * shape.edge_lengths[2].get(); + new_site + }; + + let near_left = r[0] < min[0] + range; + let near_right = r[0] > max[0] - range; + let near_top = r[1] > max[1] - range; + let near_bottom = r[1] < min[1] + range; + let near_front = r[2] > max[2] - range; + let near_back = r[2] < min[2] + range; + + if near_right { + result.push(new_site(-1.0, 0.0, 0.0)); + } + if near_left { + result.push(new_site(1.0, 0.0, 0.0)); + } + if near_top { + result.push(new_site(0.0, -1.0, 0.0)); + } + if near_bottom { + result.push(new_site(0.0, 1.0, 0.0)); + } + if near_front { + result.push(new_site(0.0, 0.0, -1.0)); + } + if near_back { + result.push(new_site(0.0, 0.0, 1.0)); + } + + if near_right && near_top { + result.push(new_site(-1.0, -1.0, 0.0)); + } + if near_right && near_bottom { + result.push(new_site(-1.0, 1.0, 0.0)); + } + if near_right && near_front { + result.push(new_site(-1.0, 0.0, -1.0)); + } + if near_right && near_back { + result.push(new_site(-1.0, 0.0, 1.0)); + } + if near_left && near_top { + result.push(new_site(1.0, -1.0, 0.0)); + } + if near_left && near_bottom { + result.push(new_site(1.0, 1.0, 0.0)); + } + if near_left && near_front { + result.push(new_site(1.0, 0.0, -1.0)); + } + if near_left && near_back { + result.push(new_site(1.0, 0.0, 1.0)); + } + + if near_top && near_front { + result.push(new_site(0.0, -1.0, -1.0)); + } + if near_bottom && near_front { + result.push(new_site(0.0, 1.0, -1.0)); + } + if near_top && near_back { + result.push(new_site(0.0, -1.0, 1.0)); + } + if near_bottom && near_back { + result.push(new_site(0.0, 1.0, 1.0)); + } + + if near_right && near_top && near_front { + result.push(new_site(-1.0, -1.0, -1.0)); + } + if near_right && near_top && near_back { + result.push(new_site(-1.0, -1.0, 1.0)); + } + if near_right && near_bottom && near_front { + result.push(new_site(-1.0, 1.0, -1.0)); + } + if near_right && near_bottom && near_back { + result.push(new_site(-1.0, 1.0, 1.0)); + } + if near_left && near_top && near_front { + result.push(new_site(1.0, -1.0, -1.0)); + } + if near_left && near_top && near_back { + result.push(new_site(1.0, -1.0, 1.0)); + } + if near_left && near_bottom && near_front { + result.push(new_site(1.0, 1.0, -1.0)); + } + if near_left && near_bottom && near_back { + result.push(new_site(1.0, 1.0, 1.0)); + } + + result +} + +// ------------------------------------------------------------------------------------------------- +// Benchmarks. +// ------------------------------------------------------------------------------------------------- + +/// The shipped `generate_ghosts` (the new generic bitmask algorithm). +#[divan::bench_group] +mod trait_method { + use super::*; + + macro_rules! corner { + ($($name:ident at $n:literal),+ $(,)?) => { + $(#[divan::bench] + fn $name(bencher: Bencher) { + let periodic = periodic_box::<$n>(EDGE, RANGE); + let coord = EDGE / 2.0 - RANGE / 2.0; + bencher + .with_inputs(|| Cartesian::<$n>::from([coord; $n])) + .bench_local_values(|pt| { + black_box(periodic.generate_ghosts(&Point::new(pt))) + }); + })+ + }; + } + + macro_rules! interior { + ($($name:ident at $n:literal),+ $(,)?) => { + $(#[divan::bench] + fn $name(bencher: Bencher) { + let periodic = periodic_box::<$n>(EDGE, RANGE); + bencher + .with_inputs(|| Cartesian::<$n>::from([0.0; $n])) + .bench_local_values(|pt| { + black_box(periodic.generate_ghosts(&Point::new(pt))) + }); + })+ + }; + } + + corner! { + corner_1d at 1, corner_2d at 2, corner_3d at 3, corner_4d at 4, corner_5d at 5, + } + + interior! { + interior_1d at 1, interior_2d at 2, interior_3d at 3, interior_4d at 4, interior_5d at 5, + } +} + +/// The original hand-unrolled dimension-specialized implementations. +#[divan::bench_group] +mod specific { + use super::*; + + macro_rules! corner { + ($name:ident, $n:literal, $alg:path) => { + #[divan::bench] + fn $name(bencher: Bencher) { + let shape = Hypercuboid::<$n>::with_equal_edges(pos(EDGE)); + let coord = EDGE / 2.0 - RANGE / 2.0; + bencher + .with_inputs(|| Cartesian::<$n>::from([coord; $n])) + .bench_local_values(|pt| black_box($alg(&shape, RANGE, &Point::new(pt)))); + } + }; + } + + macro_rules! interior { + ($name:ident, $n:literal, $alg:path) => { + #[divan::bench] + fn $name(bencher: Bencher) { + let shape = Hypercuboid::<$n>::with_equal_edges(pos(EDGE)); + bencher + .with_inputs(|| Cartesian::<$n>::from([0.0; $n])) + .bench_local_values(|pt| black_box($alg(&shape, RANGE, &Point::new(pt)))); + } + }; + } + + corner! { corner_2d, 2, specific_2d } + corner! { corner_3d, 3, specific_3d } + interior! { interior_2d, 2, specific_2d } + interior! { interior_3d, 3, specific_3d } +} From c474ce9befd24f9a6373439baa3e47da82489bcd Mon Sep 17 00:00:00 2001 From: janbridley Date: Fri, 26 Jun 2026 16:40:51 -0400 Subject: [PATCH 3/8] bench in cargo.toml --- hoomd-microstate/Cargo.toml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/hoomd-microstate/Cargo.toml b/hoomd-microstate/Cargo.toml index 66e5c990f..ede8332e2 100644 --- a/hoomd-microstate/Cargo.toml +++ b/hoomd-microstate/Cargo.toml @@ -43,3 +43,7 @@ log.workspace = true rand.workspace = true rstest.workspace = true tempfile.workspace = true + +[[bench]] +name = "cuboid" +harness = false From b378f9b51ac5ea23fbeee9bf995f3fcd2858db73 Mon Sep 17 00:00:00 2001 From: janbridley Date: Fri, 26 Jun 2026 16:52:48 -0400 Subject: [PATCH 4/8] lint --- hoomd-microstate/benches/cuboid.rs | 4 ---- 1 file changed, 4 deletions(-) diff --git a/hoomd-microstate/benches/cuboid.rs b/hoomd-microstate/benches/cuboid.rs index bcb411aaf..a11f3a25c 100644 --- a/hoomd-microstate/benches/cuboid.rs +++ b/hoomd-microstate/benches/cuboid.rs @@ -6,10 +6,6 @@ reason = "benches don't need public documentation" )] #![expect(clippy::wildcard_imports, reason = "simplifies code")] -#![expect( - clippy::needless_pass_by_value, - reason = "divan takes Bencher by value" -)] //! Benchmark periodic ghost generation for [`Hypercuboid`]. //! From a652c7fd1664dac8cd19a83c772f7b183e49c1a5 Mon Sep 17 00:00:00 2001 From: janbridley Date: Fri, 26 Jun 2026 16:54:00 -0400 Subject: [PATCH 5/8] test against previous behavior --- .../src/boundary/periodic/cuboid.rs | 261 ++++++++++++++++++ 1 file changed, 261 insertions(+) diff --git a/hoomd-microstate/src/boundary/periodic/cuboid.rs b/hoomd-microstate/src/boundary/periodic/cuboid.rs index 5cc8e0ffb..be8da86af 100644 --- a/hoomd-microstate/src/boundary/periodic/cuboid.rs +++ b/hoomd-microstate/src/boundary/periodic/cuboid.rs @@ -958,4 +958,265 @@ mod tests { } } } + + /// The generic `generate_ghosts` must produce the same *set* of ghosts as the + /// hand-unrolled implementations it replaced, for every input + mod differential { + use super::*; + use rand::RngExt; + + fn pos(value: f64) -> PositiveReal { + value + .try_into() + .expect("hard-coded constant should be positive") + } + + /// The pre-generalization hand-unrolled 2D implementation, kept as a + /// reference for the generic [`GenerateGhosts::generate_ghosts`]. + fn reference_2d( + shape: &Hypercuboid<2>, + range: f64, + site: &Point>, + ) -> ArrayVec>, MAX_GHOSTS> { + let mut result = ArrayVec::new(); + let r = site.position; + if !shape.is_point_inside(&r) { + return result; + } + let max = shape.maximal_extents(); + let min = shape.minimal_extents(); + let new_site = |x: f64, y: f64| { + let mut new_site = *site; + new_site.position[0] += x * shape.edge_lengths[0].get(); + new_site.position[1] += y * shape.edge_lengths[1].get(); + new_site + }; + let near_left = r[0] < min[0] + range; + let near_right = r[0] > max[0] - range; + let near_top = r[1] > max[1] - range; + let near_bottom = r[1] < min[1] + range; + if near_right { + result.push(new_site(-1.0, 0.0)); + } + if near_left { + result.push(new_site(1.0, 0.0)); + } + if near_top { + result.push(new_site(0.0, -1.0)); + } + if near_bottom { + result.push(new_site(0.0, 1.0)); + } + if near_right && near_top { + result.push(new_site(-1.0, -1.0)); + } + if near_right && near_bottom { + result.push(new_site(-1.0, 1.0)); + } + if near_left && near_top { + result.push(new_site(1.0, -1.0)); + } + if near_left && near_bottom { + result.push(new_site(1.0, 1.0)); + } + result + } + + /// The pre-generalization hand-unrolled 3D implementation, kept as a + /// reference for the generic [`GenerateGhosts::generate_ghosts`]. + #[allow( + clippy::too_many_lines, + reason = "mirrors the original hand-unrolled code" + )] + fn reference_3d( + shape: &Hypercuboid<3>, + range: f64, + site: &Point>, + ) -> ArrayVec>, MAX_GHOSTS> { + let mut result = ArrayVec::new(); + let r = site.position; + if !shape.is_point_inside(&r) { + return result; + } + let max = shape.maximal_extents(); + let min = shape.minimal_extents(); + let new_site = |x: f64, y: f64, z: f64| { + let mut new_site = *site; + new_site.position[0] += x * shape.edge_lengths[0].get(); + new_site.position[1] += y * shape.edge_lengths[1].get(); + new_site.position[2] += z * shape.edge_lengths[2].get(); + new_site + }; + let near_left = r[0] < min[0] + range; + let near_right = r[0] > max[0] - range; + let near_top = r[1] > max[1] - range; + let near_bottom = r[1] < min[1] + range; + let near_front = r[2] > max[2] - range; + let near_back = r[2] < min[2] + range; + if near_right { + result.push(new_site(-1.0, 0.0, 0.0)); + } + if near_left { + result.push(new_site(1.0, 0.0, 0.0)); + } + if near_top { + result.push(new_site(0.0, -1.0, 0.0)); + } + if near_bottom { + result.push(new_site(0.0, 1.0, 0.0)); + } + if near_front { + result.push(new_site(0.0, 0.0, -1.0)); + } + if near_back { + result.push(new_site(0.0, 0.0, 1.0)); + } + if near_right && near_top { + result.push(new_site(-1.0, -1.0, 0.0)); + } + if near_right && near_bottom { + result.push(new_site(-1.0, 1.0, 0.0)); + } + if near_right && near_front { + result.push(new_site(-1.0, 0.0, -1.0)); + } + if near_right && near_back { + result.push(new_site(-1.0, 0.0, 1.0)); + } + if near_left && near_top { + result.push(new_site(1.0, -1.0, 0.0)); + } + if near_left && near_bottom { + result.push(new_site(1.0, 1.0, 0.0)); + } + if near_left && near_front { + result.push(new_site(1.0, 0.0, -1.0)); + } + if near_left && near_back { + result.push(new_site(1.0, 0.0, 1.0)); + } + if near_top && near_front { + result.push(new_site(0.0, -1.0, -1.0)); + } + if near_bottom && near_front { + result.push(new_site(0.0, 1.0, -1.0)); + } + if near_top && near_back { + result.push(new_site(0.0, -1.0, 1.0)); + } + if near_bottom && near_back { + result.push(new_site(0.0, 1.0, 1.0)); + } + if near_right && near_top && near_front { + result.push(new_site(-1.0, -1.0, -1.0)); + } + if near_right && near_top && near_back { + result.push(new_site(-1.0, -1.0, 1.0)); + } + if near_right && near_bottom && near_front { + result.push(new_site(-1.0, 1.0, -1.0)); + } + if near_right && near_bottom && near_back { + result.push(new_site(-1.0, 1.0, 1.0)); + } + if near_left && near_top && near_front { + result.push(new_site(1.0, -1.0, -1.0)); + } + if near_left && near_top && near_back { + result.push(new_site(1.0, -1.0, 1.0)); + } + if near_left && near_bottom && near_front { + result.push(new_site(1.0, 1.0, -1.0)); + } + if near_left && near_bottom && near_back { + result.push(new_site(1.0, 1.0, 1.0)); + } + result + } + + /// Sample a point whose coordinates each fall in one of five zones: the + /// lower boundary band, the upper boundary band, the interior, just + /// outside the upper face, or just outside the lower face. This exercises + /// interior, face, edge, corner, and out-of-bounds inputs uniformly. + fn sample_point( + min: [f64; N], + max: [f64; N], + range: f64, + rng: &mut StdRng, + ) -> Cartesian { + let mut coords = [0.0_f64; N]; + for i in 0..N { + let (lo, hi) = (min[i], max[i]); + let t = rng.random::(); + coords[i] = match rng.random::() % 5 { + 0 => lo + t * range, + 1 => hi - t * range, + 2 => lo + range + t * (hi - lo - 2.0 * range).max(0.0), + 3 => hi + t * range, + _ => lo - t * range, + }; + } + Cartesian::from(coords) + } + + /// Assert two ghost collections hold the same positions, ignoring order. + fn assert_same_set( + actual: &ArrayVec>, MAX_GHOSTS>, + expected: &ArrayVec>, MAX_GHOSTS>, + ) { + assert_eq!(actual.len(), expected.len()); + let mut actual: Vec> = actual.iter().map(|ghost| ghost.position).collect(); + let mut expected: Vec> = + expected.iter().map(|ghost| ghost.position).collect(); + actual.sort_by(|a, b| a.coordinates.partial_cmp(&b.coordinates).unwrap()); + expected.sort_by(|a, b| a.coordinates.partial_cmp(&b.coordinates).unwrap()); + for (actual, expected) in actual.iter().zip(expected.iter()) { + assert_relative_eq!(actual, expected); + } + } + + #[test] + fn matches_reference_2d() { + let mut rng = StdRng::seed_from_u64(0xC0DE); + for (edge_lengths, range) in [ + ([pos(20.0), pos(10.0)], 1.0), + ([pos(5.0), pos(5.0)], 2.0), + ([pos(100.0), pos(3.0)], 1.0), + ] { + let cuboid = Hypercuboid { edge_lengths }; + let periodic = + Periodic::new(range, cuboid.clone()).expect("hard-coded range should be valid"); + let min = cuboid.minimal_extents(); + let max = cuboid.maximal_extents(); + for _ in 0..4096 { + let site = Point::new(sample_point(min, max, range, &mut rng)); + let actual = periodic.generate_ghosts(&site); + let expected = reference_2d(&cuboid, range, &site); + assert_same_set(&actual, &expected); + } + } + } + + #[test] + fn matches_reference_3d() { + let mut rng = StdRng::seed_from_u64(0xBEEF); + for (edge_lengths, range) in [ + ([pos(20.0), pos(10.0), pos(40.0)], 1.0), + ([pos(6.0), pos(6.0), pos(6.0)], 2.0), + ([pos(80.0), pos(5.0), pos(30.0)], 1.0), + ] { + let cuboid = Hypercuboid { edge_lengths }; + let periodic = + Periodic::new(range, cuboid.clone()).expect("hard-coded range should be valid"); + let min = cuboid.minimal_extents(); + let max = cuboid.maximal_extents(); + for _ in 0..4096 { + let site = Point::new(sample_point(min, max, range, &mut rng)); + let actual = periodic.generate_ghosts(&site); + let expected = reference_3d(&cuboid, range, &site); + assert_same_set(&actual, &expected); + } + } + } + } } From 479b497c8f104c673a2e8b57dd4ad4c0cf1560b2 Mon Sep 17 00:00:00 2001 From: janbridley Date: Fri, 26 Jun 2026 17:05:03 -0400 Subject: [PATCH 6/8] fmt +nightly --- hoomd-microstate/src/boundary/periodic/cuboid.rs | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/hoomd-microstate/src/boundary/periodic/cuboid.rs b/hoomd-microstate/src/boundary/periodic/cuboid.rs index be8da86af..60c385b63 100644 --- a/hoomd-microstate/src/boundary/periodic/cuboid.rs +++ b/hoomd-microstate/src/boundary/periodic/cuboid.rs @@ -119,7 +119,7 @@ where /// ``` /// use hoomd_geometry::shape::Rectangle; /// use hoomd_microstate::{ - /// boundary::{Periodic, GenerateGhosts}, + /// boundary::{GenerateGhosts, Periodic}, /// property::Point, /// }; /// use hoomd_vector::Cartesian; @@ -128,7 +128,8 @@ where /// let periodic = /// Periodic::new(1.0, Rectangle::with_equal_edges(10.0.try_into()?))?; /// // A site near the right edge produces an image shifted across it. - /// let ghosts = periodic.generate_ghosts(&Point::new(Cartesian::from([4.6, 0.0]))); + /// let ghosts = + /// periodic.generate_ghosts(&Point::new(Cartesian::from([4.6, 0.0]))); /// assert_eq!(ghosts.len(), 1); /// # Ok(()) /// # } @@ -1177,7 +1178,7 @@ mod tests { #[test] fn matches_reference_2d() { - let mut rng = StdRng::seed_from_u64(0xC0DE); + let mut rng = StdRng::seed_from_u64(0xc0de); for (edge_lengths, range) in [ ([pos(20.0), pos(10.0)], 1.0), ([pos(5.0), pos(5.0)], 2.0), @@ -1199,7 +1200,7 @@ mod tests { #[test] fn matches_reference_3d() { - let mut rng = StdRng::seed_from_u64(0xBEEF); + let mut rng = StdRng::seed_from_u64(0xbeef); for (edge_lengths, range) in [ ([pos(20.0), pos(10.0), pos(40.0)], 1.0), ([pos(6.0), pos(6.0), pos(6.0)], 2.0), From fd4ebb6566d2d48105457dee8723d93fc46d036d Mon Sep 17 00:00:00 2001 From: janbridley Date: Fri, 26 Jun 2026 17:09:06 -0400 Subject: [PATCH 7/8] fix comment --- hoomd-microstate/src/boundary/periodic/cuboid.rs | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/hoomd-microstate/src/boundary/periodic/cuboid.rs b/hoomd-microstate/src/boundary/periodic/cuboid.rs index 60c385b63..d3b4cf38a 100644 --- a/hoomd-microstate/src/boundary/periodic/cuboid.rs +++ b/hoomd-microstate/src/boundary/periodic/cuboid.rs @@ -110,9 +110,8 @@ where /// site near the middle of an `(N - 1)`-face produces one; a /// site near an edge (where two boundaries meet) produces three; and so on, up to /// `2.pow(N) - 1` ghosts for a site near all `2 * N` boundaries at once. - /// The images are emitted in descending order of the set of active directions. - /// Images are currently emitted in descending order of the subset bitmask (with - /// higher bit positions corresponding to folds across a higher-dimensional facet) + /// Images are currently emitted in descending numeric order of the subset bitmask, + /// where bit i indicates the image is folded across the facet normal to dimension i /// /// # Example /// From 3bf1026feae9584d241c7ca9fa66a995cb4c4849 Mon Sep 17 00:00:00 2001 From: janbridley Date: Fri, 26 Jun 2026 18:12:02 -0400 Subject: [PATCH 8/8] remove bench, does not work on trunk --- hoomd-microstate/benches/cuboid.rs | 339 ----------------------------- 1 file changed, 339 deletions(-) delete mode 100644 hoomd-microstate/benches/cuboid.rs diff --git a/hoomd-microstate/benches/cuboid.rs b/hoomd-microstate/benches/cuboid.rs deleted file mode 100644 index a11f3a25c..000000000 --- a/hoomd-microstate/benches/cuboid.rs +++ /dev/null @@ -1,339 +0,0 @@ -// Copyright (c) 2024-2026 The Regents of the University of Michigan. -// Part of hoomd-rs, released under the BSD 3-Clause License. - -#![expect( - clippy::missing_docs_in_private_items, - reason = "benches don't need public documentation" -)] -#![expect(clippy::wildcard_imports, reason = "simplifies code")] - -//! Benchmark periodic ghost generation for [`Hypercuboid`]. -//! -//! The crate ships a single generic `GenerateGhosts` impl for -//! `Periodic>` that enumerates periodic images with a bitmask -//! loop. This benchmark compares it to the hand-unrolled, dimension-specialized -//! 2D/3D implementations it replaced. -//! -//! The buffer capacity is the crate's real `MAX_GHOSTS` (set at build time via -//! `HOOMD_MAX_GHOSTS`), read here from the same environment variable so the -//! self-contained `specific_*` variants allocate the *same* buffer the shipped -//! `generate_ghosts` does. Because `MAX_GHOSTS` is a single global constant, -//! each dimension is measured at its own optimal capacity by running the bench -//! separately with `HOOMD_MAX_GHOSTS = 2^N - 1` and filtering to that -//! dimension's rows (e.g. `HOOMD_MAX_GHOSTS=7 cargo bench --bench cuboid -- 3d`). -//! -//! * `trait_method` — the shipped `generate_ghosts` (the new generic bitmask -//! algorithm). Run for the dimension matching the configured `MAX_GHOSTS`. -//! * `specific` — the original hand-unrolled dimension-specialized algorithm, -//! reimplemented here for direct comparison. -//! -//! Each is measured on a `corner` site (the worst case, `2^N - 1` ghosts) and an -//! `interior` site (the common case, no ghosts). - -use arrayvec::ArrayVec; -use divan::{Bencher, black_box}; -use hoomd_geometry::{IsPointInside, shape::Hypercuboid}; -use hoomd_microstate::{ - boundary::{GenerateGhosts, Periodic}, - property::Point, -}; -use hoomd_utility::valid::PositiveReal; -use hoomd_vector::Cartesian; - -fn main() { - divan::main(); -} - -/// The crate's `MAX_GHOSTS`, read from the same `HOOMD_MAX_GHOSTS` build-time -/// environment variable so the bench and the shipped code agree on the buffer -/// size. -const MAX_GHOSTS: usize = match option_env!("HOOMD_MAX_GHOSTS") { - // `from_str` is not const, so manually parse like the crate does. - Some(val) => match usize::from_str_radix(val, 10) { - Ok(n) => n, - Err(_) => panic!("HOOMD_MAX_GHOSTS must be a non-negative integer"), - }, - None => 12, -}; - -/// Edge length and interaction range used by every benchmark. -const EDGE: f64 = 20.0; -const RANGE: f64 = 1.0; - -/// Build a `PositiveReal` from a hard-coded constant. -fn pos(value: f64) -> PositiveReal { - value - .try_into() - .expect("hard-coded constant should be positive") -} - -/// Build a cubic periodic box. -fn periodic_box(edge: f64, range: f64) -> Periodic> { - let shape = Hypercuboid::with_equal_edges(pos(edge)); - Periodic::new(range, shape).expect("hard-coded range should be valid") -} - -// ------------------------------------------------------------------------------------------------- -// The original hand-unrolled, dimension-specialized implementations. -// -// These are verbatim reimplementations of the pre-generalization `generate_ghosts` -// bodies, so they measure exactly the algorithm the generic impl replaced. -// `min`/`max` are hoisted once, matching the originals. -// ------------------------------------------------------------------------------------------------- - -/// The original hand-unrolled 2D implementation. -fn specific_2d( - shape: &Hypercuboid<2>, - range: f64, - site: &Point>, -) -> ArrayVec>, MAX_GHOSTS> { - let mut result = ArrayVec::new(); - let r = site.position; - if !shape.is_point_inside(&r) { - return result; - } - - let max = shape.maximal_extents(); - let min = shape.minimal_extents(); - - let new_site = |x: f64, y: f64| { - let mut new_site = *site; - new_site.position[0] += x * shape.edge_lengths[0].get(); - new_site.position[1] += y * shape.edge_lengths[1].get(); - new_site - }; - - let near_left = r[0] < min[0] + range; - let near_right = r[0] > max[0] - range; - let near_top = r[1] > max[1] - range; - let near_bottom = r[1] < min[1] + range; - - if near_right { - result.push(new_site(-1.0, 0.0)); - } - if near_left { - result.push(new_site(1.0, 0.0)); - } - if near_top { - result.push(new_site(0.0, -1.0)); - } - if near_bottom { - result.push(new_site(0.0, 1.0)); - } - if near_right && near_top { - result.push(new_site(-1.0, -1.0)); - } - if near_right && near_bottom { - result.push(new_site(-1.0, 1.0)); - } - if near_left && near_top { - result.push(new_site(1.0, -1.0)); - } - if near_left && near_bottom { - result.push(new_site(1.0, 1.0)); - } - - result -} - -/// The original hand-unrolled 3D implementation. -#[allow( - clippy::too_many_lines, - reason = "mirrors the original hand-unrolled code" -)] -fn specific_3d( - shape: &Hypercuboid<3>, - range: f64, - site: &Point>, -) -> ArrayVec>, MAX_GHOSTS> { - let mut result = ArrayVec::new(); - let r = site.position; - if !shape.is_point_inside(&r) { - return result; - } - - let max = shape.maximal_extents(); - let min = shape.minimal_extents(); - - let new_site = |x: f64, y: f64, z: f64| { - let mut new_site = *site; - new_site.position[0] += x * shape.edge_lengths[0].get(); - new_site.position[1] += y * shape.edge_lengths[1].get(); - new_site.position[2] += z * shape.edge_lengths[2].get(); - new_site - }; - - let near_left = r[0] < min[0] + range; - let near_right = r[0] > max[0] - range; - let near_top = r[1] > max[1] - range; - let near_bottom = r[1] < min[1] + range; - let near_front = r[2] > max[2] - range; - let near_back = r[2] < min[2] + range; - - if near_right { - result.push(new_site(-1.0, 0.0, 0.0)); - } - if near_left { - result.push(new_site(1.0, 0.0, 0.0)); - } - if near_top { - result.push(new_site(0.0, -1.0, 0.0)); - } - if near_bottom { - result.push(new_site(0.0, 1.0, 0.0)); - } - if near_front { - result.push(new_site(0.0, 0.0, -1.0)); - } - if near_back { - result.push(new_site(0.0, 0.0, 1.0)); - } - - if near_right && near_top { - result.push(new_site(-1.0, -1.0, 0.0)); - } - if near_right && near_bottom { - result.push(new_site(-1.0, 1.0, 0.0)); - } - if near_right && near_front { - result.push(new_site(-1.0, 0.0, -1.0)); - } - if near_right && near_back { - result.push(new_site(-1.0, 0.0, 1.0)); - } - if near_left && near_top { - result.push(new_site(1.0, -1.0, 0.0)); - } - if near_left && near_bottom { - result.push(new_site(1.0, 1.0, 0.0)); - } - if near_left && near_front { - result.push(new_site(1.0, 0.0, -1.0)); - } - if near_left && near_back { - result.push(new_site(1.0, 0.0, 1.0)); - } - - if near_top && near_front { - result.push(new_site(0.0, -1.0, -1.0)); - } - if near_bottom && near_front { - result.push(new_site(0.0, 1.0, -1.0)); - } - if near_top && near_back { - result.push(new_site(0.0, -1.0, 1.0)); - } - if near_bottom && near_back { - result.push(new_site(0.0, 1.0, 1.0)); - } - - if near_right && near_top && near_front { - result.push(new_site(-1.0, -1.0, -1.0)); - } - if near_right && near_top && near_back { - result.push(new_site(-1.0, -1.0, 1.0)); - } - if near_right && near_bottom && near_front { - result.push(new_site(-1.0, 1.0, -1.0)); - } - if near_right && near_bottom && near_back { - result.push(new_site(-1.0, 1.0, 1.0)); - } - if near_left && near_top && near_front { - result.push(new_site(1.0, -1.0, -1.0)); - } - if near_left && near_top && near_back { - result.push(new_site(1.0, -1.0, 1.0)); - } - if near_left && near_bottom && near_front { - result.push(new_site(1.0, 1.0, -1.0)); - } - if near_left && near_bottom && near_back { - result.push(new_site(1.0, 1.0, 1.0)); - } - - result -} - -// ------------------------------------------------------------------------------------------------- -// Benchmarks. -// ------------------------------------------------------------------------------------------------- - -/// The shipped `generate_ghosts` (the new generic bitmask algorithm). -#[divan::bench_group] -mod trait_method { - use super::*; - - macro_rules! corner { - ($($name:ident at $n:literal),+ $(,)?) => { - $(#[divan::bench] - fn $name(bencher: Bencher) { - let periodic = periodic_box::<$n>(EDGE, RANGE); - let coord = EDGE / 2.0 - RANGE / 2.0; - bencher - .with_inputs(|| Cartesian::<$n>::from([coord; $n])) - .bench_local_values(|pt| { - black_box(periodic.generate_ghosts(&Point::new(pt))) - }); - })+ - }; - } - - macro_rules! interior { - ($($name:ident at $n:literal),+ $(,)?) => { - $(#[divan::bench] - fn $name(bencher: Bencher) { - let periodic = periodic_box::<$n>(EDGE, RANGE); - bencher - .with_inputs(|| Cartesian::<$n>::from([0.0; $n])) - .bench_local_values(|pt| { - black_box(periodic.generate_ghosts(&Point::new(pt))) - }); - })+ - }; - } - - corner! { - corner_1d at 1, corner_2d at 2, corner_3d at 3, corner_4d at 4, corner_5d at 5, - } - - interior! { - interior_1d at 1, interior_2d at 2, interior_3d at 3, interior_4d at 4, interior_5d at 5, - } -} - -/// The original hand-unrolled dimension-specialized implementations. -#[divan::bench_group] -mod specific { - use super::*; - - macro_rules! corner { - ($name:ident, $n:literal, $alg:path) => { - #[divan::bench] - fn $name(bencher: Bencher) { - let shape = Hypercuboid::<$n>::with_equal_edges(pos(EDGE)); - let coord = EDGE / 2.0 - RANGE / 2.0; - bencher - .with_inputs(|| Cartesian::<$n>::from([coord; $n])) - .bench_local_values(|pt| black_box($alg(&shape, RANGE, &Point::new(pt)))); - } - }; - } - - macro_rules! interior { - ($name:ident, $n:literal, $alg:path) => { - #[divan::bench] - fn $name(bencher: Bencher) { - let shape = Hypercuboid::<$n>::with_equal_edges(pos(EDGE)); - bencher - .with_inputs(|| Cartesian::<$n>::from([0.0; $n])) - .bench_local_values(|pt| black_box($alg(&shape, RANGE, &Point::new(pt)))); - } - }; - } - - corner! { corner_2d, 2, specific_2d } - corner! { corner_3d, 3, specific_3d } - interior! { interior_2d, 2, specific_2d } - interior! { interior_3d, 3, specific_3d } -}