rcb: avoid re-computations of min/max/sum

instead, keep track of the bounding boxes and the sum of the parts.

This reduces the load of high-iter-count runs.

Author: hhirtz

src/algorithms/recursive_bisection.rs

@@ -1,6 +1,7 @@
 use super::Error;
 use crate::geometry::OrientedBoundingBox;
 use crate::geometry::PointND;
+use crate::BoundingBox;
 use nalgebra::allocator::Allocator;
 use nalgebra::ArrayStorage;
 use nalgebra::Const;
@@ -25,17 +26,30 @@ struct Item<'p, const D: usize, W> {
     part: &'p AtomicUsize,
 }
 
+/// Return value of [rcb_split].
+struct SplitResult<W> {
+    /// Index of the first item in the right part in the array of [Item]s.
+    split_idx: usize,
+    /// Weight of the left part, used to compute the sum for the next iteration.
+    weight_left: W,
+    /// Coordinate value of the split, used to compute the [BoundingBox]es.
+    split_pos: f64,
+}
+
 fn rcb_split<const D: usize, W>(
     items: &mut [Item<D, W>],
     coord: usize,
     tolerance: f64,
     mut min: f64,
     mut max: f64,
     sum: W,
-) -> usize
+) -> SplitResult<W>
 where
     W: RcbWeight,
 {
+    let span = tracing::info_span!("rcb_split");
+    let _enter = span.enter();
+
     let mut prev_count_left = usize::MAX;
     loop {
         let split_target = (min + max) / 2.0;
@@ -57,7 +71,11 @@ where
             f64::abs((weight_left - ideal_weight_left) / ideal_weight_left)
         };
         if count_left == prev_count_left || imbalance < tolerance {
-            return count_left;
+            return SplitResult {
+                split_idx: count_left,
+                weight_left,
+                split_pos: split_target,
+            };
         }
         prev_count_left = count_left;
 
@@ -76,6 +94,8 @@ fn rcb_recurse<const D: usize, W>(
     iter_id: usize,
     coord: usize,
     tolerance: f64,
+    sum: W,
+    bb: BoundingBox<D>,
 ) where
     W: RcbWeight,
 {
@@ -84,33 +104,35 @@ fn rcb_recurse<const D: usize, W>(
         return;
     }
     if iter_count == 0 {
+        let span = tracing::info_span!(
+            "rcb_recurse: apply_part_id",
+            item_count = items.len(),
+            iter_id,
+        );
+        let _enter = span.enter();
+
         items
             .into_par_iter()
             .for_each(|item| item.part.store(iter_id, Ordering::Relaxed));
         return;
     }
 
-    let sum: W = items.par_iter().map(|item| item.weight).sum();
-    let (min, max) = items
-        .par_iter()
-        .fold(
-            || (f64::INFINITY, f64::NEG_INFINITY),
-            |(min, max), item| {
-                (
-                    f64::min(min, item.point[coord]),
-                    f64::max(max, item.point[coord]),
-                )
-            },
-        )
-        .reduce(
-            || (f64::INFINITY, f64::NEG_INFINITY),
-            |(min0, max0), (min1, max1)| (f64::min(min0, min1), f64::max(max0, max1)),
-        );
+    let span = tracing::info_span!("rcb_recurse");
+    let enter = span.enter();
 
-    let split_idx = rcb_split(items, coord, tolerance, min, max, sum);
+    let min = bb.p_min[coord];
+    let max = bb.p_max[coord];
+    let SplitResult {
+        split_idx,
+        weight_left,
+        split_pos,
+    } = rcb_split(items, coord, tolerance, min, max, sum);
     let (left, right) = if split_idx == items.len() {
         items.split_at_mut(items.len())
     } else {
+        let span = tracing::info_span!("select_nth_unstable");
+        let _enter = span.enter();
+
         let (left, _, _right_minus_one) = items
             .select_nth_unstable_by(split_idx, |item1, item2| {
                 f64::partial_cmp(&item1.point[coord], &item2.point[coord]).unwrap()
@@ -119,6 +141,13 @@ fn rcb_recurse<const D: usize, W>(
         items.split_at_mut(left_len)
     };
 
+    let mut bb_left = bb.clone();
+    bb_left.p_max[coord] = split_pos;
+    let mut bb_right = bb;
+    bb_right.p_min[coord] = split_pos;
+
+    mem::drop(enter);
+
     rayon::join(
         || {
             rcb_recurse(
@@ -127,6 +156,8 @@ fn rcb_recurse<const D: usize, W>(
                 2 * iter_id + 1,
                 (coord + 1) % D,
                 tolerance,
+                weight_left,
+                bb_left,
             )
         },
         || {
@@ -136,6 +167,8 @@ fn rcb_recurse<const D: usize, W>(
                 2 * iter_id + 2,
                 (coord + 1) % D,
                 tolerance,
+                sum - weight_left,
+                bb_right,
             )
         },
     );
@@ -171,7 +204,7 @@ where
         });
     }
     if partition.is_empty() {
-        // Would make the partition.min() at the end panic.
+        // Would make BoundingBox::from_points and partition.min() panic.
         return Ok(());
     }
 
@@ -184,9 +217,12 @@ where
             part,
         })
         .collect();
+    let sum = items.par_iter().map(|item| item.weight).sum();
+    let bb = BoundingBox::from_points(items.par_iter().map(|item| item.point)).unwrap();
 
-    rcb_recurse(&mut items, iter_count, 0, 0, tolerance);
+    rcb_recurse(&mut items, iter_count, 0, 0, tolerance, sum, bb);
 
+    // Part IDs must start from zero.
     let part_id_offset = *partition.par_iter().min().unwrap();
     partition
         .par_iter_mut()