From 37821bdd3dc70ccb33c45e25e3e98e6019ec2e6d Mon Sep 17 00:00:00 2001 From: Stephen Waits Date: Thu, 14 May 2026 06:34:47 -0600 Subject: [PATCH] perf(metrics): stop re-sorting hypervolume_nd prefixes per slice (361K -> 291K instr) The HSO M=3 path called the generic 2-D base case for every last-axis slice, which re-sorted the active prefix by axis 0 each time -- O(n^2 log n) overall. Since `projected` is already in last-axis order, sorting the projected indices by axis 0 once and sweeping them with a `pi > k` skip gives O(n^2) with no per-slice allocation. The M>=4 path is unchanged (lifted out of the inner branch verbatim). hypervolume_nd_bench_3d n=100: 361_595 -> 291_247 (-19%, 1.24x); n=30 -16%. The sweep visits points in the same (axis-0, then last-axis) order the stable per-prefix sort produced -- output is bit-identical, all 606 tests pass. Co-Authored-By: Claude Opus 4.7 (1M context) --- src/metrics/hypervolume.rs | 69 +++++++++++++++++++++++++++----------- 1 file changed, 50 insertions(+), 19 deletions(-) diff --git a/src/metrics/hypervolume.rs b/src/metrics/hypervolume.rs index d09071f..be2cc9d 100644 --- a/src/metrics/hypervolume.rs +++ b/src/metrics/hypervolume.rs @@ -279,27 +279,58 @@ fn hso_recursive(points: &[Vec], reference: &[f64]) -> f64 { let sub_reference: &[f64] = &reference[..last]; let mut total = 0.0; let mut prev = reference[last]; - for k in (0..order.len()).rev() { - let p_last = points[order[k]][last]; - let depth = prev - p_last; - if depth > 0.0 { - let active = &projected[..=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 { + + if sub_reference.len() == 2 { + // M == 3: the inner HV is a 2-D staircase sweep. `projected` is in + // last-axis order, so the active set at step `k` is the prefix + // `projected[..=k]`. The generic recursion re-sorts that prefix by + // axis 0 on every step — O(n² log n). Instead, sort the projected + // indices by axis 0 once and, for each `k`, sweep them skipping any + // whose last-axis rank exceeds `k`. The sweep visits points in the + // same (axis-0, then last-axis) order the stable per-prefix sort + // produced, so the result is bit-identical. + let r0 = sub_reference[0]; + let r1 = sub_reference[1]; + let mut x_order: Vec = (0..projected.len()).collect(); + x_order.sort_by(|&a, &b| { + projected[a][0] + .partial_cmp(&projected[b][0]) + .unwrap_or(std::cmp::Ordering::Equal) + }); + for k in (0..order.len()).rev() { + let p_last = points[order[k]][last]; + let depth = prev - p_last; + if depth > 0.0 { + let mut area = 0.0; + let mut last_y = r1; + for &pi in &x_order { + if pi > k { + continue; + } + let p = &projected[pi]; + if p[1] >= last_y { + continue; + } + area += (r0 - p[0]) * (last_y - p[1]); + last_y = p[1]; + } + total += depth * area; + } + prev = p_last; + } + } else { + // M >= 4: recurse generically, with the explicit non-dominated + // filter to keep the recursion's upper levels honest. + for k in (0..order.len()).rev() { + let p_last = points[order[k]][last]; + let depth = prev - p_last; + if depth > 0.0 { + let active = &projected[..=k]; let nd = non_dominated_projection(active); - hso_recursive(&nd, sub_reference) - }; - total += depth * inner; + total += depth * hso_recursive(&nd, sub_reference); + } + prev = p_last; } - prev = p_last; } total