Space group: a screw zone's evidence no longer hangs on its largest absence

A screw's predicted-absent class is one axial row - half a dozen to a few dozen
reflections - and its evidence is a SUM over them, so it is decided by its
largest member. The file's own LIMIT comment said so; 7n2s is that limit firing
on real data. Between two scaling passes that differed only in which weak frames
were rejected, one of eight dead 0k0 moved from 14 +- 9 to 99 +- 10 while the
other seven did not move at all, and the zone fell from 30.1 nats to 17.1 and
lost the 2(1) under a bound of 20. That reflection was never measured to the
precision its sigma claimed: its two half-set merges read 198 and 2.5.

The zone's sum is now taken with its single largest member dropped and rescaled
for the trim - divided by n - H_n, the expected sum of the other n-1 under the
null, and multiplied back by n. ScrewZoneEvidence reads the result exactly as
before: same statistic, same floor, same bound, same calibration, with a robust
estimate of the zone's deadness in place of a fragile one. One member only,
whatever the zone holds: a zone with two strong absences is a zone that is not
extinct. On a uniformly dead zone the rescale under-states by 2.2 nats at eight
absences and 3.9 at sixty-four - it only ever refuses, never claims. Glide zones
keep the untrimmed sum: a plane holds hundreds to thousands of reflections and
no single one can carry the verdict.

7n2s -> P 1 21 1 (zone 27.6 nats, set by the seven reflections that did not
move), matching its deposit; a second monoclinic crystal decided six nats under
the bound (7 absent, 1 violation, 13.9 nats) reaches 21.9 and its 2(1) as well.
Unchanged on 7mzt, 7k1l, 11if, 9hs7, 9zlo and four in-house reference sets.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_013nW6FNRP1bBJJ8pfHiByAT

image_analysis/scale_merge/SearchSpaceGroup.cpp

@@ -960,6 +960,38 @@ double ScrewZoneEvidence(double sum_u, int n_absent) {
     return -a * std::log(u) + std::lgamma(a + 1);
 }
 
