diff --git a/rugnux/CMakeLists.txt b/rugnux/CMakeLists.txt index d21fa2c66..3748f9808 100644 --- a/rugnux/CMakeLists.txt +++ b/rugnux/CMakeLists.txt @@ -9,6 +9,7 @@ ADD_LIBRARY(JFJochRugnux STATIC RugnuxPreScan.cpp RugnuxPasses.cpp RugnuxImagePass.cpp + RugnuxFarEnd.cpp RugnuxFirstPass.h RugnuxFirstPass.cpp RugnuxFirstPassRescues.cpp diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 1b4aa19da..659039ccc 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -88,7 +88,8 @@ using namespace rugnux_internal; -bool ValidationEvidencePrefers(const ValidationSpotEvidence ¤t, const ValidationSpotEvidence &candidate) { +std::pair ValidationEvidenceDifference(const ValidationSpotEvidence ¤t, + const ValidationSpotEvidence &candidate) { const auto excess = [](const ValidationSpotEvidence &e) { return static_cast(e.on_lattice - e.by_chance) / static_cast(std::max(1, e.spots)); }; @@ -97,8 +98,12 @@ bool ValidationEvidencePrefers(const ValidationSpotEvidence ¤t, const Vali const double p = static_cast(e.on_lattice) / n; return p * (1.0 - p) / n; }; - return excess(candidate) - excess(current) - > SPOT_BUDGET_SIGNIFICANCE_Z * std::sqrt(rate_var(candidate) + rate_var(current)); + return {excess(candidate) - excess(current), std::sqrt(rate_var(candidate) + rate_var(current))}; +} + +bool ValidationEvidencePrefers(const ValidationSpotEvidence ¤t, const ValidationSpotEvidence &candidate) { + const auto [difference, sigma] = ValidationEvidenceDifference(current, candidate); + return difference > SPOT_BUDGET_SIGNIFICANCE_Z * sigma; } bool ValidationEvidenceBeatsChance(const ValidationSpotEvidence &e) { diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 8526efac0..0f2806303 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -239,6 +239,10 @@ bool ValidationEvidenceBeatsChance(const ValidationSpotEvidence &e); // shares (SPOT_BUDGET_SIGNIFICANCE_Z). Each is measured against its own null, so a denser lattice is // not credited for the spots it catches by accident. bool ValidationEvidencePrefers(const ValidationSpotEvidence ¤t, const ValidationSpotEvidence &candidate); +// What ValidationEvidencePrefers compares: the candidate's excess share minus the current one's, and +// the binomial standard error of that difference. +std::pair ValidationEvidenceDifference(const ValidationSpotEvidence ¤t, + const ValidationSpotEvidence &candidate); // The goniometer rotation scale - the factor by which the stage turned relative to the angles stored // in the file - walked to the fit's fixed point and decided on the whole sweep. The fit only sees the @@ -974,6 +978,9 @@ class Rugnux { // Its steps and the state they share (RugnuxFirstPass.h). struct FirstPassRun; void ImagePass(PipelineLocals &p); + // After the image pass of a rotation sweep: the decisions taken on the sweep's first part + // against what pass 2 measured on the rest of it. Log only (RugnuxFarEnd.cpp). + void FarEndCheck(PipelineLocals &p, const RotationIndexerResult &rot); // Scaling and merging, the space-group search, the analyses of the merge and the merged output. // True where the pass ends early without merging: a post-refinement probe, or a starved pass. bool ScaleMergeAndSymmetry(PipelineLocals &p); diff --git a/rugnux/RugnuxFarEnd.cpp b/rugnux/RugnuxFarEnd.cpp new file mode 100644 index 000000000..4d5dac496 --- /dev/null +++ b/rugnux/RugnuxFarEnd.cpp @@ -0,0 +1,255 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include +#include +#include +#include + +#include "Rugnux.h" +#include "RugnuxPipeline.h" +#include "../image_analysis/indexing/AnalyzeIndexing.h" // SPOT_BUDGET_SIGNIFICANCE_Z + +// The far-end check. LOG ONLY: it decides nothing and changes nothing the run writes. +// +// The decisions before pass 2 - the lattice, the cell, the rotation scale, the orientation the +// sweep is predicted from, the shadow mask - are taken on the first part of the sweep +// (PrepassEnd), so that pass 1 can run during collection. Pass 2 then integrates the whole sweep, +// and what it measured on the rest is evidence those decisions never saw. Here the sweep is cut +// into blocks of FAR_END_BLOCK_DEG, a handful of per-frame quantities pass 2 already has in memory +// are averaged per block, and each quantity's level over the rest of the sweep is compared with its +// level on the first part, the block means of each part being its samples (Welch's t): what a part's +// blocks scatter by - orientation, dose, counting noise alike - is that part's noise, so one odd +// block is not a shift. A quantity is out of line when its two-sided p, times the number of +// quantities compared, is below the significance the run tests everything else at +// (SPOT_BUDGET_SIGNIFICANCE_Z, 0.1 %): an unchanged sweep then says "hold" in all but one run in a +// thousand. No number in it is set by hand beyond the block width. +// +// The quantities, and the decision each one speaks to: +// - the share of frames that index on their own, and the share of each frame's spots on the +// lattice: the lattice (a cell that grows with the dose, a crystal that slips); +// - the cell scale the strong reflections ask for, d at the observed centroid over d predicted: +// the cell; +// - the residual rotation between each frame's refined orientation and the sweep's lattice turned +// to that frame, about the spindle and off it: the rotation scale and the orientation; +// - the background: the shadow mask (hardware that shadows more, or less, as the sample turns). + +namespace { + constexpr double FAR_END_BLOCK_DEG = 10.0; + + // P(|T| > t) for Student's t on an integer number of degrees of freedom, in closed form + // (Abramowitz and Stegun 26.7.3 and 26.7.4). + double TwoSidedP(double t, int dof) { + if (!std::isfinite(t)) + return std::isnan(t) ? 1.0 : 0.0; + const double theta = std::atan(std::abs(t) / std::sqrt(static_cast(dof))); + const double c2 = std::cos(theta) * std::cos(theta); + double sum = 1.0, term = 1.0; + if (dof % 2 == 1) { + for (int j = 1; 2 * j + 1 <= dof - 2; j++) + sum += term *= c2 * (2.0 * j) / (2.0 * j + 1.0); + const double a = dof == 1 ? theta : theta + std::sin(theta) * std::cos(theta) * sum; + return std::max(0.0, 1.0 - 2.0 * a / M_PI); + } + for (int j = 1; 2 * j <= dof - 2; j++) + sum += term *= c2 * (2.0 * j - 1.0) / (2.0 * j); + return std::max(0.0, 1.0 - std::sin(theta) * sum); + } + + struct Block { + double mean = NAN; + double se = NAN; + int frames = 0; + }; + + Block Summarise(const std::vector &v, int first, int last) { + Block b; + double sum = 0.0, sum2 = 0.0; + for (int i = first; i < last; i++) + if (std::isfinite(v[i])) { + sum += v[i]; + sum2 += v[i] * v[i]; + b.frames++; + } + if (b.frames < 2) + return b; + b.mean = sum / b.frames; + const double var = std::max(0.0, (sum2 - b.frames * b.mean * b.mean) / (b.frames - 1)); + b.se = std::sqrt(var / b.frames); + return b; + } +} + +void Rugnux::FarEndCheck(PipelineLocals &p, const RotationIndexerResult &rot) { + auto &logger = p.logger; + const int n = p.images_to_process; + const int early = p.prepass_images; + if (!rot.axis || !p.indexer) + return; + if (early >= n) { + logger.Info("far-end check: the decisions read the whole sweep"); + return; + } + const GoniometerAxis &gon = *rot.axis; + const double step = std::abs(gon.GetAngle_deg(1.0f) - gon.GetAngle_deg(0.0f)); + if (!(step > 0.0)) + return; + const int block = std::max(1, static_cast(std::lround(FAR_END_BLOCK_DEG / step))); + + std::vector indexed, spots, spots_indexed, bkg; + p.plots.GetPlotRaw(indexed, PlotType::IndexingRate, ""); + p.plots.GetPlotRaw(spots, PlotType::SpotCount, ""); + p.plots.GetPlotRaw(spots_indexed, PlotType::SpotCountIndexed, ""); + p.plots.GetPlotRaw(bkg, PlotType::BkgEstimate, ""); + const auto at = [](const std::vector &v, int i) { + return i < static_cast(v.size()) ? static_cast(v[i]) : NAN; + }; + + enum { INDEXED, ON_LATTICE, CELL, ABOUT_SPINDLE, OFF_SPINDLE, BACKGROUND, N_STATS }; + struct Stat { const char *what; const char *decision; double unit; const char *suffix; const char *tag; }; + const Stat stats[N_STATS] = { + {"frames indexing on their own", "the lattice", 100.0, "%", "indexed"}, + {"spots on the lattice", "the lattice", 100.0, "%", "on-lattice"}, + {"cell scale of the strong reflections (d observed / d predicted - 1)", "the cell", 100.0, "%", "cell"}, + {"orientation residual about the spindle", "the rotation scale", 1.0, " deg", "about-spindle"}, + {"orientation residual off the spindle", "the orientation", 1.0, " deg", "off-spindle"}, + {"background", "the shadow mask", 1.0, "", "background"}, + }; + std::vector> frame(N_STATS, std::vector(n, NAN)); + const Coord spindle = gon.GetAxis().Normalize(); + const auto &outcomes = p.indexer->GetIntegrationOutcome(); + for (int o = 0; o < n; o++) { + frame[INDEXED][o] = at(indexed, o); + if (at(spots, o) > 0.0) + frame[ON_LATTICE][o] = at(spots_indexed, o) / at(spots, o); + frame[BACKGROUND][o] = at(bkg, o); + if (o >= static_cast(outcomes.size()) || outcomes[o].reflections.empty()) + continue; + const IntegrationOutcome &out = outcomes[o]; + + // The rotation R taking the sweep's lattice, turned to this frame as the per-image path + // turns it, onto the frame's refined lattice: R = A P^-1, with P^-1's rows the reciprocal + // vectors of P. Its antisymmetric part is the small residual rotation. + const float angle = gon.GetAngle_deg(static_cast(o)) + gon.GetWedge_deg() / 2.0f; + const CrystalLattice pred = rot.lattice.Multiply(gon.GetTransformationAngle(-angle)); + const Coord pv[3] = {pred.Vec0(), pred.Vec1(), pred.Vec2()}; + Coord av[3] = {out.latt.Vec0(), out.latt.Vec1(), out.latt.Vec2()}; + for (int i = 0; i < 3; i++) + if (av[i] * pv[i] < 0.0f) // a basis sign convention, not a rotation + av[i] = -av[i]; + const float det = pv[0] * (pv[1] % pv[2]); + if (std::abs(det) > 0.0f) { + const Coord rows[3] = {(pv[1] % pv[2]) / det, (pv[2] % pv[0]) / det, (pv[0] % pv[1]) / det}; + double R[3][3]; + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) + R[j][k] = av[0][j] * rows[0][k] + av[1][j] * rows[1][k] + av[2][j] * rows[2][k]; + const Coord w(static_cast(0.5 * (R[2][1] - R[1][2])), + static_cast(0.5 * (R[0][2] - R[2][0])), + static_cast(0.5 * (R[1][0] - R[0][1]))); + const double about = w * spindle; + frame[ABOUT_SPINDLE][o] = about * 180.0 / M_PI; + frame[OFF_SPINDLE][o] = (w - spindle * static_cast(about)).Length() * 180.0 / M_PI; + } + + // The cell the strong reflections ask for: d at the observed centroid against d predicted. + double sum = 0.0; + int cnt = 0; + for (const auto &r : out.reflections) { + if (!r.observed || r.on_ice_ring || !(r.sigma > 0.0f) || r.I < 3.0f * r.sigma || !(r.d > 0.0f)) + continue; + const float q = out.geom.DetectorToRecip(r.observed_x, r.observed_y).Length(); + if (!(q > 0.0f)) + continue; + sum += 1.0 / (q * r.d) - 1.0; + cnt++; + } + if (cnt > 0) + frame[CELL][o] = sum / cnt; + } + + // The blocks wholly inside the first part, and those wholly past it; a block that straddles + // the boundary is neither, and a last block shorter than half a block is left out. + struct Range { int first, last; }; + std::vector early_blocks, late_blocks; + for (int first = 0; first < n; first += block) { + const int last = std::min(n, first + block); + if (last - first < (block + 1) / 2) + continue; + if (last <= early) + early_blocks.push_back({first, last}); + else if (first >= early) + late_blocks.push_back({first, last}); + } + const auto deg = [&](int ordinal) { return gon.GetAngle_deg(static_cast(ordinal)); }; + + // Each quantity's level over the rest of the sweep against its level on the first part: the + // block means of each part are the samples, so what one part's blocks scatter by - orientation, + // dose, counting noise alike - is that part's noise, and a single block out of line is not a + // shift of the part (Welch's t, its degrees of freedom rounded down). + struct Comparison { int stat; double late, early, noise, t, p; }; + std::vector comparisons; + for (int s = 0; s < N_STATS; s++) { + const auto sample = [&](const std::vector &ranges) { + std::vector means; + for (const auto &r : ranges) + if (const Block b = Summarise(frame[s], r.first, r.last); std::isfinite(b.mean)) + means.push_back(b.mean); + double m = 0.0, v = 0.0; + for (const double x : means) + m += x; + m /= std::max(1, means.size()); + for (const double x : means) + v += (x - m) * (x - m); + v /= std::max(2, means.size()) - 1; + return std::tuple{m, v, static_cast(means.size())}; + }; + const auto [me, ve, ne] = sample(early_blocks); + const auto [ml, vl, nl] = sample(late_blocks); + if (ne < 3 || nl < 3) + continue; + const double ae = ve / ne, al = vl / nl; + const double noise = std::sqrt(ae + al); + const double t = noise > 0.0 ? (ml - me) / noise : (ml == me ? 0.0 : INFINITY); + const double dof = noise > 0.0 ? (ae + al) * (ae + al) / (ae * ae / (ne - 1) + al * al / (nl - 1)) + : ne + nl - 2; + comparisons.push_back({s, ml, me, noise, t, TwoSidedP(t, std::max(1, static_cast(dof)))}); + } + const double alpha = std::erfc(SPOT_BUDGET_SIGNIFICANCE_Z / std::sqrt(2.0)); + const double family = static_cast(comparisons.size()); + const auto describe = [&](const Comparison &c) { + const Stat &st = stats[c.stat]; + return fmt::format("{} is {:.3g}{} over the rest of the sweep against {:.3g}{} on the first part " + "(+- {:.2g}, t = {:.1f}, p = {:.2g} in a family of {})", + st.what, c.late * st.unit, st.suffix, c.early * st.unit, st.suffix, + c.noise * st.unit, c.t, c.p, comparisons.size()); + }; + std::vector revisits; + std::string by_quantity; // each quantity's p times the family, for the record + const Comparison *closest = nullptr; + for (const auto &c : comparisons) { + by_quantity += fmt::format("{}{} {:.2g}", by_quantity.empty() ? "" : ", ", stats[c.stat].tag, + std::min(1.0, c.p * family)); + if (c.p * family < alpha) + revisits.push_back(fmt::format("{} would be revisited because {}", stats[c.stat].decision, + describe(c))); + if (!closest || c.p < closest->p) + closest = &c; + } + + const std::string blocks = fmt::format("{} early and {} late blocks of {:.0f} deg", early_blocks.size(), + late_blocks.size(), block * step); + if (!closest) + logger.Info("far-end check: nothing to compare ({}; the first part ends at {:.0f} deg)", blocks, + deg(early)); + else if (revisits.empty()) + logger.Info("far-end check: early decisions hold ({}; closest: {}; p x family by quantity: {})", + blocks, describe(*closest), by_quantity); + else { + std::string line; + for (const auto &r : revisits) + line += (line.empty() ? "" : "; ") + r; + logger.Info("far-end check: {} ({}; p x family by quantity: {})", line, blocks, by_quantity); + } +} diff --git a/rugnux/RugnuxFirstPass.cpp b/rugnux/RugnuxFirstPass.cpp index cc9f06b22..d2d8e1b5e 100644 --- a/rugnux/RugnuxFirstPass.cpp +++ b/rugnux/RugnuxFirstPass.cpp @@ -393,6 +393,7 @@ bool Rugnux::FirstPassRotationIndexing(PipelineLocals &p) { // not pay for it - measured, an unguarded call spent 41-76 % of the wall clock on a // healthy dataset evaluating rungs that arithmetically could not be adopted (one crystal // 41 s guarded against 175 s unguarded, same answer). + run.NearTieGate("ladder floor", "standing pass", best.score, n_validation / 6); if (best.score < n_validation / 6) { if (auto rescued = run.FirstPassLadder(best.score, n_validation / 6, 0)) best = *rescued; @@ -840,6 +841,8 @@ auto Rugnux::FirstPassRun::ScoreSchemes(const SchemeIndexers &ris, IndexAndRefin if (result.unconstrained && score < majority) { RotationIndexerResult alt = *result.unconstrained; const int alt_score = CountIndexed(idx, alt); + NearTieGate("metric-symmetry demotion", fmt::format("scheme '{}' unconstrained", name), + alt_score, majority); if (alt_score > majority) { logger.Info("Scheme '{}': {}-centred {} indexes {}/{} frames but its unconstrained " "cell indexes {}/{} - the metric symmetry is a false promotion, dropping it", @@ -871,6 +874,9 @@ auto Rugnux::FirstPassRun::ScoreSchemes(const SchemeIndexers &ris, IndexAndRefin // A later scheme wins if it indexes clearly more frames (>10%). const bool clearly_more = static_cast(score) > bp.score * 1.1f + 0.5f; + if (bp.result.has_value() && !SameLattice(result, *bp.result)) + NearTieFrames("scheme choice", fmt::format("'{}'", name), score, + fmt::format("'{}'", bp.name), bp.score); // Integer-supercell tie-break. The validation-frame count saturates - a spurious axis // multiple (2x/3x...) indexes every frame its true cell does, so both schemes reach the @@ -923,6 +929,10 @@ auto Rugnux::FirstPassRun::ScoreSchemes(const SchemeIndexers &ris, IndexAndRefin const CrystalLattice sub = SubLattice(larger_p, ev, static_cast(nearest)).NiggliReduce(); const size_t n_sub = IndexedMillerIndices(sub, cloud, index_tol).size(); const bool larger_is_real = hkl.size() > n_smaller && hkl_p.size() > n_sub; + NearTieCounts("axis-harmonic arbiter", "larger cell", static_cast(hkl.size()), + "smaller cell", static_cast(n_smaller)); + NearTieCounts("axis-harmonic arbiter (own sub-lattice)", "larger cell", + static_cast(hkl_p.size()), "sub-lattice", static_cast(n_sub)); logger.Info("Axis-harmonic arbiter: '{}' ({:.0f} A^3) vs '{}' ({:.0f} A^3), {:.2f}x - " "the larger cell accounts for {} validation spots against the smaller " "cell's {}; on its own geometry {} against its ({},{},{}) sub-lattice's " @@ -1188,6 +1198,39 @@ auto Rugnux::FirstPassRun::FirstPassLadder(int score_before, int min_gain, int m won = static_cast(i); } + // Near ties (log only): the rung that indexes the most against the bars it has to clear, and + // against the best rung that found a different lattice. + { + int top = -1; + for (size_t i = 0; i < rungs.size(); i++) + if (rung_score[i] >= 0 && (top < 0 || rung_score[i] > rung_score[top])) + top = static_cast(i); + if (top >= 0) { + const double n = static_cast(validation.size()); + const auto var = [n](int k) { + const double p = (k + 0.5) / (n + 1.0); + return n * p * (1.0 - p); + }; + NearTie("ladder adoption", fmt::format("best rung gains {}/{}", rung_score[top] - score_before, + validation.size()), + fmt::format("min gain {}", min_gain), rung_score[top] - score_before - min_gain, + std::sqrt(var(rung_score[top]) + var(std::max(0, score_before))), + "binomial on validation frames"); + if (min_score > 0) + NearTieGate("ladder adoption (bar)", "best rung", rung_score[top], min_score); + } + if (won >= 0) { + int rival = -1; + for (size_t i = 0; i < rungs.size(); i++) + if (static_cast(i) != won && rung_score[i] >= 0 && rung_score[i] - score_before >= min_gain + && rung_score[i] >= min_score && !SameLattice(*rung_pass[i].result, *rung_pass[won].result) + && (rival < 0 || rung_score[i] > rung_score[rival])) + rival = static_cast(i); + if (rival >= 0) + NearTieFrames("ladder rung choice", "winner", rung_score[won], "other lattice", rung_score[rival]); + } + } + // ...and it must not be an axis SUB-MULTIPLE of a cell another rung found. The frame // count the choice above is made on cannot arbitrate an axis harmonic - a cell twice as // long has to place every spot twice as accurately to score the same, which the scheme @@ -1243,6 +1286,8 @@ auto Rugnux::FirstPassRun::FirstPassLadder(int score_before, int min_gain, int m won_vol, rung_score[i], static_cast(validation.size()), rung_pass[i].vol, ratio, hkl.size(), won_spots, 100.0 * ev.occupancy, ev.u, ev.v, ev.w); + NearTieCounts("ladder axis-harmonic arbiter", "larger rung", static_cast(hkl.size()), + "winner", static_cast(won_spots)); if (hkl.size() > won_spots) { logger.Warning("First-pass ladder: the rung that indexes the most frames does " "so on a {:.0f}x sub-multiple of an axis another rung found, " @@ -1524,3 +1569,59 @@ auto Rugnux::FirstPassRun::TwinCompositeHypothesis(FirstPass best) -> FirstPass } return best; } + +// Which first pass a near-tie line comes from: the geometry pre-pass, the canonical pass, or a probe +// that only scores a lattice. +std::string Rugnux::FirstPassRun::PassTag() const { + return geometry_prepass ? "pass 1" : rx.indexing_probe_only_ ? "probe" : "pass 2"; +} + +void Rugnux::FirstPassRun::NearTie(const std::string &decision, const std::string &a, const std::string &b, + double margin, double noise, const char *model) const { + if (std::abs(margin) <= NEAR_TIE_SIGMA * noise) + logger.Info("near tie: {} {} vs {} margin {:.3g} (noise {:.3g}, {}; {}) - would carry both", + decision, a, b, margin, noise, model, PassTag()); +} + +void Rugnux::FirstPassRun::NearTieFrames(const std::string &decision, const std::string &a, int ka, + const std::string &b, int kb) const { + // Two lattices that both index fewer frames than the bar the run keeps a lattice on (a sixth of + // the validation frames, see FirstPassRotationIndexing) are not two hypotheses about the crystal. + if (std::max(ka, kb) < static_cast(validation.size()) / 6) + return; + const double n = static_cast(validation.size()); + const auto var = [n](int k) { + const double p = (std::max(k, 0) + 0.5) / (n + 1.0); + return n * p * (1.0 - p); + }; + NearTie(decision, fmt::format("{} {}/{}", a, ka, validation.size()), + fmt::format("{} {}/{}", b, kb, validation.size()), + ka - kb, std::sqrt(var(ka) + var(kb)), "binomial on validation frames"); +} + +void Rugnux::FirstPassRun::NearTieGate(const std::string &decision, const std::string &what, int score, + double gate) const { + const double n = static_cast(validation.size()); + const double p = (std::max(score, 0) + 0.5) / (n + 1.0); + NearTie(decision, fmt::format("{} {}/{}", what, score, validation.size()), + fmt::format("gate {:.1f}", gate), score - gate, std::sqrt(n * p * (1.0 - p)), + "binomial on validation frames"); +} + +void Rugnux::FirstPassRun::NearTieCounts(const std::string &decision, const std::string &a, int64_t na, + const std::string &b, int64_t nb) const { + NearTie(decision, fmt::format("{} {}", a, na), fmt::format("{} {}", b, nb), + static_cast(na - nb), std::sqrt(static_cast(na + nb)), "Poisson"); +} + +void Rugnux::FirstPassRun::NearTiePooled(const std::string &decision, const std::string &a, + const PooledEvidence &ea, const std::string &b, + const PooledEvidence &eb) const { + const auto excess = [](const PooledEvidence &e) { + return 100.0 * static_cast(e.on_lattice - e.by_chance) + / static_cast(std::max(1, e.spots)); + }; + const auto [difference, sigma] = ValidationEvidenceDifference(ea, eb); + NearTie(decision, fmt::format("{} {:.1f}%", a, excess(ea)), fmt::format("{} {:.1f}%", b, excess(eb)), + -difference, sigma, "binomial on pooled validation spots"); +} diff --git a/rugnux/RugnuxFirstPass.h b/rugnux/RugnuxFirstPass.h index f44732c72..c6fdb2d88 100644 --- a/rugnux/RugnuxFirstPass.h +++ b/rugnux/RugnuxFirstPass.h @@ -13,6 +13,9 @@ // RugnuxFirstPassRescues.cpp: the rescues and hypotheses that may replace the standing lattice. // RugnuxFirstPassAccept.cpp: the no-crystal test and what the run takes from the accepted lattice. +#include +#include +#include #include #include #include @@ -32,6 +35,22 @@ namespace rugnux_internal { // the scheme comparison established it on the spots. struct FirstPass { std::optional result; int score = -1; std::string name; double vol = 0.0; std::optional indexer_error; bool larger_settled = false; }; + + // Two results are one lattice when their Niggli-reduced primitive edges agree within 2 % (the + // battery's lattice-identity test), whatever setting or class each is held in. + inline bool SameLattice(const RotationIndexerResult &a, const RotationIndexerResult &b) { + const auto edges = [](const RotationIndexerResult &r) { + const auto uc = r.lattice.ToPrimitive(r.search_result.centering).NiggliReduce().GetUnitCell(); + std::array e = {uc.a, uc.b, uc.c}; + std::sort(e.begin(), e.end()); + return e; + }; + const auto ea = edges(a), eb = edges(b); + for (int i = 0; i < 3; i++) + if (std::abs(eb[i] - ea[i]) > 0.02 * ea[i]) + return false; + return true; + } } struct Rugnux::FirstPassRun { @@ -123,6 +142,29 @@ struct Rugnux::FirstPassRun { void StartCentreSpeculation(); void StartShortAxisSpeculation(); + // --- Near ties (log only; nothing here changes a decision) --- + // + // A decision whose margin lies within NEAR_TIE_SIGMA standard deviations of its own noise was + // taken on noise: the line names the alternative that would have to be carried forward beside + // it. The noise is the comparison's own - binomial for a count of validation frames (on + // n = validation.size(), p = (k + 1/2) / (n + 1)), Poisson for counts of spots or reflections, + // and for two pooled spot shares the binomial standard error ValidationEvidencePrefers uses. + static constexpr double NEAR_TIE_SIGMA = 2.0; + std::string PassTag() const; + void NearTie(const std::string &decision, const std::string &a, const std::string &b, + double margin, double noise, const char *model) const; + // Two frame counts, a against b; nothing where neither reaches the run's floor of a sixth. + void NearTieFrames(const std::string &decision, const std::string &a, int ka, + const std::string &b, int kb) const; + // A frame count against the gate it is compared with. + void NearTieGate(const std::string &decision, const std::string &what, int score, double gate) const; + // Two counts of spots or reflections. + void NearTieCounts(const std::string &decision, const std::string &a, int64_t na, + const std::string &b, int64_t nb) const; + // Two lattices' pooled validation evidence, a the standing one. + void NearTiePooled(const std::string &decision, const std::string &a, const PooledEvidence &ea, + const std::string &b, const PooledEvidence &eb) const; + // --- State --- Rugnux ℞ diff --git a/rugnux/RugnuxFirstPassAccept.cpp b/rugnux/RugnuxFirstPassAccept.cpp index 02bc68706..a78b2c645 100644 --- a/rugnux/RugnuxFirstPassAccept.cpp +++ b/rugnux/RugnuxFirstPassAccept.cpp @@ -147,6 +147,17 @@ 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. + { + 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"); + } if (!ValidationEvidenceBeatsChance(evidence)) { { // Name the cell and Bravais class that was rejected. The commonest cause is a metric diff --git a/rugnux/RugnuxFirstPassRescues.cpp b/rugnux/RugnuxFirstPassRescues.cpp index b6f2872ce..4863f6116 100644 --- a/rugnux/RugnuxFirstPassRescues.cpp +++ b/rugnux/RugnuxFirstPassRescues.cpp @@ -115,12 +115,15 @@ namespace { // // Only after a poor pass, so a correctly-signed file costs nothing; and spot finding is not auto Rugnux::FirstPassRun::AxisSignRescue(FirstPass best) -> FirstPass { + NearTieGate("axis-sign rescue gate", "standing pass", best.score, 0.5 * static_cast(validation.size())); if (!rx.cancelled_ && best.score < 0.5 * static_cast(validation.size())) { if (const auto gon = rx.experiment_.GetGoniometer(); gon.has_value() && gon->IsScanning()) { GoniometerAxis flipped = *gon; flipped.Axis(-gon->GetAxis()); rx.experiment_.Goniometer(flipped); FirstPass alt = PickBest(*indexer_pool, *indexer); + if (alt.result.has_value()) + NearTieFrames("axis sign", "opposite sign", alt.score, "file's sign", best.score); // More frames is not enough on its own while both counts sit at the floor: 1/60 against // 0/60 is noise, and adopting it puts every rescue below on the wrong sign (measured: a // run whose file sign, kept, went on to index found nothing once flipped on 1/60). So @@ -216,6 +219,8 @@ auto Rugnux::FirstPassRun::LongAxisRescue(const FirstPass &standing) -> std::opt constrained = PickBest(*indexer_pool, *indexer); } rx.experiment_.SetUnitCell(std::nullopt); // leave the space-group / cell determination de-novo + if (constrained.result.has_value()) + NearTieFrames("long-axis rescue", "re-index", constrained.score, "standing", standing.score); if (constrained.result.has_value() && constrained.score > standing.score) return constrained; } @@ -283,6 +288,8 @@ auto Rugnux::FirstPassRun::BeamCentreCheck(FirstPass best, bool &long_axis_asked FirstPass alt = IndexAtMeasuredCentre(); const bool header_indexes = best.result.has_value() && best.score > majority; const bool measured_indexes = alt.result.has_value() && alt.score > majority; + NearTieGate("beam-centre majority gate", "file centre", best.score, majority); + NearTieGate("beam-centre majority gate", "measured centre", alt.score, majority); if (!header_indexes) return CentreWhenFileFails(std::move(best), std::move(alt), c, long_axis_asked); RestoreBeamCenter(header_x, header_y); @@ -408,6 +415,8 @@ auto Rugnux::FirstPassRun::CentreWhenFileFails(FirstPass best, FirstPass alt, co && home->result->search_result.centering == away->result->search_result.centering && std::max(home->vol, away->vol) < 1.02 * std::min(home->vol, away->vol); + if (home && away && !same_lattice) + NearTieFrames("beam centre (lean ladder)", "file centre", home->score, "measured centre", away->score); if (home && (!away || same_lattice || home->score >= away->score)) { committed = home_state; RestoreBeamCenter(header_x, header_y); @@ -498,6 +507,8 @@ auto Rugnux::FirstPassRun::CentreByPooledSpots(FirstPass best, FirstPass alt, co if (best.result.has_value()) at_header = MeasurePooled(*indexer, *best.result); long_axis_asked = true; + if (alt.result.has_value() && best.result.has_value()) + NearTiePooled("beam centre (pooled spots)", "file centre", at_header, "measured centre", at_measured); if (alt.result.has_value() && ValidationEvidenceBeatsChance(at_measured) && ValidationEvidencePrefers(at_header, at_measured)) { TryBeamCenter(measured_x, measured_y); @@ -674,6 +685,7 @@ auto Rugnux::FirstPassRun::CentresDisagree(FirstPass best, const FirstPass &alt, const PooledEvidence at_measured = MeasurePooled(*indexer, *alt.result); adopted_measured = ValidationEvidenceBeatsChance(at_measured) && ValidationEvidencePrefers(at_header, at_measured); + NearTiePooled("beam centre (harmonic pair)", "file centre", at_header, "measured centre", at_measured); logger.Info("Beam centre check: pooled validation spots on the lattice " "{:.1f}% at the file's centre and {:.1f}% at the measured " "(chance {:.1f}% and {:.1f}%)", @@ -713,6 +725,7 @@ auto Rugnux::FirstPassRun::CentresDisagree(FirstPass best, const FirstPass &alt, // counted in the pass a rescue is gated on, it lifts that pass past the gate and the rescue // that would have indexed every frame is never asked. auto Rugnux::FirstPassRun::AddEndWedge(FirstPass bp, IndexerThreadPool &pool, IndexAndRefine &idx) -> FirstPass { + NearTieGate("end-wedge gate", "standing pass", bp.score, static_cast(validation.size()) / 2); if (rx.cancelled_ || bp.score > static_cast(validation.size()) / 2) return bp; SchemeIndexers ris; @@ -840,6 +853,7 @@ auto Rugnux::FirstPassRun::ShortAxisHypothesis(FirstPass best) -> FirstPass { } const int alt_score = CountIndexed(*indexer, alt); + NearTieFrames("short-axis hypothesis", "sub-cell", alt_score, "standing", best.score); // A pass that indexes nothing is not tied with a pass that indexes nothing. Without // this, 0 against 0 passes the tie band and the volume ratio alone decides, which is // the one thing this block must never do. @@ -936,11 +950,15 @@ auto Rugnux::FirstPassRun::SubLatticeHypothesis(FirstPass best) -> FirstPass { } bool adopt = false; std::string verdict = " - none indexes a majority of the frames and as many as the standing cell"; + if (sub_best) + NearTieGate("sub-lattice majority gate", "best sub-lattice", sub_best->score, + static_cast(validation.size()) / 2); if (sub_best && sub_best->score >= best.score && sub_best->score > static_cast(validation.size()) / 2) { const PooledEvidence at_standing = MeasurePooled(*indexer, *best.result); const PooledEvidence at_sub = MeasurePooled(*indexer, *sub_best->result); adopt = ValidationEvidencePrefers(at_standing, at_sub); + NearTiePooled("sub-lattice hypothesis", "standing cell", at_standing, "sub-lattice", at_sub); verdict = fmt::format(" - the best puts {}/{} pooled spots on itself ({} by chance) against " "the standing cell's {}/{} ({}), {}", at_sub.on_lattice, at_sub.spots, at_sub.by_chance, @@ -1017,6 +1035,9 @@ auto Rugnux::FirstPassRun::Pass1Fallback(FirstPass best) -> FirstPass { const int floor = static_cast(validation.size()) / 6; const int pass1_score = (best.score < floor || !same_lattice) ? CountIndexed(*indexer, pass1_lattice) : -1; + NearTieGate("pass-1 fallback floor", "re-index", best.score, floor); + if (pass1_score >= 0 && !same_lattice) + NearTieFrames("pass-1 fallback", "re-index", best.score, "pass-1 lattice", pass1_score); if (best.score < floor || pass1_score > best.score) { const auto &pc = pass1_lattice.search_result.conventional.GetUnitCell(); const auto &bc = best.result->search_result.conventional.GetUnitCell(); diff --git a/rugnux/RugnuxImagePass.cpp b/rugnux/RugnuxImagePass.cpp index 0f7e3e64a..f0d159506 100644 --- a/rugnux/RugnuxImagePass.cpp +++ b/rugnux/RugnuxImagePass.cpp @@ -584,6 +584,9 @@ void Rugnux::ImagePass(PipelineLocals &p) { // The consensus cell below is this lattice's conventional cell, so this is the centring // it is written in - the pair the volume comparisons in Run() need. result.consensus_centering = rot->search_result.centering; + // Log only: what the rest of the sweep says about the decisions taken on its first part. + if (!p.geometry_prepass && p.write_files && experiment_.IsRotationIndexing()) + FarEndCheck(p, *rot); } result.consensus_cell = indexer->GetConsensusUnitCell(); end_msg.unit_cell = result.consensus_cell;