From 4745a6bb1656e197252eecc081bee1e2013bf96b Mon Sep 17 00:00:00 2001
From: Stephen Waits <steve@waits.net>
Date: Tue, 5 May 2026 12:03:52 -0600
Subject: [PATCH] =?UTF-8?q?perf(hypervolume):=20cut=20HSO=20recursion=20ov?=
 =?UTF-8?q?erhead=20by=20~30=C3=97=20on=20n=3D100/3-D?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The HSO recursion in `hypervolume_nd` had three overheads that
dominated SMS-EMOA's per-generation cost on DTLZ2 (5.6 s baseline,
~30 k generations × ~40 HV calls per generation = ~1.2 M HV calls
per run):

1. `active = sorted.clone()` plus `active.iter().position(...)`
   linear scan to remove the just-processed point each band — O(N)
   per band, total O(N²) per HV call.
2. Per-band re-projection
   `active.iter().map(|q| q[..last].to_vec())` — full
   Vec<Vec<f64>> rebuild for every band, O(N·M) allocations per HV
   call.
3. `non_dominated_projection` called even when recursing into the
   M=2 base case, whose sweep already filters dominated points
   internally.

Replace (1) with prefix-slicing `projected_all[..=k]` (sort points
ascending by last axis once; the active set at each band is just a
prefix). Pre-project once outside the loop (2). Skip the explicit
non-dominance filter when the inner recursion is M=2 (3).

Bit-identical output verified by re-running the compare harness and
diffing against the v0.3.0 snapshot — every quality metric matches
to the last decimal.

gungraun (instructions):
- hypervolume_nd_3d n=30:    676 902 →    87 969  (-87 %, 7.7×)
- hypervolume_nd_3d n=100: 13 523 760 →   422 767 (-97 %, 32×)

Wall-clock (compare harness, 10-seed mean):
- SMS-EMOA / DTLZ2: 5643 ms → 1413 ms (-4230 ms, -75 %)
---
 src/metrics/hypervolume.rs | 46 ++++++++++++++++++++++----------------
 1 file changed, 27 insertions(+), 19 deletions(-)

diff --git a/src/metrics/hypervolume.rs b/src/metrics/hypervolume.rs
index 130502c..875c44c 100644
--- a/src/metrics/hypervolume.rs
+++ b/src/metrics/hypervolume.rs
@@ -219,10 +219,10 @@ fn hso_recursive(points: &[Vec<f64>], reference: &[f64]) -> f64 {
     //     corresponding sub-reference), multiplied by band thickness, is
     //     the band's HV contribution.
     //
-    // Because boxes extend from `p[last]` UP TO `reference[last]`, every
-    // point is active in the band immediately below the reference. We
-    // therefore start with `active = all points` and REMOVE the largest
-    // remaining last-axis point each iteration.
+    // We sort points ascending by the last axis once, then iterate from
+    // the largest last-axis value downward. The active set at iteration
+    // `k` is exactly the prefix `sorted[..=k]` — no allocations or
+    // linear-scan removals needed.
     let last = m - 1;
     let mut sorted: Vec<Vec<f64>> = points.to_vec();
     sorted.sort_by(|a, b| {
@@ -231,24 +231,32 @@ fn hso_recursive(points: &[Vec<f64>], reference: &[f64]) -> f64 {
             .unwrap_or(std::cmp::Ordering::Equal)
     });
 
-    let sub_reference: Vec<f64> = reference[..last].to_vec();
+    // Pre-project all points onto the first M-1 axes once. The active
+    // set at each band is the prefix `projected_all[..=k]`; we slice
+    // that prefix instead of rebuilding it per band.
+    let projected_all: Vec<Vec<f64>> = sorted.iter().map(|q| q[..last].to_vec()).collect();
+    let sub_reference: &[f64] = &reference[..last];
     let mut total = 0.0;
-    let mut active: Vec<Vec<f64>> = sorted.clone();
     let mut prev = reference[last];
-    for p in sorted.into_iter().rev() {
+    for k in (0..sorted.len()).rev() {
+        let p = &sorted[k];
         let depth = prev - p[last];
-        if depth > 0.0 && !active.is_empty() {
-            let projected: Vec<Vec<f64>> = active.iter().map(|q| q[..last].to_vec()).collect();
-            let nd = non_dominated_projection(&projected);
-            total += depth * hso_recursive(&nd, &sub_reference);
-        }
-        // Remove the just-processed point (the one with the largest
-        // remaining last-axis value).
-        let idx = active
-            .iter()
-            .position(|q| (q[last] - p[last]).abs() < 1e-15 && q[..last] == p[..last]);
-        if let Some(i) = idx {
-            active.swap_remove(i);
+        if depth > 0.0 {
+            let active = &projected_all[..=k];
+            // The 2-D base case sweeps in sorted-x order and skips any
+            // point with `y >= last_y`, which is exactly the dominance
+            // filter — so for M=3 (sub_reference len 2) we can hand
+            // `active` straight to `hso_recursive` without paying for
+            // an O(K²) `non_dominated_projection` first. For M≥4 we
+            // still need the explicit filter to keep the recursion's
+            // upper levels honest.
+            let inner = if sub_reference.len() == 2 {
+                hso_recursive(active, sub_reference)
+            } else {
+                let nd = non_dominated_projection(active);
+                hso_recursive(&nd, sub_reference)
+            };
+            total += depth * inner;
         }
         prev = p[last];
     }