Space-group search: decide a screw from the evidence, not from a count of absences

A screw's predicted-absent class was required to hold min_absent_observed = 8 reflections before the
screw could be claimed. That count is the wrong measure of evidence, and it is wrong in both
directions.

A screw extinguishes 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 it measures the geometry of the sweep. A monoclinic crystal whose 2-fold sits
7.6 deg from the spindle contributed six 0k0-odd reflections, every one of them measured between
-0.013 and 4e-5 of the shell mean with zero violations, against a 0k0 row averaging 1.44x the shell
mean - and was refused its 2_1 for being six rather than eight. XDS's own integration of the same
images finds seventeen of those reflections and every one of them is likewise dead.

The count is equally wrong the other way: a uniformly weak axial row produces no violations at all,
so with enough reflections on it a screw is claimed from no evidence whatsoever. The second new test
section demonstrates exactly that on the old gate.

Judge the class by how unlikely it would be if the screw did not exist. Under "no screw" the absent
class and the rest of its row are both Wilson-distributed with the same mean, so with each absent
intensity taken in units of its row's control mean, sum_u/(sum_u + n_control) follows Beta(n_absent,
n_control) exactly; the reported evidence is -log of that lower tail. The row's own strength cancels,
which is the property the count lacks, 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.

This follows POINTLESS (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 confidence fall away with the number of axial reflections instead of refusing outright
below a count. POINTLESS calibrates its null width from control transforms of non-axial reflections;
the Beta tail here is an analytic null in its place. XDS is not a reference for this: it "deliberately
avoids any test for the presence of screw axes as these tests would depend strongly on the
completeness of the data" (Kabsch, Acta Cryst D66, 133-144 (2010), section 6), so a screw axis in a
CORRECT.LP was supplied to it, not determined by it.

Measured over five probe crystals, genuine screw conditions read 34-800 nats and false ones - the
4_1/4_3 conditions of a cubic crystal that has no screw, whose predicted-absent class is STRONGER
than its control row - read -7 to -8.5. The bound is set at 20, in the gap, at p <= 2e-9: three
well-measured dead axial reflections clear it and two do not.

min_absent_observed keeps its job for CENTERING, where a count is a fair measure - that class is a
third to a half of every reflection in the data set and the bound is never binding on a centering
that exists.

The candidate table now prints the screw-absent count and this evidence in place of the two E^2
medians that were its raw ingredients, so a refusal can be read off the log.

Measured on the five probes: the monoclinic crystal above returns to P2_1 with every merge statistic
unchanged (R_meas 58.8 -> 58.7%, CC1/2 49.1 -> 49.3%, ISa 6.61 -> 6.59 - P2 and P2_1 share a point
group, so only the symbol and the absent reflections differ). The other four are untouched, space
group included, and the two-pass branch fingerprint (indexed frames, distance, mosaicity) is
identical on all five. The full battery has not been run.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

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.

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<AxialRow, double> row_median;
+        std::map<AxialRow, double> 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<int>(screw_absent_refl.size());
         int screw_violations = 0;
-        std::vector<double> 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<AxialRow, size_t> 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<int>(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) << "<I/s>abs" << std::setw(11) << "<I/s>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/<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 <I/s> 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 <I/s> 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();

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/<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 <I/s> 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 {

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());
-}
+}
+// 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");
+    }
+}