diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index b4637fad..191546b6 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -3,6 +3,7 @@ ### 1.0.0-rc.161 This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use. +* Space-group search: a screw axis is now claimed from how decisively its predicted-absent reflections are weaker than the rest of their own axial row, instead of from a minimum count of them, so a screw survives a sweep that recorded few axial reflections and is refused on a row too weak to decide either way. * rugnux: The error-model **a** and **b** are reported in XDS's convention, and `_reflns.jfjoch_diffrn_ISa` now carries the whole-range `1/sqrt(a*b)` that XDS's ISa denotes; the strong-reflection asymptote moves to `_reflns.jfjoch_diffrn_ISa_asymptotic`. **A file written by an earlier version carries the asymptote under the old name.** * Bragg integration: the background ring's outer radius default changes from 10 px to **13 px**, which roughly doubles the pixels behind each background estimate; the signal disk is unchanged. * Scaling: when too few reflections are strong enough to constrain the error model's systematic term **b**, it is now held at zero and **ISa is reported as unmeasured** rather than as the spurious value the fit would otherwise extrapolate. diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index 8f264ca3..312fba78 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -95,6 +95,26 @@ namespace { return v[mid]; } + // How unlikely the predicted-absent class would be if the screw did not exist, in nats. + // + // Under "no screw" the absent class and the rest of its axial row are both Wilson-distributed with + // the SAME mean, so with each absent intensity expressed in units of its row's control mean, the + // fraction T = sum_u / (sum_u + n_control) follows Beta(n_absent, n_control) exactly. The row's own + // strength cancels out of it - which is the property a violation count does not have, and the + // reason a uniformly weak axial row decides nothing here instead of deciding "absent". Returns + // -log of that Beta lower tail. + // + // Only the leading term of the regularized incomplete beta is kept. It is exact as T -> 0, which is + // where a screw is claimed, and dropping the (1-T)^n_control factor only ever UNDER-states the + // evidence, which is the safe direction for a test that has to clear a bound. + double ScrewAbsenceEvidence(double sum_u, int n_absent, int n_control) { + if (n_absent <= 0 || n_control <= 0) + return 0.0; + const double a = n_absent, b = n_control; + const double T = std::max(sum_u / (sum_u + b), 1e-300); + return -(a * std::log(T) + std::lgamma(a + b) - std::lgamma(a + 1) - std::lgamma(b)); + } + std::string FormatDouble(double v, int decimals) { std::ostringstream o; o << std::fixed << std::setprecision(decimals) << v; @@ -807,14 +827,22 @@ SearchSpaceGroupResult SearchSpaceGroup( // observed only at l=4n, strongly, would read its l=4n+2 reflections as absent and be reported as // 4_1/4_3. Below this count the row falls back to the plain cut, i.e. no rescue. constexpr size_t MIN_ROW_CONTROL_REFLECTIONS = 3; - std::map row_median; + std::map row_median, row_mean; for (auto& [row, esq] : row_present_esq) - if (esq.size() >= MIN_ROW_CONTROL_REFLECTIONS) + if (esq.size() >= MIN_ROW_CONTROL_REFLECTIONS) { + row_mean[row] = std::accumulate(esq.begin(), esq.end(), 0.0) / esq.size(); 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; + // The absent intensities in units of their own row's control mean, and the rows that supplied + // one - the ingredients of ScrewAbsenceEvidence. A row with no control class of its own cannot + // say whether its absences are weak or its whole row is, so it contributes no evidence either + // way rather than being judged against the shell mean. + double screw_sum_u = 0.0; + int screw_scored = 0; + std::map screw_control_rows; 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); @@ -825,14 +853,20 @@ SearchSpaceGroupResult SearchSpaceGroup( if (a.i_over_sigma > present_cut && (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); + const auto mit = row_mean.find(a.row); + if (mit != row_mean.end() && mit->second > 0.0) { + screw_sum_u += std::max(0.0, a.e_squared) / mit->second; + ++screw_scored; + screw_control_rows[a.row] = row_present_esq.at(a.row).size(); + } } + int screw_control = 0; + for (const auto& [row, count] : screw_control_rows) + screw_control += static_cast(count); 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); + s.centering_absent = centering_absent; + s.screw_absent = screw_absent; + s.screw_absence_evidence = ScrewAbsenceEvidence(screw_sum_u, screw_scored, screw_control); if (s.absent_observed > 0) s.absent_mean_i_over_sigma = absent_sum / s.absent_observed; @@ -895,8 +929,18 @@ SearchSpaceGroupResult SearchSpaceGroup( auto net_absent = [](const SpaceGroupCandidateScore& s) { return s.absent_observed - s.absent_violations; }; + // A candidate also needs enough EVIDENCE behind its absences before they may be claimed, and the two + // kinds of absence need different measures of it. A centering class is a third to a half of every + // reflection in the data set, so a count is a fair measure and min_absent_observed is never the + // binding constraint. A screw class is a handful of axial reflections BY CONSTRUCTION - one row of + // reciprocal space, often lying near the spindle where a rotation sweep records least - and there a + // count measures the sweep's geometry, not the evidence: six axial reflections measured at zero + // against a row that averages 1.4x the shell mean settle the question, while twenty uniformly weak + // ones settle nothing. So the screw class is judged by ScrewAbsenceEvidence instead. auto eligible = [&](const SpaceGroupCandidateScore& s) { - return s.consistent && (s.absent_observed == 0 || s.absent_observed >= opt.min_absent_observed); + return s.consistent + && (s.centering_absent == 0 || s.centering_absent >= opt.min_absent_observed) + && (s.screw_absent == 0 || s.screw_absence_evidence >= opt.min_screw_absence_evidence); }; std::sort(result.candidates.begin(), result.candidates.end(), [&](const SpaceGroupCandidateScore& a, const SpaceGroupCandidateScore& b) { @@ -958,7 +1002,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(9) << "screw" << std::setw(11) << "screw evid" << std::setw(6) << "OK" << "\n"; const size_t count = std::min(max_candidates_to_print, result.candidates.size()); @@ -969,14 +1013,16 @@ 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(9) << c.screw_absent + << std::setw(11) << std::fixed << std::setprecision(1) << c.screw_absence_evidence << 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.\n"; + " screw = how many of those lie on an axial row, i.e. are extinguished by a screw rather than\n" + " by the centering; screw evid = how unlikely that class would be if the screw did not exist,\n" + " judged against the rest of its own axial row (in nats - a real screw reads tens to hundreds,\n" + " a false one at or below zero). The columns are the centering evidence; they say little\n" + " about screws, because the merged sigma shrinks with I on absent and present alike.\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 438b178b..e240f3b9 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.h +++ b/image_analysis/scale_merge/SearchSpaceGroup.h @@ -43,14 +43,14 @@ struct SpaceGroupCandidateScore { gemmi::SpaceGroup space_group; int absent_observed = 0; // observed reflections this SG predicts systematically absent int absent_violations = 0; // of those, how many are nonetheless strongly present + int centering_absent = 0; // of absent_observed, the ones extinguished by the CENTERING + int screw_absent = 0; // ... and the ones extinguished by a SCREW (i.e. on axial rows) + // How unlikely the screw-absent class would be if the screw did not exist, in nats - the quantity + // that decides whether a screw may be claimed at all (see min_screw_absence_evidence). 0 when there + // are no screw absences, or when no axial row carried a control class to judge them against. + double screw_absence_evidence = 0.0; 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 error model's - // intensity-proportional term, so they read much 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) }; @@ -240,9 +240,40 @@ struct SearchSpaceGroupOptions { // pushed genuinely-absent reflections over I/sigma>3 (a true R3 was lost at 13.5% violations). double max_absent_present_ratio = 0.5; - // Need at least this many observed reflections in the predicted-absent class before a - // screw/centering is claimed (guards against deciding from a handful of reflections). + // Need at least this many observed reflections in the CENTERING-absent class before a centering is + // claimed (guards against deciding from a handful of reflections). A real centering extinguishes a + // third to a half of every reflection in the data set, so this is never the binding constraint on a + // centering that exists - it only refuses one the data barely sampled. + // + // It does NOT apply to screws. A screw's absent class is one row of reciprocal space by + // construction, and that row is often the one a rotation sweep records least: it lies near the + // spindle, where the blind cusp maps onto itself and symmetry cannot fill it in. Counting it + // measures the geometry of the sweep, not the strength of the evidence. Applied to screws this + // bound cost a monoclinic crystal its 2_1 - six 0k0-odd reflections, every one of them measured at + // |E^2| <= 0.013 with zero violations, against a 0k0 row averaging 1.44x the shell mean, refused + // for being six rather than eight. Screws answer to min_screw_absence_evidence below instead. int min_absent_observed = 8; + + // Evidence, in nats, that a screw's predicted-absent class really is absent, required before the + // screw may be claimed. In the spirit of the POINTLESS zone test (Evans, Acta Cryst D67, 282-292 + // (2011), Appendix A3), which likewise scores an absence against the rest of its own axial row + // rather than against a global mean or a fixed cut, and likewise lets the confidence fall away with + // the number of axial reflections instead of refusing outright below a count. + // + // The statistic (see ScrewAbsenceEvidence) is the -log Beta tail of the absent class against its + // row's control class. Two properties are what a violation count lacks: the row's own strength + // cancels, so a uniformly weak axial row decides nothing rather than reading as "absent"; and the + // scale is set by the number of reflections, so few-but-decisive and many-but-marginal are told + // apart. It is sigma-free by design - the merged sigma carries the error model's + // intensity-proportional term and so shrinks with I, reading much the same on an absent reflection + // as on a present one (the columns of the candidate table show this directly). + // + // Measured over the probe crystals, genuine screw conditions read 34-800 nats and the false ones + // (the 4_1/4_3 conditions of a cubic crystal that has no screw at all, whose predicted-absent class + // is STRONGER than its control row) read -7 to -8.5. The gap is ~40 nats wide and this bound sits + // in it, at p <= 2e-9 - loose enough that three well-measured dead axial reflections clear it, + // tight enough that two do not. + double min_screw_absence_evidence = 20.0; }; struct SearchSpaceGroupResult { diff --git a/tests/SearchSpaceGroupTest.cpp b/tests/SearchSpaceGroupTest.cpp index bcacbb70..3d3e7ba3 100644 --- a/tests/SearchSpaceGroupTest.cpp +++ b/tests/SearchSpaceGroupTest.cpp @@ -234,4 +234,63 @@ TEST_CASE("SearchSpaceGroup finds a screw axis whose absent class is weak only w 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 +} +// Regression: a screw's predicted-absent class is one row of reciprocal space, and that row is often +// the one a rotation sweep records least - it lies near the spindle, where the blind cusp maps onto +// itself and symmetry cannot fill it in. Counting the class therefore measures the geometry of the +// sweep, not the strength of the evidence, and a count gate refused a monoclinic crystal its 2_1 for +// having six 0k0-odd reflections rather than eight, every one of them measured at a thousandth of the +// row beside them. The class is judged by ScrewAbsenceEvidence instead, which reads the contrast +// against the row - so few-but-decisive is accepted and many-but-marginal is not. +TEST_CASE("SearchSpaceGroup weighs a screw's absences by evidence, not by how many were recorded") { + const gemmi::SpaceGroup& sg = gemmi::get_spacegroup_by_name("P 1 21 1"); + const gemmi::GroupOps gops = sg.operations(); + + SearchSpaceGroupOptions opt; + opt.merge_friedel = true; + + SECTION("five decisive absences, below min_absent_observed: the screw is still found") { + auto merged = GenerateMergedReflectionsForSpaceGroup(sg, 18); + // Keep five of the 0k0-odd reflections, at a thousandth of their row, and drop the rest - as a + // sweep along the 2-fold does, leaving too few to satisfy a count but plenty to decide. + int kept = 0; + std::erase_if(merged, [&](MergedReflection& r) { + if (!gops.is_systematically_absent(gemmi::Op::Miller{{r.h, r.k, r.l}})) + return false; + if (kept >= 5) + return true; + ++kept; + r.I = 0.5; + return false; + }); + REQUIRE(kept == 5); + REQUIRE(kept < opt.min_absent_observed); + + 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("a uniformly weak axial row decides nothing, however many absences it holds") { + // The whole 0k0 row badly measured: the predicted-absent reflections are weak, but so is the + // rest of their row, so there is no contrast and no screw to claim. A violation count cannot + // see this - nothing on the row clears an absolute cut, so it reads zero violations and, with + // enough reflections to satisfy the count, would claim the 2_1 from no evidence at all. + auto merged = GenerateMergedReflectionsForSpaceGroup(sg, 18); + int absent_on_row = 0; + for (auto& r : merged) { + if (r.h != 0 || r.l != 0) + continue; + const bool absent = gops.is_systematically_absent(gemmi::Op::Miller{{r.h, r.k, r.l}}); + r.I = absent ? 4.0 : 5.0; + absent_on_row += absent ? 1 : 0; + } + REQUIRE(absent_on_row >= opt.min_absent_observed); + + const auto result = SearchSpaceGroup(merged, opt); + INFO(SearchSpaceGroupResultToText(result)); + REQUIRE(result.best_space_group.has_value()); + CHECK(result.best_space_group->short_name() == "P2"); + } +}