From cf4fbcf96a9467c6dc1512318272846b3cbb84db Mon Sep 17 00:00:00 2001
From: Filip Leonarski <filip.leonarski@psi.ch>
Date: Mon, 24 Aug 2026 13:11:24 +0200
Subject: [PATCH] Reduce the anomalous split once per ASU group, not once per
 observation

ComputeAsuGroups states the rule for itself - "one ASU reduction per distinct
raw hkl (not per observation)" - and the anomalous split then did a gemmi ASU
reduction and an unordered_map lookup for every one of the millions of fulls.

Both things it wants are properties of the observation's ASU GROUP rather than
of the observation: group_h/k/l is the group's SIGNED representative, so the
same reduction applied to it returns the Friedel-merged key and the hand
together. Reduce once per group into a dense accumulator indexed from there.
The hand only follows the group when the merge distinguishes the hands; a
Friedel-merged run holds both in one group and still has to ask per observation.

SigAno and the merged statistics are unchanged (2.96 over 53303 acentric pairs,
merge table byte-identical).

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 .../scale_merge/RotationScaleMerge.cpp        | 37 +++++++++++++++----
 1 file changed, 30 insertions(+), 7 deletions(-)

diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp
index a14de6e9..53c7f5a8 100644
--- a/image_analysis/scale_merge/RotationScaleMerge.cpp
+++ b/image_analysis/scale_merge/RotationScaleMerge.cpp
@@ -3146,22 +3146,45 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
         const HKLKeyGenerator anom_keygen(/*merge_friedel=*/false, sg_num);
         const gemmi::GroupOps gops = gemmi::find_spacegroup_by_number(sg_num)->operations();
         struct AnomAcc { double swI[2] = {}; double sw[2] = {}; int32_t h = 0, k = 0, l = 0; float d = NAN; };  // [hand] 0=I(+) 1=I(-)
-        std::unordered_map<uint64_t, AnomAcc> anom;
-        anom.reserve(result.merged.size() * 2 + 1);
+        // The Friedel-merged key an observation lands on, and which hand it is, are properties of its
+        // ASU GROUP, not of the observation - group_h/k/l is that group's SIGNED representative, so the
+        // same reduction applied to it returns both. Reduce once per group (a few hundred thousand)
+        // instead of once per observation (millions), and index a dense accumulator from there, which
+        // is the rule ComputeAsuGroups states for itself. The hand only follows the group when the
+        // merge distinguishes the hands; a Friedel-merged run holds both in one group and must ask.
+        std::vector<int32_t> pair_of_group(n_groups, -1);
+        std::vector<uint8_t> hand_of_group(merge_friedel ? 0 : n_groups, 0);
+        std::vector<HKLKey> pair_hkl;
+        {
+            std::unordered_map<uint64_t, int32_t> pair_id;
+            pair_id.reserve(static_cast<size_t>(n_groups) + 1);
+            for (int g = 0; g < n_groups; ++g) {
+                const HKLKey gk = anom_keygen(group_h[g], group_k[g], group_l[g]);
+                if (!merge_friedel) hand_of_group[g] = gk.plus ? 0 : 1;
+                const uint64_t key = HKLKey{gk.h, gk.k, gk.l, true}.pack();
+                const auto [it, fresh] = pair_id.emplace(key, static_cast<int32_t>(pair_hkl.size()));
+                if (fresh) pair_hkl.push_back(HKLKey{gk.h, gk.k, gk.l, true});
+                pair_of_group[g] = it->second;
+            }
+        }
+        std::vector<AnomAcc> anom(pair_hkl.size());
         for (const auto &o : fulls) {
             if (!usable_merge(o)) continue;
             if (rejected_obs[&o - fulls.data()]) continue;   // outlier-rejected in the merge
-            const HKLKey ak = anom_keygen(o.h, o.k, o.l);
-            const int hand = ak.plus ? 0 : 1;
+            const int hand = merge_friedel ? (anom_keygen(o.h, o.k, o.l).plus ? 0 : 1)
+                                           : hand_of_group[o.group];
             const float I_corr = o.I * o.corr;
             const float sigma_corr = corrected_sigma(o, I_corr, o.sigma * o.corr);
             if (!(sigma_corr > 0.0f) || !std::isfinite(sigma_corr)) continue;
             const double w = 1.0 / (static_cast<double>(sigma_corr) * sigma_corr);
-            AnomAcc &a = anom[HKLKey{ak.h, ak.k, ak.l, true}.pack()];
-            a.h = ak.h; a.k = ak.k; a.l = ak.l; a.d = o.d;
+            AnomAcc &a = anom[pair_of_group[o.group]];
+            a.d = o.d;
             a.swI[hand] += w * static_cast<double>(I_corr); a.sw[hand] += w;
         }
-        for (const auto &[fkey, a] : anom) {
+        for (size_t p = 0; p < anom.size(); ++p) {
+            AnomAcc a = anom[p];
+            a.h = pair_hkl[p].h; a.k = pair_hkl[p].k; a.l = pair_hkl[p].l;
+            const uint64_t fkey = pair_hkl[p].pack();
             // Centrics have I(+)=I(-) by symmetry; leave them without an anomalous split.
             if (gops.is_reflection_centric(gemmi::Op::Miller{a.h, a.k, a.l})) continue;
             AnomExport ex;