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