+// The LIMIT above, answered where it arises: the zone's sum with its single largest member dropped,
+// rescaled so the result still estimates the same quantity. Under the null the u are Exp(1) and
+// E[max] = H_n = 1 + 1/2 + ... + 1/n, so the remaining n-1 of them sum to n - H_n in expectation and
+// dividing by that recovers the mean the whole sum would have given. Multiplying back by n leaves a
+// number ScrewZoneEvidence reads exactly as before - same statistic, same floor, same calibration,
+// with a robust estimate of the zone's deadness in place of a fragile one.
+//
+// A screw row is half a dozen to a few dozen reflections, so ONE of them at a third of its row moves
+// the zone by more than the other seven put together. Measured on a monoclinic crystal whose eight
+// absent 0k0 are dead in every run: between two scaling passes that differed only in which weak
+// frames were rejected, one of those eight moved from 14 +- 9 to 99 +- 10 while the other seven did
+// not move at all, and the zone went from 30.1 nats to 17.1 and lost the 2(1) under a bound of 20.
+// The two half-set merges of that one reflection read 198 and 2.5 in the second pass, i.e. it was
+// never measured to the precision its sigma claimed - but no sigma is read here, on purpose (see
+// min_screw_absence_evidence), so the defence has to be the estimator. Trimmed, that zone reads 27.6
+// nats on the pass that had lost it, set by the seven reflections that did not move.
+//
+// ONE member, whatever the zone's size: a zone with two strong absences is not a zone with a bad
+// reflection in it, it is a zone that is not extinct, and trimming further would hide that. What
+// this costs a genuinely uniform zone is that the estimate is biased high by (n-1)/(n - H_n), i.e.
+// the evidence is UNDER-stated by n*log of that - 2.2 nats on a zone of eight, 3.9 on one of
+// sixty-four. It only ever refuses, never claims.
+double TrimmedZoneSum(double sum_u, double max_u, int n_absent) {
+    if (n_absent < 2)
+        return sum_u;
+    double harmonic = 0.0;
+    for (int i = 1; i <= n_absent; ++i)
+        harmonic += 1.0 / i;
+    const double expected_remaining = n_absent - harmonic;   // > 0 for every n >= 2
+    return n_absent * std::max(0.0, sum_u - max_u) / expected_remaining;
+}
+
 SearchSpaceGroupResult SearchSpaceGroup(
     const std::vector<MergedReflection>& merged,
     const SearchSpaceGroupOptions& opt) {
@@ -1937,7 +1969,8 @@ SearchSpaceGroupResult SearchSpaceGroup(
         // row with no control class of its own cannot say whether its absences are weak or its whole
         // row is, so it contributes no evidence either way rather than being judged against the shell
         // mean; it is carried as an undetermined zone so the report can say so.
-        struct ZoneSums { double sum_u = 0; int n_absent = 0; int n_control = 0; int n_violations = 0; };
+        struct ZoneSums { double sum_u = 0; double max_u = 0; int n_absent = 0; int n_control = 0;
+                          int n_violations = 0; };
         std::map<AxialRow, ZoneSums> zones;
         for (const auto& a : screw_absent_refl) {
             const auto it = row_median.find(a.row);
@@ -1964,7 +1997,9 @@ SearchSpaceGroupResult SearchSpaceGroup(
             }
             const auto mit = row_mean.find(a.row);
             if (mit != row_mean.end() && mit->second > 0.0) {
-                z.sum_u += std::max(0.0, a.e_squared_expected) / mit->second;
+                const double u = std::max(0.0, a.e_squared_expected) / mit->second;
+                z.sum_u += u;
+                z.max_u = std::max(z.max_u, u);
                 z.n_control = static_cast<int>(row_present_esq.at(a.row).size());
             }
         }
@@ -2008,10 +2043,11 @@ SearchSpaceGroupResult SearchSpaceGroup(
         constexpr double SCREW_DEFERRAL_MIN_EVIDENCE_PER_ABSENT = 1.0;
         bool screw_deferral_ok = true;
         for (const auto& [row, z] : zones) {
-            const double ev = z.n_control > 0 ? ScrewZoneEvidence(z.sum_u, z.n_absent) : 0.0;
+            const double sum_u = TrimmedZoneSum(z.sum_u, z.max_u, z.n_absent);
+            const double ev = z.n_control > 0 ? ScrewZoneEvidence(sum_u, z.n_absent) : 0.0;
             s.screw_zones.push_back({row, z.n_absent, z.n_control, ev});
             if (z.n_control > 0) {
-                screw_sum_u += z.sum_u;
+                screw_sum_u += sum_u;
                 screw_scored += z.n_absent;
                 best_zone_evidence = std::max(best_zone_evidence, ev);
             }

image_analysis/scale_merge/SearchSpaceGroup.h

@@ -155,6 +155,13 @@ double AbsenceEvidence(double sum_u, int n_absent, int n_control);
 // candidate's absences could score lower. See SearchSpaceGroup.cpp.
 double ScrewZoneEvidence(double sum_u, int n_absent);
 
+// A screw zone's sum with its single largest member dropped and rescaled to estimate the same
+// quantity (dividing by n - H_n, the expected sum of the other n-1 under the null, and multiplying
+// by n). ScrewZoneEvidence is a sum, so it is dominated by the largest one or two members of a class
+// that holds only a handful of reflections, and a zone's verdict could hang on one reflection whose
+// merged value moved with the scaling. See SearchSpaceGroup.cpp.
+double TrimmedZoneSum(double sum_u, double max_u, int n_absent);
+
 struct SearchSpaceGroupOptions {
     // Lattice (metric) symmetry from LatticeSearch. When set, the point-group search is limited to
     // the Sohncke subgroups of this system's holohedry - the metric is an upper bound on the
@@ -527,6 +534,10 @@ struct SearchSpaceGroupOptions {
     // with I, reading much the same on an absent reflection as on a present one (the <I/s> columns of
     // the candidate table show this directly).
     //
+    // The zone's sum is taken with its largest member trimmed (TrimmedZoneSum), because a sum over
+    // half a dozen reflections is otherwise decided by its largest one and moved tens of nats by any
+    // upstream change that moves that one reflection.
+    //
     // Measured over the probe crystals, an adopted screw reads 24-600 nats and the false ones (the
     // 4_1/4_3 conditions of a crystal that has no screw at all, whose predicted-absent class is as
     // strong as its control row) read -3.8 to -38.5. The bound sits in that gap - loose enough that

tests/SearchSpaceGroupTest.cpp

@@ -499,6 +499,37 @@ TEST_CASE("SearchSpaceGroup does not promote triclinic data on a pseudo-orthorho
     CHECK(result.point_group_order == 1);
 }
 
+// A screw zone is a handful of axial reflections and its evidence is a SUM over them, so it is
+// decided by its largest member. Measured on a monoclinic crystal whose eight absent 0k0 are dead in
+// every run: between two scaling passes that differed only in which weak frames were rejected, ONE
+// of the eight moved from 14 +- 9 to 99 +- 10 (its two half-set merges reading 198 and 2.5, so it
+// was never measured to the precision its sigma claimed) while the other seven did not move at all -
+// and the zone went from 30.1 nats to 17.1 and lost the 2(1) under a bound of 20. Trimming the
+// largest member and rescaling for the trim makes the two passes agree.
+TEST_CASE("A screw zone's evidence does not hang on its largest absence", "[SearchSpaceGroup]") {
+    // Seven reflections at a hundredth of their row, and one that moved between the two passes.
+    const double dead_seven = 7 * 0.01;
+    const double before = ScrewZoneEvidence(TrimmedZoneSum(dead_seven + 0.03, 0.03, 8), 8);
+    const double after = ScrewZoneEvidence(TrimmedZoneSum(dead_seven + 0.39, 0.39, 8), 8);
+
+    CHECK(before == Catch::Approx(after).margin(0.01));   // the same seven reflections, the same verdict
+    CHECK(after > 20.0);                                  // and the screw survives the move
+    // Untrimmed, that one reflection is the whole difference and it crosses the bound.
+    CHECK(ScrewZoneEvidence(dead_seven + 0.03, 8) > 20.0);
+    CHECK(ScrewZoneEvidence(dead_seven + 0.39, 8) < 20.0);
+
+    // Only ONE member is trimmed, whatever the zone holds: a zone with two strong absences is not a
+    // zone with a bad reflection in it, it is a zone that is not extinct.
+    CHECK(ScrewZoneEvidence(TrimmedZoneSum(dead_seven + 0.39 + 0.39, 0.39, 8), 8) < 20.0);
+
+    // On a uniform zone the rescale under-states rather than over-states - the safe direction.
+    CHECK(ScrewZoneEvidence(TrimmedZoneSum(8 * 0.01, 0.01, 8), 8) <
+          ScrewZoneEvidence(8 * 0.01, 8));
+    // A class that merged non-positive throughout is unchanged: sum and max are both zero, so the
+    // floor in ScrewZoneEvidence is what answers, exactly as before.
+    CHECK(TrimmedZoneSum(0.0, 0.0, 8) == 0.0);
+}
+
 // A zone whose predicted absences were never measurable must not outscore a zone that is genuinely
 // dead. sum_u is a sum of max(0, E^2)/row_mean, so it is EXACTLY zero when every absent reflection in
 // the zone merged non-positive - and the Beta tail then diverges, worth ~690 nats per reflection. That