From 4700dd57339fcdf659be26ab611296b8acca01ea Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 9 Oct 2026 21:04:51 +0200 Subject: [PATCH 1/4] rugnux: log-only far-end check and near-tie lines on the first-pass decisions Part B takes the decisions before pass 2 from the first part of a rotation sweep. Two diagnostics now say, in the log only, where that may matter; nothing the run writes changes. Far-end check (RugnuxFarEnd.cpp), after the image pass of the canonical pass: the sweep in 10 deg blocks, six per-frame quantities pass 2 already holds - frames indexing on their own, share of spots on the lattice, the cell scale the strong reflections ask for (d at the observed centroid over d predicted), the residual rotation of each frame's refined orientation against the sweep's lattice turned to that frame (about the spindle and off it), and the background - each late block against the early blocks as a Student-t prediction from their own spread, family-wise at the run's 0.1 % significance over every late-block comparison. One line: "far-end check: early decisions hold" or "far-end check: would be revisited because ". Near ties (FirstPassRun::NearTie*): where a first-pass decision compares two candidates or a count against a gate - scheme choice, metric-symmetry demotion, axis-harmonic arbiters, axis sign, long-axis rescue, beam-centre gates and choices, end-wedge gate, short-axis and sub-lattice hypotheses, ladder floor/adoption/rung choice, pass-1 fallback, no-crystal test - a margin within two standard deviations of the comparison's own noise (binomial on validation frames, Poisson on spot counts, the binomial standard error ValidationEvidencePrefers uses on pooled spots) is logged as "near tie: vs margin x (noise y) - would carry both". ValidationEvidencePrefers now reads its difference and standard error from ValidationEvidenceDifference, which the near-tie lines share; same arithmetic. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi --- rugnux/CMakeLists.txt | 1 + rugnux/Rugnux.cpp | 11 +- rugnux/Rugnux.h | 7 + rugnux/RugnuxFarEnd.cpp | 265 ++++++++++++++++++++++++++++++ rugnux/RugnuxFirstPass.cpp | 97 +++++++++++ rugnux/RugnuxFirstPass.h | 42 +++++ rugnux/RugnuxFirstPassAccept.cpp | 11 ++ rugnux/RugnuxFirstPassRescues.cpp | 21 +++ rugnux/RugnuxImagePass.cpp | 3 + 9 files changed, 455 insertions(+), 3 deletions(-) create mode 100644 rugnux/RugnuxFarEnd.cpp 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 3c482e382..4ee51216c 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 1147a40f1..ce17d9d62 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -233,6 +233,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 @@ -923,6 +927,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..720c616f9 --- /dev/null +++ b/rugnux/RugnuxFarEnd.cpp @@ -0,0 +1,265 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#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 every block past the first part is compared with the blocks inside +// it. A late block is measured against the early blocks' mean in units of the spread of the early +// blocks themselves, as a prediction for one more block (Student's t on the early blocks' degrees +// of freedom). A block is out of line when its two-sided p, times the number of late-block +// comparisons the run made over all the quantities, 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, however many blocks and quantities there are. 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 Stat stats[N_STATS] = { + {"frames indexing on their own", "the lattice", 100.0, "%"}, + {"spots on the lattice", "the lattice", 100.0, "%"}, + {"cell scale of the strong reflections (d observed / d predicted - 1)", "the cell", 100.0, "%"}, + {"orientation residual about the spindle", "the rotation scale", 1.0, " deg"}, + {"orientation residual off the spindle", "the orientation", 1.0, " deg"}, + {"background", "the shadow mask", 1.0, ""}, + }; + 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)); }; + + // Every late block of every quantity against the early blocks of the same quantity. + struct Comparison { int stat; Range range; double value, mean, noise, t, p; }; + std::vector comparisons; + for (int s = 0; s < N_STATS; s++) { + std::vector means; + for (const auto &r : early_blocks) + if (const Block b = Summarise(frame[s], r.first, r.last); std::isfinite(b.mean)) + means.push_back(b.mean); + const int ne = static_cast(means.size()); + if (ne < 3) + continue; + double me = 0.0; + for (const double m : means) + me += m; + me /= ne; + double ve = 0.0; + for (const double m : means) + ve += (m - me) * (m - me); + ve /= ne - 1; + for (const auto &r : late_blocks) { + const Block b = Summarise(frame[s], r.first, r.last); + if (!std::isfinite(b.mean)) + continue; + const double noise = std::sqrt(ve * (1.0 + 1.0 / ne) + b.se * b.se); + const double t = noise > 0.0 ? (b.mean - me) / noise : (b.mean == me ? 0.0 : INFINITY); + comparisons.push_back({s, r, b.mean, me, noise, t, TwoSidedP(t, ne - 1)}); + } + } + 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, int out_of_line, int late) { + const Stat &st = stats[c.stat]; + return fmt::format("{} over {:.0f}-{:.0f} deg is {:.3g}{} against {:.3g}{} +- {:.2g} on the first part " + "(t = {:.1f}, p = {:.2g} in a family of {}; {} of {} late blocks out of line)", + st.what, deg(c.range.first), deg(c.range.last), c.value * st.unit, st.suffix, + c.mean * st.unit, st.suffix, c.noise * st.unit, c.t, c.p, comparisons.size(), + out_of_line, late); + }; + std::vector revisits; + const Comparison *closest = nullptr; + for (int s = 0; s < N_STATS; s++) { + const Comparison *worst = nullptr; + int out_of_line = 0, late = 0; + for (const auto &c : comparisons) + if (c.stat == s) { + late++; + if (c.p * family < alpha) + out_of_line++; + if (!worst || c.p < worst->p) + worst = &c; + } + if (!worst) + continue; + if (out_of_line > 0) + revisits.push_back(fmt::format("{} would be revisited because {}", stats[s].decision, + describe(*worst, out_of_line, late))); + if (!closest || worst->p < closest->p) + closest = worst; + } + const auto late_of = [&](int s) { + return static_cast(std::count_if(comparisons.begin(), comparisons.end(), + [s](const Comparison &c) { return c.stat == s; })); + }; + + 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: {})", blocks, + describe(*closest, 0, late_of(closest->stat))); + else { + std::string line; + for (const auto &r : revisits) + line += (line.empty() ? "" : "; ") + r; + logger.Info("far-end check: {} ({})", line, blocks); + } +} diff --git a/rugnux/RugnuxFirstPass.cpp b/rugnux/RugnuxFirstPass.cpp index dade24cc5..580861825 100644 --- a/rugnux/RugnuxFirstPass.cpp +++ b/rugnux/RugnuxFirstPass.cpp @@ -392,6 +392,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; @@ -839,6 +840,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", @@ -870,6 +873,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 @@ -922,6 +928,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 " @@ -1184,6 +1194,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 @@ -1239,6 +1282,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, " @@ -1411,3 +1456,55 @@ auto Rugnux::FirstPassRun::IndexSchemes() -> FirstPass { StartShortAxisSpeculation(); }); } + +// 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 { + 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 a268c97eb..87bd1b85b 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 @@ -30,6 +33,22 @@ namespace rugnux_internal { // verdict on the crystal; an empty result with one is a verdict on the machine. struct FirstPass { std::optional result; int score = -1; std::string name; double vol = 0.0; std::optional indexer_error; }; + + // 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 { @@ -120,6 +139,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. + 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 abca5c1bf..60d826f7a 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; From a2babf7e3d68ab0dede315e8da6dd44a47aca5c4 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 9 Oct 2026 22:58:32 +0200 Subject: [PATCH 2/4] rugnux far-end check: both directions, the rotation rate, every quantity's p in the line The first part can be the odd one (a crystal that settles or drifts over the opening degrees), so every early block is also tested against the late ones. The rotation scale is also checked as a rate: the slope of the residual about the spindle over the late blocks against that over the early ones. Every far-end line now ends with each quantity's smallest p times the family size, so the population can be read off the logs. Still log only. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi --- rugnux/RugnuxFarEnd.cpp | 144 ++++++++++++++++++++++++++++------------ 1 file changed, 103 insertions(+), 41 deletions(-) diff --git a/rugnux/RugnuxFarEnd.cpp b/rugnux/RugnuxFarEnd.cpp index 720c616f9..f379d5539 100644 --- a/rugnux/RugnuxFarEnd.cpp +++ b/rugnux/RugnuxFarEnd.cpp @@ -18,9 +18,9 @@ // 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 every block past the first part is compared with the blocks inside -// it. A late block is measured against the early blocks' mean in units of the spread of the early -// blocks themselves, as a prediction for one more block (Student's t on the early blocks' degrees -// of freedom). A block is out of line when its two-sided p, times the number of late-block +// it, and every block inside it with those past it. A block is measured against the other part's +// mean in units of the spread of that part's blocks themselves, as a prediction for one more block +// (Student's t on that part's degrees of freedom). A block is out of line when its two-sided p, times the number of late-block // comparisons the run made over all the quantities, 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, however many blocks and quantities there are. No number in it is set @@ -32,7 +32,9 @@ // - 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; +// to that frame, about the spindle and off it: the rotation scale and the orientation - and for +// the rotation scale also its rate, the slope of the residual about the spindle against the +// angle, late blocks against early ones; // - the background: the shadow mask (hardware that shadows more, or less, as the sample turns). namespace { @@ -106,15 +108,17 @@ void Rugnux::FarEndCheck(PipelineLocals &p, const RotationIndexerResult &rot) { 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 Stat stats[N_STATS] = { - {"frames indexing on their own", "the lattice", 100.0, "%"}, - {"spots on the lattice", "the lattice", 100.0, "%"}, - {"cell scale of the strong reflections (d observed / d predicted - 1)", "the cell", 100.0, "%"}, - {"orientation residual about the spindle", "the rotation scale", 1.0, " deg"}, - {"orientation residual off the spindle", "the orientation", 1.0, " deg"}, - {"background", "the shadow mask", 1.0, ""}, + enum { INDEXED, ON_LATTICE, CELL, ABOUT_SPINDLE, OFF_SPINDLE, BACKGROUND, N_STATS, RATE = N_STATS }; + struct Stat { const char *what; const char *decision; double unit; const char *suffix; const char *tag; }; + const Stat stats[N_STATS + 1] = { + {"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"}, + {"rate of the residual about the spindle (per 100 deg of rotation)", "the rotation scale", 100.0, "%", + "rate"}, }; std::vector> frame(N_STATS, std::vector(n, NAN)); const Coord spindle = gon.GetAxis().Normalize(); @@ -184,52 +188,108 @@ void Rugnux::FarEndCheck(PipelineLocals &p, const RotationIndexerResult &rot) { } const auto deg = [&](int ordinal) { return gon.GetAngle_deg(static_cast(ordinal)); }; - // Every late block of every quantity against the early blocks of the same quantity. - struct Comparison { int stat; Range range; double value, mean, noise, t, p; }; + // Every late block of every quantity against the early blocks of the same quantity - and every + // early block against the late ones, because the part that is out of line can be the first: a + // crystal that settles, or an orientation that drifts, over the opening degrees makes the first + // part the odd one, and its spread then hides nothing the late blocks show. + struct Comparison { int stat; Range range; double value, mean, noise, t, p; bool first_part; }; std::vector comparisons; - for (int s = 0; s < N_STATS; s++) { + const auto compare = [&](int s, const std::vector &reference, const std::vector &tested, + bool first_part) { std::vector means; - for (const auto &r : early_blocks) + for (const auto &r : reference) if (const Block b = Summarise(frame[s], r.first, r.last); std::isfinite(b.mean)) means.push_back(b.mean); - const int ne = static_cast(means.size()); - if (ne < 3) - continue; - double me = 0.0; + const int nr = static_cast(means.size()); + if (nr < 3) + return; + double mr = 0.0; for (const double m : means) - me += m; - me /= ne; - double ve = 0.0; + mr += m; + mr /= nr; + double vr = 0.0; for (const double m : means) - ve += (m - me) * (m - me); - ve /= ne - 1; - for (const auto &r : late_blocks) { + vr += (m - mr) * (m - mr); + vr /= nr - 1; + for (const auto &r : tested) { const Block b = Summarise(frame[s], r.first, r.last); if (!std::isfinite(b.mean)) continue; - const double noise = std::sqrt(ve * (1.0 + 1.0 / ne) + b.se * b.se); - const double t = noise > 0.0 ? (b.mean - me) / noise : (b.mean == me ? 0.0 : INFINITY); - comparisons.push_back({s, r, b.mean, me, noise, t, TwoSidedP(t, ne - 1)}); + const double noise = std::sqrt(vr * (1.0 + 1.0 / nr) + b.se * b.se); + const double t = noise > 0.0 ? (b.mean - mr) / noise : (b.mean == mr ? 0.0 : INFINITY); + comparisons.push_back({s, r, b.mean, mr, noise, t, TwoSidedP(t, nr - 1), first_part}); + } + }; + for (int s = 0; s < N_STATS; s++) { + compare(s, early_blocks, late_blocks, false); + compare(s, late_blocks, early_blocks, true); + } + // The rotation scale is a RATE: a stage turning at the wrong rate leaves a residual about the + // spindle that grows with the angle, and the scale decided on the first part zeroes its slope + // there. So the slope of the late blocks is compared with that of the early ones, each a + // straight-line fit to its blocks' means with its own standard error. + { + struct Fit { double slope = NAN, se = NAN; int blocks = 0; }; + const auto fit = [&](const std::vector &ranges) { + std::vector x, y; + for (const auto &r : ranges) + if (const Block b = Summarise(frame[ABOUT_SPINDLE], r.first, r.last); std::isfinite(b.mean)) { + x.push_back(0.5 * (deg(r.first) + deg(r.last))); + y.push_back(b.mean); + } + Fit f; + f.blocks = static_cast(x.size()); + if (f.blocks < 3) + return f; + double mx = 0.0, my = 0.0; + for (int i = 0; i < f.blocks; i++) { + mx += x[i]; + my += y[i]; + } + mx /= f.blocks; + my /= f.blocks; + double sxx = 0.0, sxy = 0.0; + for (int i = 0; i < f.blocks; i++) { + sxx += (x[i] - mx) * (x[i] - mx); + sxy += (x[i] - mx) * (y[i] - my); + } + f.slope = sxy / sxx; + double rss = 0.0; + for (int i = 0; i < f.blocks; i++) { + const double e = y[i] - my - f.slope * (x[i] - mx); + rss += e * e; + } + f.se = std::sqrt(rss / (f.blocks - 2) / sxx); + return f; + }; + const Fit fe = fit(early_blocks), fl = fit(late_blocks); + if (fe.blocks >= 3 && fl.blocks >= 3) { + const double noise = std::sqrt(fe.se * fe.se + fl.se * fl.se); + const double t = noise > 0.0 ? (fl.slope - fe.slope) / noise : 0.0; + comparisons.push_back({RATE, {late_blocks.front().first, late_blocks.back().last}, fl.slope, + fe.slope, noise, t, TwoSidedP(t, fe.blocks + fl.blocks - 4), false}); } } 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, int out_of_line, int late) { const Stat &st = stats[c.stat]; - return fmt::format("{} over {:.0f}-{:.0f} deg is {:.3g}{} against {:.3g}{} +- {:.2g} on the first part " - "(t = {:.1f}, p = {:.2g} in a family of {}; {} of {} late blocks out of line)", + return fmt::format("{} over {:.0f}-{:.0f} deg is {:.3g}{} against {:.3g}{} +- {:.2g} on {} " + "(t = {:.1f}, p = {:.2g} in a family of {}; {} of {} blocks out of line)", st.what, deg(c.range.first), deg(c.range.last), c.value * st.unit, st.suffix, - c.mean * st.unit, st.suffix, c.noise * st.unit, c.t, c.p, comparisons.size(), - out_of_line, late); + c.mean * st.unit, st.suffix, c.noise * st.unit, + c.first_part ? "the rest of the sweep" : "the first part", c.t, c.p, + comparisons.size(), out_of_line, late); }; std::vector revisits; + std::string by_quantity; // each quantity's smallest p times the family, for the record const Comparison *closest = nullptr; - for (int s = 0; s < N_STATS; s++) { + for (int s = 0; s <= RATE; s++) { const Comparison *worst = nullptr; - int out_of_line = 0, late = 0; + int out_of_line = 0, tested = 0; for (const auto &c : comparisons) if (c.stat == s) { - late++; + tested++; if (c.p * family < alpha) out_of_line++; if (!worst || c.p < worst->p) @@ -237,9 +297,11 @@ void Rugnux::FarEndCheck(PipelineLocals &p, const RotationIndexerResult &rot) { } if (!worst) continue; + by_quantity += fmt::format("{}{} {:.2g}", by_quantity.empty() ? "" : ", ", stats[s].tag, + std::min(1.0, worst->p * family)); if (out_of_line > 0) revisits.push_back(fmt::format("{} would be revisited because {}", stats[s].decision, - describe(*worst, out_of_line, late))); + describe(*worst, out_of_line, tested))); if (!closest || worst->p < closest->p) closest = worst; } @@ -254,12 +316,12 @@ void Rugnux::FarEndCheck(PipelineLocals &p, const RotationIndexerResult &rot) { 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: {})", blocks, - describe(*closest, 0, late_of(closest->stat))); + logger.Info("far-end check: early decisions hold ({}; closest: {}; p x family by quantity: {})", + blocks, describe(*closest, 0, late_of(closest->stat)), by_quantity); else { std::string line; for (const auto &r : revisits) line += (line.empty() ? "" : "; ") + r; - logger.Info("far-end check: {} ({})", line, blocks); + logger.Info("far-end check: {} ({}; p x family by quantity: {})", line, blocks, by_quantity); } } From 663cb2b276dbe8e99013504428ffa1f655e4bc09 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 9 Oct 2026 23:10:07 +0200 Subject: [PATCH 3/4] rugnux near ties: no frame-count tie between two lattices that both miss the floor Two candidates that both index fewer than a sixth of the validation frames (the bar the run keeps a lattice on) are not two hypotheses about the crystal; on the first 34 battery sets most scheme "ties" were 0/60 against 0/60. Log only, as before. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi --- rugnux/RugnuxFirstPass.cpp | 4 ++++ rugnux/RugnuxFirstPass.h | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/rugnux/RugnuxFirstPass.cpp b/rugnux/RugnuxFirstPass.cpp index 580861825..bc6e1cd84 100644 --- a/rugnux/RugnuxFirstPass.cpp +++ b/rugnux/RugnuxFirstPass.cpp @@ -1472,6 +1472,10 @@ void Rugnux::FirstPassRun::NearTie(const std::string &decision, const std::strin 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); diff --git a/rugnux/RugnuxFirstPass.h b/rugnux/RugnuxFirstPass.h index 87bd1b85b..12710dbf3 100644 --- a/rugnux/RugnuxFirstPass.h +++ b/rugnux/RugnuxFirstPass.h @@ -150,7 +150,7 @@ struct Rugnux::FirstPassRun { 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. + // 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. From 7facb71cf16458ac55eaf079ec90b8f3e16e472f Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 9 Oct 2026 23:52:23 +0200 Subject: [PATCH 4/4] rugnux far-end check: each quantity's level, rest of the sweep against first part Testing every block against the other part's blocks fired on single odd blocks: a block's mean varies with the orientation, which a handful of reference blocks cannot represent (one block of 36 on a control, a 0.04 deg residual at t = 23 on another). The rotation rate, a straight-line fit per part, failed the same way, its standard error blind to the smooth orientation-dependent wiggle, and it did not see the target it was added for (the residual there runs at the same rate through both parts). Now each quantity is one comparison: the block means of the rest of the sweep against those of the first part, Welch's t, family-wise at 0.1 % over the quantities compared. Still log only. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi --- rugnux/RugnuxFarEnd.cpp | 178 ++++++++++++---------------------------- 1 file changed, 53 insertions(+), 125 deletions(-) diff --git a/rugnux/RugnuxFarEnd.cpp b/rugnux/RugnuxFarEnd.cpp index f379d5539..4d5dac496 100644 --- a/rugnux/RugnuxFarEnd.cpp +++ b/rugnux/RugnuxFarEnd.cpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include "Rugnux.h" @@ -17,14 +18,13 @@ // (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 every block past the first part is compared with the blocks inside -// it, and every block inside it with those past it. A block is measured against the other part's -// mean in units of the spread of that part's blocks themselves, as a prediction for one more block -// (Student's t on that part's degrees of freedom). A block is out of line when its two-sided p, times the number of late-block -// comparisons the run made over all the quantities, 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, however many blocks and quantities there are. No number in it is set -// by hand beyond the block width. +// 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 @@ -32,9 +32,7 @@ // - 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 - and for -// the rotation scale also its rate, the slope of the residual about the spindle against the -// angle, late blocks against early ones; +// 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 { @@ -108,17 +106,15 @@ void Rugnux::FarEndCheck(PipelineLocals &p, const RotationIndexerResult &rot) { return i < static_cast(v.size()) ? static_cast(v[i]) : NAN; }; - enum { INDEXED, ON_LATTICE, CELL, ABOUT_SPINDLE, OFF_SPINDLE, BACKGROUND, N_STATS, RATE = N_STATS }; + 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 + 1] = { + 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"}, - {"rate of the residual about the spindle (per 100 deg of rotation)", "the rotation scale", 100.0, "%", - "rate"}, }; std::vector> frame(N_STATS, std::vector(n, NAN)); const Coord spindle = gon.GetAxis().Normalize(); @@ -188,127 +184,59 @@ void Rugnux::FarEndCheck(PipelineLocals &p, const RotationIndexerResult &rot) { } const auto deg = [&](int ordinal) { return gon.GetAngle_deg(static_cast(ordinal)); }; - // Every late block of every quantity against the early blocks of the same quantity - and every - // early block against the late ones, because the part that is out of line can be the first: a - // crystal that settles, or an orientation that drifts, over the opening degrees makes the first - // part the odd one, and its spread then hides nothing the late blocks show. - struct Comparison { int stat; Range range; double value, mean, noise, t, p; bool first_part; }; + // 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; - const auto compare = [&](int s, const std::vector &reference, const std::vector &tested, - bool first_part) { - std::vector means; - for (const auto &r : reference) - if (const Block b = Summarise(frame[s], r.first, r.last); std::isfinite(b.mean)) - means.push_back(b.mean); - const int nr = static_cast(means.size()); - if (nr < 3) - return; - double mr = 0.0; - for (const double m : means) - mr += m; - mr /= nr; - double vr = 0.0; - for (const double m : means) - vr += (m - mr) * (m - mr); - vr /= nr - 1; - for (const auto &r : tested) { - const Block b = Summarise(frame[s], r.first, r.last); - if (!std::isfinite(b.mean)) - continue; - const double noise = std::sqrt(vr * (1.0 + 1.0 / nr) + b.se * b.se); - const double t = noise > 0.0 ? (b.mean - mr) / noise : (b.mean == mr ? 0.0 : INFINITY); - comparisons.push_back({s, r, b.mean, mr, noise, t, TwoSidedP(t, nr - 1), first_part}); - } - }; for (int s = 0; s < N_STATS; s++) { - compare(s, early_blocks, late_blocks, false); - compare(s, late_blocks, early_blocks, true); - } - // The rotation scale is a RATE: a stage turning at the wrong rate leaves a residual about the - // spindle that grows with the angle, and the scale decided on the first part zeroes its slope - // there. So the slope of the late blocks is compared with that of the early ones, each a - // straight-line fit to its blocks' means with its own standard error. - { - struct Fit { double slope = NAN, se = NAN; int blocks = 0; }; - const auto fit = [&](const std::vector &ranges) { - std::vector x, y; + const auto sample = [&](const std::vector &ranges) { + std::vector means; for (const auto &r : ranges) - if (const Block b = Summarise(frame[ABOUT_SPINDLE], r.first, r.last); std::isfinite(b.mean)) { - x.push_back(0.5 * (deg(r.first) + deg(r.last))); - y.push_back(b.mean); - } - Fit f; - f.blocks = static_cast(x.size()); - if (f.blocks < 3) - return f; - double mx = 0.0, my = 0.0; - for (int i = 0; i < f.blocks; i++) { - mx += x[i]; - my += y[i]; - } - mx /= f.blocks; - my /= f.blocks; - double sxx = 0.0, sxy = 0.0; - for (int i = 0; i < f.blocks; i++) { - sxx += (x[i] - mx) * (x[i] - mx); - sxy += (x[i] - mx) * (y[i] - my); - } - f.slope = sxy / sxx; - double rss = 0.0; - for (int i = 0; i < f.blocks; i++) { - const double e = y[i] - my - f.slope * (x[i] - mx); - rss += e * e; - } - f.se = std::sqrt(rss / (f.blocks - 2) / sxx); - return f; + 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 Fit fe = fit(early_blocks), fl = fit(late_blocks); - if (fe.blocks >= 3 && fl.blocks >= 3) { - const double noise = std::sqrt(fe.se * fe.se + fl.se * fl.se); - const double t = noise > 0.0 ? (fl.slope - fe.slope) / noise : 0.0; - comparisons.push_back({RATE, {late_blocks.front().first, late_blocks.back().last}, fl.slope, - fe.slope, noise, t, TwoSidedP(t, fe.blocks + fl.blocks - 4), false}); - } + 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, int out_of_line, int late) { + const auto describe = [&](const Comparison &c) { const Stat &st = stats[c.stat]; - return fmt::format("{} over {:.0f}-{:.0f} deg is {:.3g}{} against {:.3g}{} +- {:.2g} on {} " - "(t = {:.1f}, p = {:.2g} in a family of {}; {} of {} blocks out of line)", - st.what, deg(c.range.first), deg(c.range.last), c.value * st.unit, st.suffix, - c.mean * st.unit, st.suffix, c.noise * st.unit, - c.first_part ? "the rest of the sweep" : "the first part", c.t, c.p, - comparisons.size(), out_of_line, late); + 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 smallest p times the family, for the record + std::string by_quantity; // each quantity's p times the family, for the record const Comparison *closest = nullptr; - for (int s = 0; s <= RATE; s++) { - const Comparison *worst = nullptr; - int out_of_line = 0, tested = 0; - for (const auto &c : comparisons) - if (c.stat == s) { - tested++; - if (c.p * family < alpha) - out_of_line++; - if (!worst || c.p < worst->p) - worst = &c; - } - if (!worst) - continue; - by_quantity += fmt::format("{}{} {:.2g}", by_quantity.empty() ? "" : ", ", stats[s].tag, - std::min(1.0, worst->p * family)); - if (out_of_line > 0) - revisits.push_back(fmt::format("{} would be revisited because {}", stats[s].decision, - describe(*worst, out_of_line, tested))); - if (!closest || worst->p < closest->p) - closest = worst; + 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 auto late_of = [&](int s) { - return static_cast(std::count_if(comparisons.begin(), comparisons.end(), - [s](const Comparison &c) { return c.stat == s; })); - }; const std::string blocks = fmt::format("{} early and {} late blocks of {:.0f} deg", early_blocks.size(), late_blocks.size(), block * step); @@ -317,7 +245,7 @@ void Rugnux::FarEndCheck(PipelineLocals &p, const RotationIndexerResult &rot) { deg(early)); else if (revisits.empty()) logger.Info("far-end check: early decisions hold ({}; closest: {}; p x family by quantity: {})", - blocks, describe(*closest, 0, late_of(closest->stat)), by_quantity); + blocks, describe(*closest), by_quantity); else { std::string line; for (const auto &r : revisits)