From a29edb69cfed5717f9b7c0a1df54cd754d5c0178 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 8 Oct 2026 21:00:40 +0200 Subject: [PATCH] rugnux: decide whether there is a lattice on distinct reflections, against every wrong angle The first pass's last gate - is the lattice this crystal's at all - compared the validation spots on the lattice with the median of five displaced-spindle-angle nulls, against the null's binomial noise floored at one spot. On the two in-house no-crystal controls it is fed a handful of genuine stray reflections from something rotating with the sample (the same detector positions, a pair symmetric about the beam 198 frames apart, in both sweeps, offset by the 15-frame start difference), plus one reflection near the spindle axis that stays in diffraction for ~46 frames. With the decisions taken from the first part of the sweep the 60 validation frames are 3 deg apart instead of 6, so that one reflection is counted on three validation frames, and the gate passed on 4 spots of 151 (null median 0) and 7 of 72, where the whole-sweep flow had seen 3 and 3. Now the evidence unit is the distinct Miller index - a reflection on several validation frames is one piece of evidence - and the test uses every displaced angle's count: with no crystal the recorded angles are one of six exchangeable arrangements, so the gate is the binomial tail of the recorded arrangement's share at p = 1/6, at the same 0.001 significance (SPOT_BUDGET_SIGNIFICANCE_Z). No threshold is added; the median and the one-spot floor go from the gate. The comparisons between lattices (ValidationEvidencePrefers) are unchanged. no-crystal control 1: 4 distinct vs 2 over 5 wrong angles, tail 0.009 -> no lattice (was P 1) no-crystal control 2: 5 distinct vs 2, tail 0.002 -> no lattice (was P 1) p.mtz md5 identical to rc175 7c9a6c7d1 (GPU) on 5nw5, 5reo, 7atg, 8v4o, 9qw8, 9rci, 5ky6, insu_I_x06da_low_isa, myob_x06da_powder_2; cytidine and myob_x06da_split pass. The weakest real lattice in the open/in-house/private battery logs has 117 spots against a null of 0. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi --- rugnux/Rugnux.cpp | 10 +++++++ rugnux/Rugnux.h | 18 +++++++++++++ rugnux/RugnuxFirstPass.cpp | 48 ++++++++++++++++++++------------- tests/RotationScaleWalkTest.cpp | 13 +++++++++ 4 files changed, 71 insertions(+), 18 deletions(-) diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index a271fe370..804f78bdd 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -101,6 +101,16 @@ bool ValidationEvidencePrefers(const ValidationSpotEvidence ¤t, const Vali > SPOT_BUDGET_SIGNIFICANCE_Z * std::sqrt(rate_var(candidate) + rate_var(current)); } +bool ValidationEvidenceBeatsChance(const ValidationSpotEvidence &e) { + const int64_t n = e.reflections + e.reflections_by_chance; + const double p = 1.0 / (1.0 + VALIDATION_NULL_DISPLACEMENTS); + double tail = 0.0; + for (int64_t j = e.reflections; j <= n; j++) + tail += std::exp(std::lgamma(n + 1.0) - std::lgamma(j + 1.0) - std::lgamma(n - j + 1.0) + + j * std::log(p) + (n - j) * std::log(1.0 - p)); + return tail < std::erfc(SPOT_BUDGET_SIGNIFICANCE_Z / std::sqrt(2.0)); +} + RotationScaleWalk WalkRotationScale(double first_fit, const std::function &index_at, const std::function(float)> &refit_at, diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 355e1c0b3..4df80b6b8 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -201,12 +201,30 @@ struct ProcessConfig { // the ice rings unless those are indexed), how many lie on the lattice, and how many still do at a // wrong spindle angle - the median over displaced angles, which is what chance and the per-frame // orientation polish give for free. +// +// Whether there is a lattice at all is asked of distinct reflections rather than spots: a reflection +// that stays in diffraction across several validation frames (one near the spindle axis, or any one +// when the frames are a few degrees apart) is one piece of evidence however many frames it lies on. +// `reflections` counts the distinct Miller indices the lattice takes; `reflections_by_chance` sums +// the same count over every displaced angle (VALIDATION_NULL_DISPLACEMENTS of them). struct ValidationSpotEvidence { int64_t spots = 0; int64_t on_lattice = 0; int64_t by_chance = 0; + int64_t reflections = 0; + int64_t reflections_by_chance = 0; }; +constexpr int VALIDATION_NULL_DISPLACEMENTS = 5; + +// Whether the lattice puts more distinct reflections on itself at the recorded spindle angles than +// a wrong angle does. With no crystal, the recorded angles are just one more of the +// 1 + VALIDATION_NULL_DISPLACEMENTS arrangements, so each of the reflections the arrangements take +// between them falls on the recorded one with probability 1 / (1 + VALIDATION_NULL_DISPLACEMENTS); +// the lattice is this crystal's when the recorded arrangement takes so many that the binomial tail +// is below the significance of SPOT_BUDGET_SIGNIFICANCE_Z. +bool ValidationEvidenceBeatsChance(const ValidationSpotEvidence &e); + // Whether `candidate` puts a larger share of its validation spots on its lattice, over and above what // a wrong spindle angle puts there, than `current` does - by more than the binomial noise of the two // shares (SPOT_BUDGET_SIGNIFICANCE_Z). Each is measured against its own null, so a denser lattice is diff --git a/rugnux/RugnuxFirstPass.cpp b/rugnux/RugnuxFirstPass.cpp index 599a59494..1b3014cfc 100644 --- a/rugnux/RugnuxFirstPass.cpp +++ b/rugnux/RugnuxFirstPass.cpp @@ -580,9 +580,12 @@ bool Rugnux::FirstPassRotationIndexing(PipelineLocals &p) { // the bar the real count is measured against, and it is measured rather than chosen, so no // spots-per-frame number survives anywhere in the test. auto pooled_on_lattice = [&](IndexAndRefine &idx, const RotationIndexerResult &r, - int shift) -> int64_t { + int shift) -> std::pair { idx.ForceRotationIndexerResult(r); prefetch_spots(validation); + const bool index_ice_rings = experiment_.GetIndexingSettings().GetIndexIceRings(); + std::mutex reflections_mutex; + std::set> reflections; std::atomic on_lattice = 0; ParallelFor(static_cast(validation.size()), std::min(validation.size(), config_.nthreads), [&](int i) { @@ -590,8 +593,13 @@ bool Rugnux::FirstPassRotationIndexing(PipelineLocals &p) { m.number = validation[(i + shift) % validation.size()]; m.spots = spot_cache.at(validation[i]); on_lattice += idx.CountSpotsOnLattice(m, validation_settings); + // The spots the count took are the ones it gave a Miller index (none is 0,0,0). + std::unique_lock ul(reflections_mutex); + for (const auto &sp : m.spots) + if ((index_ice_rings || !sp.ice_ring) && (sp.h != 0 || sp.k != 0 || sp.l != 0)) + reflections.insert({sp.h, sp.k, sp.l}); }); - return on_lattice.load(); + return {on_lattice.load(), static_cast(reflections.size())}; }; // The pooled evidence for one candidate lattice: the validation frames' non-ice spots, how many @@ -605,7 +613,7 @@ bool Rugnux::FirstPassRotationIndexing(PipelineLocals &p) { // unlucky displacement can do. auto pooled_evidence = [&](IndexAndRefine &idx, const RotationIndexerResult &r) { PooledEvidence e; - e.on_lattice = pooled_on_lattice(idx, r, 0); + std::tie(e.on_lattice, e.reflections) = pooled_on_lattice(idx, r, 0); const bool index_ice_rings = experiment_.GetIndexingSettings().GetIndexIceRings(); for (const int ordinal : validation) for (const auto &sp : spot_cache.at(ordinal)) @@ -613,23 +621,22 @@ bool Rugnux::FirstPassRotationIndexing(PipelineLocals &p) { e.spots++; const int n = static_cast(validation.size()); std::vector chance; - for (int k = 1; k <= 5; k++) - chance.push_back(pooled_on_lattice(idx, r, std::max(1, k * n / 7))); + for (int k = 1; k <= VALIDATION_NULL_DISPLACEMENTS; k++) { + const auto [spots_on, reflections_on] = pooled_on_lattice(idx, r, std::max(1, k * n / 7)); + chance.push_back(spots_on); + e.reflections_by_chance += reflections_on; + } std::sort(chance.begin(), chance.end()); e.by_chance = chance[chance.size() / 2]; return e; }; - // This crystal's lattice is the one that puts materially more spots on itself than a wrong - // angle does. The excess is counted against the binomial noise of the null's own rate over the - // same spots, so the bar scales with the evidence in hand and carries no absolute count; the - // significance is the one-false-positive-in-a-thousand bar the spot-budget measurement is - // already set at. A null that indexed nothing is still credited with one spot of noise, as the - // next displacement could have found one. + // This crystal's lattice is the one that puts materially more reflections on itself than a + // wrong angle does (ValidationEvidenceBeatsChance): every displaced angle's count is used, so + // the bar carries the null's own noise, and it carries no absolute count; the significance is + // the one-false-positive-in-a-thousand bar the spot-budget measurement is already set at. auto beats_chance = [](const PooledEvidence &e) { - const double p = static_cast(e.by_chance) / static_cast(std::max(1, e.spots)); - const double sigma = std::sqrt(std::max(1.0, static_cast(e.spots) * p * (1.0 - p))); - return static_cast(e.on_lattice - e.by_chance) > SPOT_BUDGET_SIGNIFICANCE_Z * sigma; + return ValidationEvidenceBeatsChance(e); }; // Two lattices are COMPARED on the same quantity, ValidationEvidencePrefers: the share of the @@ -2268,10 +2275,13 @@ bool Rugnux::FirstPassRotationIndexing(PipelineLocals &p) { const auto &c = best.result->search_result.conventional.GetUnitCell(); throw NoLatticeFound(fmt::format("Two-pass rotation indexing found only a lattice that " "takes {:.1f}% of the validation spots against the " - "{:.1f}% a wrong spindle angle takes - it is not this " - "crystal's lattice. It was {}-centred {}, " + "{:.1f}% a wrong spindle angle takes ({} distinct " + "reflections against {} over {} wrong angles) - it is " + "not this crystal's lattice. It was {}-centred {}, " "{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}.{}", 100.0 * pooled, 100.0 * pooled_chance, + evidence.reflections, evidence.reflections_by_chance, + VALIDATION_NULL_DISPLACEMENTS, best.result->search_result.centering, gemmi::crystal_system_str(best.result->search_result.system), c.a, c.b, c.c, c.alpha, c.beta, c.gamma, next_steps)); @@ -2300,9 +2310,11 @@ bool Rugnux::FirstPassRotationIndexing(PipelineLocals &p) { logger.Info("First-pass spot finding: {} frames in {} batches, {:.2f} s", spot_cache.size(), prefetch_batches, prefetch_time_s); logger.Info("Two-pass rotation indexing found lattice (scheme '{}': {}/{} validation frames, " - "{}/{} validation spots = {:.1f}% against {:.1f}% at a wrong spindle angle)", + "{}/{} validation spots = {:.1f}% against {:.1f}% at a wrong spindle angle; {} " + "distinct reflections against {} over {} wrong angles)", best.name, best.score, static_cast(validation.size()), - evidence.on_lattice, evidence.spots, 100.0 * pooled, 100.0 * pooled_chance); + evidence.on_lattice, evidence.spots, 100.0 * pooled, 100.0 * pooled_chance, + evidence.reflections, evidence.reflections_by_chance, VALIDATION_NULL_DISPLACEMENTS); { const auto &c = best.result->search_result.conventional.GetUnitCell(); Note(fmt::format("lattice {} {}: {:.2f} {:.2f} {:.2f} {:.1f} {:.1f} {:.1f}", diff --git a/tests/RotationScaleWalkTest.cpp b/tests/RotationScaleWalkTest.cpp index 248bc9ecd..7fc4e1abb 100644 --- a/tests/RotationScaleWalkTest.cpp +++ b/tests/RotationScaleWalkTest.cpp @@ -134,3 +134,16 @@ TEST_CASE("ChooseLogScaleSmoothing_SmoothsNoiseFollowsSignal", "[RotationScale]" CHECK(wave[47] == Catch::Approx(step[47]).margin(0.01)); CHECK(l_wave < l_noise); } + +TEST_CASE("ValidationEvidenceBeatsChance_ExchangeableArrangements", "[RotationScale]") { + // Four reflections taken at the recorded angles and none at any wrong one: (1/6)^4 < 0.001. + CHECK(ValidationEvidenceBeatsChance({.reflections = 4, .reflections_by_chance = 0})); + CHECK_FALSE(ValidationEvidenceBeatsChance({.reflections = 3, .reflections_by_chance = 0})); + // The same four with two more spread over the wrong angles is what chance does once in ~100. + CHECK_FALSE(ValidationEvidenceBeatsChance({.reflections = 4, .reflections_by_chance = 2})); + CHECK_FALSE(ValidationEvidenceBeatsChance({})); + // A real lattice: thousands against a handful. + CHECK(ValidationEvidenceBeatsChance({.reflections = 3000, .reflections_by_chance = 250})); + // As many at the recorded angles as at an average wrong one is no lattice. + CHECK_FALSE(ValidationEvidenceBeatsChance({.reflections = 100, .reflections_by_chance = 500})); +}