mx/Jungfraujoch
3
2
Fork 0
Code Issues Pull Requests Actions 2 Packages Projects Releases 81 Wiki Activity

rugnux far-end check: each quantity's level, rest of the sweep against first part

Browse Source 
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) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi
This commit is contained in:
leonarski_f
2026-10-09 23:52:23 +02:00
co-authored by Claude Opus 5.5
parent 663cb2b276
commit 7facb71cf1
1 changed files with 53 additions and 125 deletions
Download Patch File Download Diff File Expand all files Collapse all files
rugnux/RugnuxFarEnd.cpp
+53 -125
View File
@@ -4,6 +4,7 @@
#include <algorithm>
#include <cmath>
#include <string>
#include <tuple>
#include <vector>
#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<int>(v.size()) ? static_cast<double>(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<std::vector<double>> frame(N_STATS, std::vector<double>(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<float>(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<Comparison> comparisons;
const auto compare = [&](int s, const std::vector<Range> &reference, const std::vector<Range> &tested,
bool first_part) {
std::vector<double> 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<int>(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<Range> &ranges) {
std::vector<double> x, y;
const auto sample = [&](const std::vector<Range> &ranges) {
std::vector<double> 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<int>(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<size_t>(1, means.size());
for (const double x : means)
v += (x - m) * (x - m);
v /= std::max<size_t>(2, means.size()) - 1;
return std::tuple<double, double, int>{m, v, static_cast<int>(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<int>(dof)))});
}
const double alpha = std::erfc(SPOT_BUDGET_SIGNIFICANCE_Z / std::sqrt(2.0));
const double family = static_cast<double>(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<std::string> 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<int>(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)