From 7facb71cf16458ac55eaf079ec90b8f3e16e472f Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 9 Oct 2026 23:52:23 +0200 Subject: [PATCH] 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)