From ca0a9d4dfe3aec58106ebb21b0bcbbca398f814f Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 20 Sep 2026 17:02:24 +0200 Subject: [PATCH] Keep an absence at its row's own strength in the screw evidence it refutes A screw zone's evidence drops its largest member to defend against one badly-measured reflection whose sigma lies about it, and such a reflection is by construction weak - a fraction of the row it sits on. A member that has passed the violation test AND stands at or above the mean of its row's own present class is a different animal: not a measurement that moved, but a reflection that is there. Trimming it removed the single datum that refutes the claim, and the violation-count deferral then read the inflated evidence to forgive the very violation that had been trimmed out of it. Nested screw ORDERS are decided entirely on this. 6_1 extinguishes l != 6n and 6_2/6_4 extinguish l != 3n, so the two differ only on l = 3n not 6n. On a hexagonal crystal whose 00l row holds five present reflections, the strongest of the whole row lay in that difference: trimmed, 6_1 read the row as perfectly dead and won on the count of absences alone - nine at 50.6 nats with one violation against seven at 43.3 with none - and the run reported the wrong screw order with the right one ranked below it. With the violation left in, 6_1 reads 6.8 and is refused. An independent POINTLESS run on the same P1 merge puts the 6_1 condition at probability 0.000 and the 3n condition at 0.998. Trim only among members that are not both flagged present and at full row strength. Both halves of the condition are needed and the corpus separates them: the reflection above stands at 1.94 of its row's mean, where two monoclinic crystals whose 0k0 are genuinely dead carry one violation each at 0.31 and 0.74 of their row - the mis-measurement the trim exists for, and one that costs a real 2_1 if it stays in. Zones with no violations are bit-identical, and so is every candidate whose absent class is clean. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_013nW6FNRP1bBJJ8pfHiByAT --- .../scale_merge/SearchSpaceGroup.cpp | 29 ++++++- tests/SearchSpaceGroupTest.cpp | 79 +++++++++++++++++++ 2 files changed, 105 insertions(+), 3 deletions(-) diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index 0a74c7e67..f180daa6c 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -2066,9 +2066,10 @@ SearchSpaceGroupResult SearchSpaceGroup( // see the deferral below - so there the count is the one instrument still standing. auto& z = zones[a.row]; z.n_absent += 1; - if (a.i_over_sigma > present_cut && + const bool violation = a.i_over_sigma > present_cut && (opt.present_e_squared <= 0.0 || - a.e_squared_expected > opt.present_e_squared * row_scale)) { + a.e_squared_expected > opt.present_e_squared * row_scale); + if (violation) { ++screw_violations; // Per zone as well: the count-veto deferral below is licensed zone by zone. z.n_violations += 1; @@ -2077,7 +2078,29 @@ SearchSpaceGroupResult SearchSpaceGroup( if (mit != row_mean.end() && mit->second > 0.0) { 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); + // AN ABSENCE AT ITS ROW'S OWN MEAN IS NOT AN OUTLIER. TrimmedZoneSum drops the + // zone's largest member to defend against one badly-measured reflection whose sigma + // lies about it, and such a reflection is by construction WEAK - a fraction of the + // row it sits on. A member that has passed the violation test AND stands at or above + // the mean of its row's own present class is a different animal: not a measurement + // that moved, but a reflection that is there. Trimming it removes the single datum + // that refutes the claim, and the inflated evidence is then read by the deferral + // below to forgive the very violation that was trimmed out of it. + // + // Nested screw ORDERS are decided entirely on this. 6_1 predicts l = 6n and 6_2/6_4 + // predict l = 3n, so the two differ only on l = 3n not 6n; where the strongest + // reflection of the whole row lies in that difference, trimming it let 6_1 read the + // row as perfectly dead and win on the count of absences alone - measured on a + // hexagonal crystal, nine absences at 50.6 nats with one violation beat seven at + // 43.3 with none, and read 6.8 with the violation left in. + // + // Both halves of the condition are needed, and the corpus separates them cleanly: + // that reflection stands at 1.94 of its row's mean, where a monoclinic crystal whose + // eight 0k0 are genuinely dead carries one violation at 0.31 of its row - the + // mis-measurement the trim exists for, and one that must stay trimmable or a real + // 2_1 is lost. Everything below the row mean is trimmed exactly as before. + if (!(violation && u >= 1.0)) + z.max_u = std::max(z.max_u, u); z.n_control = static_cast(row_present_esq.at(a.row).size()); } } diff --git a/tests/SearchSpaceGroupTest.cpp b/tests/SearchSpaceGroupTest.cpp index cd193ae87..88123bdd3 100644 --- a/tests/SearchSpaceGroupTest.cpp +++ b/tests/SearchSpaceGroupTest.cpp @@ -578,6 +578,85 @@ TEST_CASE("A screw zone's evidence does not hang on its largest absence", "[Sear CHECK(TrimmedZoneSum(0.0, 0.0, 8) == 0.0); } +// Screw ORDERS on one axial row are nested: 6_1 extinguishes l != 6n and 6_2/6_4 extinguish +// l != 3n, so 6_1's absent class is 6_2's plus the l = 3n that are not 6n, and the whole of the +// evidence between the two lies in that difference. Trimming the zone's largest member defends +// against one badly-measured reflection, but where the difference class holds the strongest +// reflection on the row it trimmed away the only datum that refutes 6_1 - which then read the row +// as dead, won on its two extra absences, and used the inflated evidence to excuse the very +// violation it had discarded. A member that is both flagged PRESENT and standing at its row's own +// mean is not an outlier, so it is not trimmed. Measured on a hexagonal crystal: nine absences at +// 50.6 nats with one violation beat seven at 43.3 with none, and read 6.8 once the violation - at +// 1.94 of its row - stayed in. The two sections here are the two sides of that bound. +TEST_CASE("SearchSpaceGroup does not trim away the reflection that refutes a screw order", + "[SearchSpaceGroup]") { + SearchSpaceGroupOptions opt; + opt.merge_friedel = true; + + SECTION("an absence at its row's own strength decides against the order that claims it") { + const gemmi::SpaceGroup& sg = gemmi::get_spacegroup_by_name("P 64"); + auto merged = GenerateMergedReflectionsForSpaceGroup(sg, 12); + + // The 00l row as such a crystal records it: l = 3n present, everything else dead, and one + // l = 3n that is NOT 6n - the class 6_1 has to call absent and 6_4 does not - by far the + // strongest reflection on the row. + int on_row = 0; + for (auto& r : merged) { + if (r.h != 0 || r.k != 0) + continue; + const int l = std::abs(r.l); + ++on_row; + if (l % 3 != 0) r.I = 0.0f; // extinguished by the 3n condition, in both candidates + else if (l % 6 == 0) r.I = 300.0f; // the control class 6_1 keeps for itself + else r.I = (l == 9) ? 4000.0f : 5.0f; + } + REQUIRE(on_row >= 8); + + const auto result = SearchSpaceGroup(merged, opt); + INFO(SearchSpaceGroupResultToText(result)); + REQUIRE(result.best_space_group.has_value()); + std::vector accepted{result.best_space_group->short_name()}; + for (const auto& alt : result.alternatives) + accepted.push_back(alt.short_name()); + // P6_2 and P6_4 are enantiomorphs and indistinguishable from intensities; P6_1/P6_5 are a + // different claim and must not be what comes out. + CHECK(std::find(accepted.begin(), accepted.end(), "P64") != accepted.end()); + CHECK(std::find(accepted.begin(), accepted.end(), "P61") == accepted.end()); + CHECK(std::find(accepted.begin(), accepted.end(), "P65") == accepted.end()); + } + + SECTION("one weak absence that moved is still an outlier, and the screw survives it") { + // The other side of the bound, and the case the trim was built for: a genuine 2_1 whose + // 0k0-odd class is dead but for one reflection at a third of its row, measured with a sigma + // that makes it read present. Trimmable as before - it is nowhere near the row's strength - + // and losing that would cost a real screw. + const gemmi::SpaceGroup& sg = gemmi::get_spacegroup_by_name("P 1 21 1"); + auto merged = GenerateMergedReflectionsForSpaceGroup(sg, 18); + + const gemmi::GroupOps gops = sg.operations(); + int absent = 0; + float row_strength = 0.0f; + for (const auto& r : merged) + if (r.h == 0 && r.l == 0 && !gops.is_systematically_absent(gemmi::Op::Miller{{r.h, r.k, r.l}})) + row_strength = std::max(row_strength, r.I); + for (auto& r : merged) { + if (r.h != 0 || r.l != 0) + continue; + if (!gops.is_systematically_absent(gemmi::Op::Miller{{r.h, r.k, r.l}})) + continue; + ++absent; + r.I = absent == 1 ? 0.3f * row_strength : 0.0f; + r.sigma = 0.3f; + } + REQUIRE(absent >= 6); + + const auto result = SearchSpaceGroup(merged, opt); + INFO(SearchSpaceGroupResultToText(result)); + REQUIRE(result.best_space_group.has_value()); + CHECK(result.best_space_group->short_name() == "P21"); + } +} + // 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