Space-group search: take the operator disagreement as a median, not a mean

A merohedral twin mixes EVERY reflection with its twin mate, so it shifts the
whole distribution of |I1-I2|/(I1+I2). A minority of badly measured
reflections shifts only the tail. The mean cannot tell those apart; the median
is blind to the second and just as sensitive to the first.

Measured on real crystals, moving the statistic from the mean to the median
leaves genuine promotions where they are and pushes every twin up:

  genuine tetragonal    1.016 -> 1.013
  genuine lysozyme      1.051 -> 1.067
  genuine tetragonal    1.238 -> 1.231
  twin (-e 1050)        1.272 -> 1.447
  twin (-e 450)         1.280 -> 1.622
  twin (full)           1.441 -> 1.522
  twin (-e 600)         1.427 -> 2.010

The margin around the 1.25 bound widens from 2.7% (genuine 1.238 against twin
1.272 - uncomfortably tight for a decision that cannot be undone downstream)
to 17.5% (1.231 against 1.447). The bound itself does not move.

Rotation battery, 33 crystals in both detection modes: no point group changed
in either (30/33 and 29/33, as before), and only one crystal's numbers move at
all - the one already documented as nondeterministic between repeat runs of
the same binary. The synthetic twin-fraction x multiplicity grid passes
unchanged. So this buys margin, not outcomes.

Found while testing a different hypothesis, which the same measurement refuted:
a tetragonal crystal whose 422 promotion is wrongly refused reads 1.484 by the
mean and 1.472 by the median, i.e. its disagreement is distribution-wide and is
NOT a badly-integrated minority. That crystal's cause is elsewhere and is not
addressed here - see the note below.

  Its indexing-ambiguity operator (-k,-h,-l) lies INSIDE 422 but OUTSIDE 222,
  so the subgroup merge the search is given mixes lattices indexed in the two
  alternative hands. That corrupts exactly the 4-fold relationships and leaves
  the 2-fold ones intact - measured, the 222 step reads 0.917 and the 422 step
  1.484 - and the corruption is indistinguishable from a twin law. Forcing the
  tetragonal group merges the two hands as equivalent and the same data give
  CC1/2 99.2% at multiplicity 10.7, matching XDS. The failure is worse the
  BETTER the frames index (99.9% vs 63.3% for the run that gets it right),
  because indexing more frames picks up more of both hands.

  So no statistic computed on a subgroup merge can arbitrate a promotion whose
  added operators include an indexing-ambiguity operator. Fixing that means
  resolving the ambiguity before the search, or detecting the coincidence and
  deciding another way; the operators needed to detect it are already computed
  (the run warns about them).

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

image_analysis/scale_merge/SearchSpaceGroup.cpp

