From bac2d01c44a36dba028bf15966aa098a054f635a Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 20 Sep 2026 08:17:48 +0200 Subject: [PATCH] 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) Claude-Session: https://claude.ai/code/session_013nW6FNRP1bBJJ8pfHiByAT --- .../scale_merge/SearchSpaceGroup.cpp | 44 +++++++++++++++++-- image_analysis/scale_merge/SearchSpaceGroup.h | 11 +++++ tests/SearchSpaceGroupTest.cpp | 31 +++++++++++++ 3 files changed, 82 insertions(+), 4 deletions(-) diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index 3de94bdc8..9f33bf808 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/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& 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 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(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); } diff --git a/image_analysis/scale_merge/SearchSpaceGroup.h b/image_analysis/scale_merge/SearchSpaceGroup.h index c9d291b04..5dcef8b99 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.h +++ b/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 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 diff --git a/tests/SearchSpaceGroupTest.cpp b/tests/SearchSpaceGroupTest.cpp index 2e6361d8e..0b1257926 100644 --- a/tests/SearchSpaceGroupTest.cpp +++ b/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