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