@@ -273,16 +273,18 @@ SearchSpaceGroupResult SearchSpaceGroup(
         s.n_pairs = static_cast<int>(x.size());
         s.cc = PearsonCC(x, y);
         // Sigma-free disagreement over the same pairs (see SpaceGroupOptions::max_operator_h_ratio).
-        double h_sum = 0.0;
-        int h_n = 0;
+        std::vector<double> hv;
+        hv.reserve(x.size());
         for (size_t p = 0; p < x.size(); ++p) {
             const double denom = x[p] + y[p];
-            if (denom > 0.0) {
-                h_sum += std::fabs(x[p] - y[p]) / denom;
-                ++h_n;
-            }
+            if (denom > 0.0)
+                hv.push_back(std::fabs(x[p] - y[p]) / denom);
+        }
+        if (!hv.empty()) {
+            const size_t mid = hv.size() / 2;
+            std::nth_element(hv.begin(), hv.begin() + mid, hv.end());
+            s.h_stat = hv[mid];
         }
-        s.h_stat = h_n > 0 ? h_sum / h_n : 0.0;
         s.present = s.n_pairs >= opt.min_pairs_per_operator && std::isfinite(s.cc) &&
                     s.cc >= opt.min_operator_cc;
         return s;
@@ -536,7 +538,7 @@ SearchSpaceGroupResult SearchSpaceGroup(
                 const bool in_parent = std::binary_search(parent->rotation_set.begin(),
                                                           parent->rotation_set.end(), RotKey(rot));
                 if (in_parent) { h_par += os.h_stat; ++n_par; pairs_par += os.n_pairs; }
-                else            { h_new += os.h_stat; ++n_new; pairs_new += os.n_pairs; }
+                else           { h_new += os.h_stat; ++n_new; pairs_new += os.n_pairs; }
             }
             if (n_new > 0 && n_par > 0 && pairs_new >= opt.min_pairs_for_h
                 && pairs_par >= opt.min_pairs_for_h && h_par > 0.0) {

image_analysis/scale_merge/SearchSpaceGroup.h

@@ -30,11 +30,12 @@ struct SpaceGroupOperatorScore {
     double cc = 0.0;              // correlation of I(h) with I(Rh)
     int n_pairs = 0;             // independent reflection pairs the CC was computed from
     bool present = false;        // operator confirmed as a real symmetry of the intensities
-    // Mean |I1-I2|/(I1+I2) over this operator's pairs - the disagreement the operator implies, with no
+    // MEDIAN |I1-I2|/(I1+I2) over this operator's pairs - the disagreement the operator implies, with no
     // sigma in it. Unlike chi^2 and the systematic-b, which are ratios to a merge error model that
     // drifts with multiplicity, this is a property of the intensities alone; compared against an
     // operator already confirmed on the same reflections it is what separates a real symmetry from a
-    // twin law (see max_operator_h_ratio).
+    // twin law (see max_operator_h_ratio). The median rather than the mean because a merohedral twin
+    // perturbs EVERY pair while a badly integrated minority perturbs only the tail.
     double h_stat = 0.0;
 };
 
@@ -131,16 +132,20 @@ struct SearchSpaceGroupOptions {
     // is a non-negligible fraction of I; below the floor the increase is treated as noise, not a twin.
     double min_systematic_b_for_veto = 0.05;
 
-    // Promotion gate on the operator disagreement H = <|I1-I2|/(I1+I2)>, taken as the ratio of the
+    // Promotion gate on the operator disagreement H = median|I1-I2|/(I1+I2), taken as the ratio of the
     // operators the promotion ADDS to the operators of the parent group already confirmed on the same
     // reflections. A real symmetry operator relates equal intensities, so its H matches the parent's
     // (ratio ~1); a merohedral twin law relates DIFFERENT reflections mixed in proportion alpha, so its
-    // H is systematically larger. Measured over 27 runs spanning 5 promotion types and 450-1800 images:
-    // genuine symmetry 0.862-1.219, merohedral twins 1.270-2.084 - a clean gap, and unlike the chi^2 and
-    // systematic-b ratios (genuine 1.00-3.47 / 1.09-3.89 vs twin 1.35-3.32 / 1.77-4.76, fully
-    // interleaved) it does not drift with multiplicity, because there is no sigma in it and the parent
-    // normalisation cancels data quality. For a partial twin <H> over the twin law is (1-2*alpha)/2, so
-    // the excess over the parent also estimates the twin fraction rather than being a tuned constant.
+    // H is systematically larger. Unlike the chi^2 and systematic-b ratios (genuine 1.00-3.47 /
+    // 1.09-3.89 vs twin 1.35-3.32 / 1.77-4.76, fully interleaved) it does not drift with multiplicity,
+    // because there is no sigma in it and the parent normalisation cancels data quality.
+    //
+    // A MEDIAN, not a mean. A merohedral twin mixes every reflection with its twin mate, so it shifts the
+    // whole distribution of |I1-I2|/(I1+I2); a minority of badly measured reflections shifts only the
+    // tail. Measured on the same real crystals, moving from the mean to the median leaves genuine
+    // promotions where they are (1.016 -> 1.013, 1.051 -> 1.067, 1.238 -> 1.231) and pushes every real
+    // twin UP (1.280 -> 1.622, 1.427 -> 2.010, 1.272 -> 1.447, 1.441 -> 1.522), widening the margin
+    // around this bound from 2.7% to 17.5%. Battery-neutral: 33 crystals, no point group changed.
     //
     // The parent normalisation is what makes this work, and it is not optional. Symmetry-related
     // reflections never agree exactly on real data - absorption, illumination and partiality differ

tests/SearchSpaceGroupTwinTest.cpp

@@ -149,9 +149,9 @@ namespace {
         return out;
     }
 
-    // The statistic the promotion is actually decided on: the mean operator disagreement
-    // H = <|I1-I2|/(I1+I2)> over the operators the promotion ADDS, divided by the mean over the parent
-    // group's own operators, measured on the same reflections. Mirrors what SearchSpaceGroup computes
+    // The statistic the promotion is actually decided on: the operator disagreement
+    // H = median|I1-I2|/(I1+I2) over the operators the promotion ADDS, divided by the same over the
+    // parent group's own operators, measured on the same reflections. Mirrors what SearchSpaceGroup computes
     // for the sub -> super step, so a test can report the margin the max_operator_h_ratio bound has.
     //
     // The parent normalisation is the whole design, not a detail. An ABSOLUTE per-operator bound cannot
@@ -172,6 +172,7 @@ namespace {
         for (const auto& s : operators) {
             if (s.n_pairs < 200)                    // SearchSpaceGroupOptions::min_pairs_for_h
                 continue;
+            // The median, which is what the promotion is gated on (SearchSpaceGroup).
             if (parent_ops.count(s.op_triplet_hkl) > 0) { h_parent += s.h_stat; ++n_parent; }
             else                                        { h_added  += s.h_stat; ++n_added; }
         }