diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 6e7b86c7a..c404441ce 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -110,14 +110,25 @@ bool ValidationEvidenceBeatsChance(const ValidationSpotEvidence &e) { return ValidationEvidenceBeatsChance(e, std::erfc(SPOT_BUDGET_SIGNIFICANCE_Z / std::sqrt(2.0))); } -bool ValidationEvidenceBeatsChance(const ValidationSpotEvidence &e, double significance) { +double ValidationEvidenceLog10Tail(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; + // Summed in logs: a real lattice's tail is far below the smallest double. + std::vector log_terms; 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 < significance; + log_terms.push_back(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)); + if (log_terms.empty()) + return 0.0; + const double top = *std::max_element(log_terms.begin(), log_terms.end()); + double sum = 0.0; + for (const double t : log_terms) + sum += std::exp(t - top); + return (top + std::log(sum)) / std::log(10.0); +} + +bool ValidationEvidenceBeatsChance(const ValidationSpotEvidence &e, double significance) { + return ValidationEvidenceLog10Tail(e) < std::log10(significance); } RotationScaleWalk WalkRotationScale(double first_fit, diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 1a49debf0..ab95b7efc 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -234,6 +234,8 @@ constexpr int VALIDATION_NULL_DISPLACEMENTS = 5; // is below `significance` - by default that of SPOT_BUDGET_SIGNIFICANCE_Z. bool ValidationEvidenceBeatsChance(const ValidationSpotEvidence &e); bool ValidationEvidenceBeatsChance(const ValidationSpotEvidence &e, double significance); +// The binomial tail ValidationEvidenceBeatsChance compares with its significance, as log10. +double ValidationEvidenceLog10Tail(const ValidationSpotEvidence &e); // The significance the run's no-crystal test (RejectChanceLattice) asks of the lattice the first // pass ended on - Z ~ 6 instead of the 0.001 the comparisons between candidate lattices use. The diff --git a/rugnux/RugnuxFirstPassAccept.cpp b/rugnux/RugnuxFirstPassAccept.cpp index 0ef4ec03b..55fd883f5 100644 --- a/rugnux/RugnuxFirstPassAccept.cpp +++ b/rugnux/RugnuxFirstPassAccept.cpp @@ -147,16 +147,21 @@ auto Rugnux::FirstPassRun::MeasureVerdict(const FirstPass &best, const std::vect void Rugnux::FirstPassRun::RejectChanceLattice(const FirstPass &best, const Verdict &v) const { const PooledEvidence &evidence = v.evidence; const double pooled = v.pooled, pooled_chance = v.pooled_chance; - // Near tie (log only): the distinct reflections against the count at which the binomial tail - // ValidationEvidenceBeatsChance tests reaches its significance, in the normal approximation. + // Near tie (log only): the exact binomial tail against the significance, both as log10. The + // noise is how far the tail moves for one binomial standard deviation of the count. { const double n = static_cast(evidence.reflections + evidence.reflections_by_chance); const double p = 1.0 / (1.0 + VALIDATION_NULL_DISPLACEMENTS); - const double sd = std::sqrt(n * p * (1.0 - p)); - const double bar = n * p + SPOT_BUDGET_SIGNIFICANCE_Z * sd; - NearTie("no-crystal test", fmt::format("{} distinct reflections", evidence.reflections), - fmt::format("bar {:.1f}", bar), static_cast(evidence.reflections) - bar, sd, - "binomial on distinct reflections"); + const double log10_tail = ValidationEvidenceLog10Tail(evidence); + PooledEvidence one_more = evidence; + one_more.reflections++; + const double per_reflection = log10_tail - ValidationEvidenceLog10Tail(one_more); + NearTie("no-crystal test", + fmt::format("{} distinct reflections vs {} by chance, log10 tail {:.2f}", + evidence.reflections, evidence.reflections_by_chance, log10_tail), + fmt::format("log10 significance {:.2f}", std::log10(NO_CRYSTAL_SIGNIFICANCE)), + std::log10(NO_CRYSTAL_SIGNIFICANCE) - log10_tail, + per_reflection * std::sqrt(n * p * (1.0 - p)), "exact binomial on distinct reflections, decades"); } if (!ValidationEvidenceBeatsChance(evidence, NO_CRYSTAL_SIGNIFICANCE)) { { diff --git a/tests/RotationScaleWalkTest.cpp b/tests/RotationScaleWalkTest.cpp index 9b03e33bd..617610485 100644 --- a/tests/RotationScaleWalkTest.cpp +++ b/tests/RotationScaleWalkTest.cpp @@ -158,3 +158,11 @@ TEST_CASE("ValidationEvidenceBeatsChance_NoCrystalSignificance", "[RotationScale // A weak real lattice on a few frames is still tens of decades past it. CHECK(ValidationEvidenceBeatsChance({.reflections = 115, .reflections_by_chance = 4}, NO_CRYSTAL_SIGNIFICANCE)); } + +TEST_CASE("ValidationEvidenceLog10Tail_Exact", "[RotationScale]") { + CHECK(ValidationEvidenceLog10Tail({.reflections = 4, .reflections_by_chance = 2}) == Catch::Approx(-2.0604).margin(1e-4)); + CHECK(ValidationEvidenceLog10Tail({.reflections = 5, .reflections_by_chance = 2}) == Catch::Approx(-2.6981).margin(1e-4)); + // Far below the smallest double, still finite. + CHECK(ValidationEvidenceLog10Tail({.reflections = 1354, .reflections_by_chance = 777}) == Catch::Approx(-509.691).margin(1e-3)); + CHECK(ValidationEvidenceLog10Tail({}) == 0.0); +}