From 6df657fe62960e1092ffde188e1ef3f4010723be Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 3 Jul 2026 15:18:27 +0200 Subject: [PATCH] SG search: gate the systematic-absence test on intensity, not just I/sigma The monoclinic 2_1 screw was missed on EP_cs_02-10 and EP_cs_01-17 (adopted P2 instead of P2_1). The 0k0-odd reflections are correctly measured as weak (~0.1-1% of 0k0-even, matching XDS), but their merged sigmas are ~2x too small, so their I/sigma clears the present_i_over_sigma cut and they count as screw-axis violations - the search then falls back to the symmorphic group. Add a resolution-normalised intensity E^2 = I / (shell), computed from equal-count resolution shells, and require a reflection to reach present_e_squared (0.3) as well as present_i_over_sigma before it counts as violating a predicted absence. This only tightens "present", so it cannot manufacture a screw whose predicted-absent class carries real intensity (a symmorphic crystal's axial reflections sit at E^2 ~ 1 and still register as violations). On the 18-crystal rotation battery this recovers the screw on EP_cs_02-10 and EP_cs_01-17 (-> P2_1) and, as a side effect, on MyoB (-> P2_1), pding4_001 (-> P4_122/P4_322) and pding4_003 (-> P2_2_2_1) - all confirmed by genuine absences in the reference intensities (absent class at 0.02-0.76% of the allowed class), which the old sigma-only test also missed. The other 13 crystals, indexing rate and ISa are unchanged; the screw-free control (Ins_I -> I23/I2_13) is unaffected. Add a regression test that reproduces the under-estimated-sigma screw and checks the gate recovers it (and that disabling the gate reproduces the miss). Co-Authored-By: Claude Opus 4.8 --- .../scale_merge/SearchSpaceGroup.cpp | 32 +++++++++++++- image_analysis/scale_merge/SearchSpaceGroup.h | 7 +++ tests/SearchSpaceGroupTest.cpp | 44 +++++++++++++++++++ 3 files changed, 82 insertions(+), 1 deletion(-) diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index 4835739d..ec434941 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -197,6 +197,31 @@ SearchSpaceGroupResult SearchSpaceGroup( (opt.min_i_over_sigma <= 0 || IoverSigma[i] >= opt.min_i_over_sigma); } + // Resolution-normalised intensity E^2 = I / (shell), from equal-count resolution shells over + // the reflections the absence test uses. Lets the absence test judge "present" by intensity + // magnitude, not by a possibly under-estimated sigma (see present_e_squared). + std::vector Esq(n, 0.0); + { + std::vector order; + order.reserve(n); + for (size_t i = 0; i < n; ++i) + if (pass_absence[i]) + order.push_back(i); + std::sort(order.begin(), order.end(), + [&](size_t a, size_t b) { return merged[a].d > merged[b].d; }); // low res -> high res + const int bins = std::clamp(static_cast(order.size() / 100), 1, 25); + const size_t per = (order.size() + bins - 1) / std::max(1, bins); + for (size_t b = 0; b * per < order.size(); ++b) { + const size_t lo = b * per, hi = std::min(order.size(), lo + per); + double sum = 0.0; + for (size_t j = lo; j < hi; ++j) + sum += I[order[j]]; + const double mean = (hi > lo) ? sum / static_cast(hi - lo) : 0.0; + for (size_t j = lo; j < hi; ++j) + Esq[order[j]] = mean > 0.0 ? I[order[j]] / mean : 0.0; + } + } + std::unordered_map key_to_index; key_to_index.reserve(n * 2); for (size_t i = 0; i < n; ++i) @@ -370,7 +395,12 @@ SearchSpaceGroupResult SearchSpaceGroup( if (!pass_absence[i]) continue; const gemmi::Op::Miller hkl{{H[i], K[i], L[i]}}; - const bool present = IoverSigma[i] > opt.present_i_over_sigma; + // Present := statistically significant AND intensity-significant. The E^2 gate keeps a + // weak axial reflection with an under-estimated sigma (fake high I/sigma) from faking a + // screw-axis violation; it only relaxes "present", so it cannot over-call a screw whose + // predicted-absent class carries real intensity. + const bool present = IoverSigma[i] > opt.present_i_over_sigma && + (opt.present_e_squared <= 0.0 || Esq[i] > opt.present_e_squared); if (CenteringAbsent(gops, hkl)) { s.absent_observed += 1; diff --git a/image_analysis/scale_merge/SearchSpaceGroup.h b/image_analysis/scale_merge/SearchSpaceGroup.h index 82109b9c..c512525f 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.h +++ b/image_analysis/scale_merge/SearchSpaceGroup.h @@ -89,6 +89,13 @@ struct SearchSpaceGroupOptions { // wrongly-assumed absence, and to keep the correlation stage from pairing near-zero reflections. double present_i_over_sigma = 3.0; + // Extra intensity gate for the systematic-absence test: a reflection also has to reach this + // resolution-normalised intensity E^2 = I / (shell) to count as violating a predicted absence. + // The error model can under-estimate sigma on weak axial reflections and fake a high I/sigma, so a + // reflection at a few percent of the shell-mean intensity is judged absent regardless of its sigma + // (e.g. the monoclinic 2_1: 0k0-odd at ~1% of 0k0-even). 0 disables the gate (I/sigma only). + double present_e_squared = 0.3; + // A candidate's absence conditions are accepted when at most this fraction of the reflections // it predicts absent are in fact strongly present. double max_absent_violation_fraction = 0.10; diff --git a/tests/SearchSpaceGroupTest.cpp b/tests/SearchSpaceGroupTest.cpp index c330f5a5..618242c4 100644 --- a/tests/SearchSpaceGroupTest.cpp +++ b/tests/SearchSpaceGroupTest.cpp @@ -146,4 +146,48 @@ TEST_CASE("SearchSpaceGroup detects synthetic space groups") { CHECK(std::find(accepted.begin(), accepted.end(), tc.expected_short_name) != accepted.end()); } } +} + +// Regression: a real screw axis whose systematically-absent reflections carry a genuinely weak +// intensity but an UNDER-estimated sigma (so their I/sigma clears the "present" cut) must still be +// found. Reproduces the monoclinic 2_1 miss on EP_cs_02-10 / EP_cs_01-17, where the merged sigmas on +// the 0k0-odd reflections were ~2x too small and faked screw-axis violations. The E^2 intensity gate +// (present_e_squared) is what keeps those reflections classified absent. +TEST_CASE("SearchSpaceGroup finds a screw axis despite under-estimated sigmas on absent reflections") { + const gemmi::SpaceGroup& sg = gemmi::get_spacegroup_by_name("P 1 21 1"); + auto merged = GenerateMergedReflectionsForSpaceGroup(sg, 18); + + // Every systematically-absent (0k0, k odd) reflection: small-but-nonzero intensity (~2% of a + // normal reflection) with a far-too-small sigma, so I/sigma ~ 27 fakes a "present" reflection. + const gemmi::GroupOps gops = sg.operations(); + int absent_count = 0; + for (auto& r : merged) { + const gemmi::Op::Miller hkl{{r.h, r.k, r.l}}; + if (gops.is_systematically_absent(hkl)) { + r.I = 8.0f; + r.sigma = 0.3f; + ++absent_count; + } + } + REQUIRE(absent_count >= 8); // enough predicted-absent reflections to be trusted + + SearchSpaceGroupOptions opt; + opt.merge_friedel = true; + + SECTION("intensity gate on (default): screw recovered") { + const auto result = SearchSpaceGroup(merged, opt); + INFO(SearchSpaceGroupResultToText(result)); + REQUIRE(result.best_space_group.has_value()); + CHECK(result.best_space_group->short_name() == "P21"); + } + + SECTION("intensity gate off (I/sigma only): the screw is missed") { + // Documents the failure the gate fixes: with I/sigma alone the too-small sigmas fake + // violations and the search falls back to the symmorphic group. + opt.present_e_squared = 0.0; + const auto result = SearchSpaceGroup(merged, opt); + INFO(SearchSpaceGroupResultToText(result)); + REQUIRE(result.best_space_group.has_value()); + CHECK(result.best_space_group->short_name() == "P2"); + } } \ No newline at end of file