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) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nW6FNRP1bBJJ8pfHiByAT
This commit is contained in:
@@ -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<std::string> 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
|
||||
|
||||
Reference in New Issue
Block a user