diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index 9d108de5..fd03e6fc 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -72,6 +73,28 @@ namespace { return false; } + // The reciprocal-space ROW a reflection lies on: its direction, reduced by the gcd and + // sign-canonicalised, so 0,0,l and 0,0,-l are one row and h,0,0 is a different one. Screw + // absences are judged against the other reflections of their own row (see below). + using AxialRow = std::array; + + AxialRow RowOf(int h, int k, int l) { + const int g = std::gcd(std::gcd(std::abs(h), std::abs(k)), std::abs(l)); + if (g > 0) { h /= g; k /= g; l /= g; } + if (std::make_tuple(-h, -k, -l) < std::make_tuple(h, k, l)) + return {-h, -k, -l}; + return {h, k, l}; + } + + // Median of an unordered set (reordered in place); 0 for an empty set. + double MedianOf(std::vector& v) { + if (v.empty()) + return 0.0; + const size_t mid = v.size() / 2; + std::nth_element(v.begin(), v.begin() + mid, v.end()); + return v[mid]; + } + std::string FormatDouble(double v, int decimals) { std::ostringstream o; o << std::fixed << std::setprecision(decimals) << v; @@ -646,8 +669,13 @@ SearchSpaceGroupResult SearchSpaceGroup( // screw axes (e.g. I4_132 on I432 data). int centering_absent = 0, centering_violations = 0; double centering_absent_sum = 0; - int screw_absent = 0, screw_violations = 0; int present_strong = 0; + // A screw axis extinguishes only the reflections that lie ON it, so its absent class and the + // rest of that same axial row are collected apart from the general reflections and judged + // against each other, ROW BY ROW (see screw_e_squared below). + struct ScrewAbsent { AxialRow row; double e_squared; double i_over_sigma; }; + std::vector screw_absent_refl; + std::map> row_present_esq; for (size_t i = 0; i < n; ++i) { if (!pass_absence[i]) @@ -669,15 +697,55 @@ SearchSpaceGroupResult SearchSpaceGroup( } else if (gops.is_systematically_absent(hkl)) { s.absent_observed += 1; absent_sum += IoverSigma[i]; - screw_absent += 1; - if (present) { s.absent_violations += 1; screw_violations += 1; } + screw_absent_refl.push_back({RowOf(H[i], K[i], L[i]), Esq[i], IoverSigma[i]}); } else { present_n += 1; present_sum += IoverSigma[i]; if (present) present_strong += 1; + // A non-identity rotation of the group maps this reflection to itself, i.e. it lies on + // a rotation axis - the control class for that axis's screw absences. + if (gops.epsilon_factor_without_centering(hkl) > 1) + row_present_esq[RowOf(H[i], K[i], L[i])].push_back(Esq[i]); } } + // "Too strong to be systematically absent" is judged, for a SCREW, against the axial row the + // screw constrains rather than against the shell mean over all reflections. An axial row can be + // far stronger than an average reflection, and (shell) falls with resolution while a + // systematically-absent reflection keeps a small non-decaying residual (background / profile + // leakage) - so at high resolution the plain E^2 cut turns those residuals into violations even + // though the reflections next to them in the same row are tens of times stronger. That cost a + // real P4_1 2_1 2 crystal its 4_1: 18 of its 47 absent 00l crossed the cut, all beyond 3.7 A, + // at 1-2% of the l=4n reflections beside them. Scaling by the row's own median E^2 removes the + // resolution dependence; the scale is floored at 1 so a row weaker than average keeps the plain + // cut, which makes this a rescue only - a screw can be recovered by it, never lost. + // + // Row by row, not pooled: a 4_1 along c and a 2_1 along a are separate conditions with separate + // control rows, and on the same crystal one row can be 20x an average reflection while another + // is half of one. Pooling them lets the weak rows set the threshold for the strong one and the + // rescue never fires (that very crystal pooled to a row median of 0.7 and stayed at P42_12). + std::map row_median; + for (auto& [row, esq] : row_present_esq) + row_median[row] = MedianOf(esq); + + const int screw_absent = static_cast(screw_absent_refl.size()); + int screw_violations = 0; + std::vector screw_absent_esq, screw_control_esq; + for (const auto& a : screw_absent_refl) { + const auto it = row_median.find(a.row); + const double row_scale = it == row_median.end() ? 1.0 : std::max(1.0, it->second); + if (a.i_over_sigma > opt.present_i_over_sigma && + (opt.present_e_squared <= 0.0 || a.e_squared > opt.present_e_squared * row_scale)) + ++screw_violations; + screw_absent_esq.push_back(a.e_squared); + if (it != row_median.end()) + screw_control_esq.push_back(it->second); + } + s.absent_violations += screw_violations; + // Reported evidence: the absent class and the rows it is judged against (see the header). + s.screw_absent_median_e_squared = MedianOf(screw_absent_esq); + s.screw_row_median_e_squared = MedianOf(screw_control_esq); + if (s.absent_observed > 0) s.absent_mean_i_over_sigma = absent_sum / s.absent_observed; if (present_n > 0) @@ -691,7 +759,8 @@ SearchSpaceGroupResult SearchSpaceGroup( // reflections randomly clear I/sigma>3 to trip the 10% bound though the class is 3-4x weaker (a // true R3 at 13.5% violations, absent 1.7 vs present 6.0). The mean is well-determined here // because a centering-absent class holds a third-to-half of all reflections. Screws keep the - // count gate: their predicted-absent class is a handful of axial reflections, too few to average. + // count gate: their predicted-absent class is a handful of axial reflections, too few to average + // - what they get instead is a row-relative threshold for counting a violation at all. const double present_mean = present_n > 0 ? present_sum / present_n : 0.0; const double centering_absent_mean = centering_absent > 0 ? centering_absent_sum / centering_absent : 0.0; @@ -796,6 +865,7 @@ std::string SearchSpaceGroupResultToText(const SearchSpaceGroupResult& result, os << " " << std::setw(10) << std::left << "SG" << std::right << std::setw(9) << "absent" << std::setw(7) << "viol" << std::setw(11) << "abs" << std::setw(11) << "pres" + << std::setw(10) << "E2 screw" << std::setw(9) << "E2 row" << std::setw(6) << "OK" << "\n"; const size_t count = std::min(max_candidates_to_print, result.candidates.size()); @@ -806,8 +876,14 @@ std::string SearchSpaceGroupResultToText(const SearchSpaceGroupResult& result, << std::setw(9) << c.absent_observed << std::setw(7) << c.absent_violations << std::setw(11) << std::fixed << std::setprecision(2) << c.absent_mean_i_over_sigma << std::setw(11) << std::fixed << std::setprecision(2) << c.present_mean_i_over_sigma + << std::setw(10) << std::fixed << std::setprecision(3) << c.screw_absent_median_e_squared + << std::setw(9) << std::fixed << std::setprecision(3) << c.screw_row_median_e_squared << std::setw(6) << (c.consistent ? "yes" : "no") << "\n"; } + os << " absent/viol = reflections the group predicts absent, and how many are nonetheless present.\n" + " E2 screw / E2 row = median I/(shell) of the SCREW-absent reflections and of the rest of\n" + " their axial rows - a real screw leaves the first far below the second. The columns are\n" + " the centering evidence; they say little about screws (the merged sigma floor pins I/sigma).\n"; if (result.best_space_group.has_value()) { os << "Best space group: " << result.best_space_group->short_name(); diff --git a/image_analysis/scale_merge/SearchSpaceGroup.h b/image_analysis/scale_merge/SearchSpaceGroup.h index ad630672..e8b93553 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.h +++ b/image_analysis/scale_merge/SearchSpaceGroup.h @@ -45,6 +45,12 @@ struct SpaceGroupCandidateScore { int absent_violations = 0; // of those, how many are nonetheless strongly present double absent_mean_i_over_sigma = 0.0; double present_mean_i_over_sigma = 0.0; + // Screw evidence, as median E^2 = I/(shell): the reflections a screw axis extinguishes against + // the rest of that same axial row. This is the comparison the screw test makes, and the only one + // that means anything for it - the I/sigma means above are dominated by the merged-sigma floor + // (sigma >= b|I| pins I/sigma at ISa), so they read the same for absent and present alike. + double screw_absent_median_e_squared = 0.0; + double screw_row_median_e_squared = 0.0; bool consistent = false; // absent class confirmed weak (few violations) bool selected = false; // chosen result (or its enantiomorph) }; @@ -181,6 +187,10 @@ struct SearchSpaceGroupOptions { // 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). + // + // For a SCREW the threshold is this fraction of the median E^2 of the axial row the screw + // constrains (floored at the plain value), not of an average reflection at that resolution - an + // axial row is often far stronger than the shell mean, and only the row is a fair comparison. double present_e_squared = 0.3; // A candidate's SCREW/glide absence conditions are accepted when at most this fraction of the diff --git a/tests/SearchSpaceGroupTest.cpp b/tests/SearchSpaceGroupTest.cpp index 1fbcc5d9..bcacbb70 100644 --- a/tests/SearchSpaceGroupTest.cpp +++ b/tests/SearchSpaceGroupTest.cpp @@ -190,4 +190,48 @@ TEST_CASE("SearchSpaceGroup finds a screw axis despite under-estimated sigmas on REQUIRE(result.best_space_group.has_value()); CHECK(result.best_space_group->short_name() == "P2"); } +} + +// Regression: the E^2 gate above compares a reflection to the mean of its RESOLUTION SHELL, which +// falls off with resolution, while a systematically-absent reflection keeps a small non-decaying +// residual (background / profile leakage). On a crystal whose axial rows are much stronger than an +// average reflection, that turns the high-resolution residuals into screw-axis violations and the +// screw is lost, although the reflections beside them in the same row are tens of times stronger. +// A tetragonal 42_12 case failed exactly this way (18 of 47 absent 00l over the cut, all beyond +// 3.7 A, at 1-2% of the l=4n reflections next to them). The threshold is therefore taken relative to +// the axial row the screw constrains, not to the shell. +TEST_CASE("SearchSpaceGroup finds a screw axis whose absent class is weak only within its own row") { + const gemmi::SpaceGroup& sg = gemmi::get_spacegroup_by_name("P 43 21 2"); + auto merged = GenerateMergedReflectionsForSpaceGroup(sg, 12); + + // Axial rows 40x stronger than a general reflection, and an absent class carrying ~2% of its own + // row - but half of a general reflection, so a threshold set against the shell calls every one of + // them a violation while a threshold set against the row calls none. + const gemmi::GroupOps gops = sg.operations(); + int absent_on_axis = 0; + for (auto& r : merged) { + const gemmi::Op::Miller hkl{{r.h, r.k, r.l}}; + if (gops.epsilon_factor_without_centering(hkl) <= 1) + continue; + if (gops.is_systematically_absent(hkl)) { + r.I = 300.0f; + r.sigma = 1.0f; + ++absent_on_axis; + } else { + r.I *= 40.0f; + } + } + REQUIRE(absent_on_axis >= 8); + + SearchSpaceGroupOptions opt; + opt.merge_friedel = true; + + const auto result = SearchSpaceGroup(merged, opt); + INFO(SearchSpaceGroupResultToText(result)); + REQUIRE(result.best_space_group.has_value()); + // P4_1 2_1 2 and P4_3 2_1 2 are enantiomorphs and indistinguishable from intensities. + std::vector accepted{result.best_space_group->short_name()}; + for (const auto& alt : result.alternatives) + accepted.push_back(alt.short_name()); + CHECK(std::find(accepted.begin(), accepted.end(), "P43212") != accepted.end()); } \ No newline at end of file