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

RotationScaleMerge: the partials' credible-scale floor is a fraction of the typical frame, not of the median (RESULT-CHANGING)

Browse Source 
The floor below which a frame's partial scale is not credible (MIN_CREDIBLE_SCALE_RATIO = 0.02) was
read off the MEDIAN fitted scale - in the loop's smoothing window, in its step, and in the guard
after it (DropCollapsedScales). On a sweep that spent most of its turn out of the beam the median is
itself a dead frame, and it walks: an empty frame's scale has no fixed point above zero, so the dead
frames slide toward zero round after round, the median and the floor slide with them, and frames on
the floor cross it, re-enter the window at their neighbours' scale (3x up in one round) and cross it
again. The loop then has no fixed point, and which frames reach the merge depends on the round it
stops on.

The floor is now read off the run's typical frame (TypicalFrameScale), which is what the loop pins
its gauge to - inside the loop it is the constant 0.02, so it stands still. The fulls' guard keeps
the median (their scales are not gauge-pinned). Dropped frames still carry the median scale, so the
sweep ledger reads them as before.

lyso_x06da_half_image (1800 frames, two thirds weak or out of the beam; median frame at 1/5.7 of the
typical one): traced, the final loop's floor fell 2.4x and every spike in its step was one frame
re-entering the window; the R_meas 287% -> 685% of the rc178 loop fixes was the same non-settling
loop ended by two stopping rules (delta-CC1/2 split the weak stretch differently, 206 frames at
1/20-1/50 of the median merged or not). The pre-pass P1 loop that crept to the cap is the geometry
pre-pass, which commits nothing here.

battery --only, before (rc178-aa-tol5 f217f180a) -> after, GPU:
- lyso_x06da_half_image: R_meas 689% -> 81%, CC1/2 0.932 -> 0.971, REFRES CC1/2 0.484 -> 0.911,
  ISa 4.13 -> 4.22, model R-free 0.3302 -> 0.3271, CC_model 0.826 -> 0.829; final loop's smallest
  step 5.6e-3 -> 1.8e-4. --no-gpu gives the same numbers.
- lcystine_x10sa_25keV: R_meas 0.396 -> 0.202, ISa 4.45 -> 8.49, CC1/2 0.9987 -> 0.9998.
- lyso_x10sa_strong: R_meas 0.238 -> 0.151, ISa 5.51 -> 5.74, R-free 0.3653 -> 0.3658,
  CC_model 0.859 -> 0.857, d_min 1.354 -> 1.344.
- insu_I_x06da_ref: R_meas 0.347 -> 0.290, R-free and ISa unchanged.
- 7mzt: CC1/2 0.967 -> 0.983, CC_model 0.532 -> 0.548, R-free 0.4183 -> 0.4197, d_min 3.18 -> 3.21.
- unchanged: lyso_x06da_ref, myob_x06da_split, myob_x06da_powder_2, thau_bl1a_4p6keV, 6o2h,
  3meb, 9z72, 9o0h. Space groups identical on all 13.
Tests: [RotationScale], AndersonMixing* pass.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01PXc1zD1yxHqaR1DGdCTHL4
This commit is contained in:
leonarski_fandClaude Opus 5.5 committed 2026-10-20 02:18:19 +02:00
1 parent 1949155c12
commit 2a5d303886
2 files changed
+41 -40
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

image_analysis/scale_merge/RotationScaleMerge.cpp
+38 -37
View File
@@ -83,15 +83,17 @@ namespace {
// Interleaved folds of the cross-validation that sets the smoothness of the fulls' scale.
constexpr int CV_FOLDS = 5;
// A fitted per-frame scale below this fraction of the run's median is not a measurement of anything:
// it would mean the frame received 2% of its neighbours' dose while still producing indexable spots.
// A fitted per-frame scale below this fraction of the run's typical frame (of the median, for the
// fulls - see DropCollapsedScales) is not a measurement of anything: it would mean the frame received
// 2% of its neighbours' dose while still producing indexable spots.
// It matters because the scale enters as 1/G, so a G that collapses toward zero multiplies every
// intensity on that frame without bound - and sigma by the identical factor, which is why no
// sigma-based outlier test can see it. Measured over 12 rotation crystals in the default
// configuration, the smallest LEGITIMATE min(G)/median(G) is 0.070; the case this bound exists for
// sat at 0.0017 (a 546x amplification that put 99% of sum(I^2) into 0.66% of the observations and
// took the merged CC1/2 from 93.7% to 17.2%). 0.02 is 3.5x below anything real and 12x above the
// failure.
// failure. Against the typical frame, over the 252 merged sets of one battery (read off their final
// per-frame scales), the smallest merged frame of a set that drops none sat at 0.036.
constexpr double MIN_CREDIBLE_SCALE_RATIO = 0.02;
// The same quantity, for the P1 merge that feeds the de-novo space-group search only. That merge
@@ -1999,9 +2001,13 @@ RotationScaleMerge::ScalingLoopOutcome RotationScaleMerge::RunScalingLoop(
// What each frame's corr is multiplied by: the fitted G over the G it ends the round with.
for (int f = 0; f < n_frames; ++f)
ratio[f] = scaled_any[f] && std::isfinite(g[f]) && g[f] > 0.0 ? g[f] : 1.0;
// The credible floor is a fraction of the run's typical frame - the frame the gauge is pinned
// to below, so the floor stands still while the loop runs (see DropCollapsedScales).
for (int f = 0; f < n_frames; ++f)
n_obs[f] = scaled_any[f] && std::isfinite(g[f]) && g[f] > 0.0 ? frame_obs_count[f] : 0;
if (window > 1) {
frame_scaled_scratch = scaled_any;
ComputeSmoothGWindow(g, window, g_smooth);
ComputeSmoothGWindow(g, window, TypicalFrameScale(g, n_obs) * MIN_CREDIBLE_SCALE_RATIO, g_smooth);
for (int f = 0; f < n_frames; ++f)
if (scaled_any[f] && std::isfinite(g_smooth[f]) && g_smooth[f] > 0.0)
g[f] = g_smooth[f];
@@ -2020,15 +2026,9 @@ RotationScaleMerge::ScalingLoopOutcome RotationScaleMerge::RunScalingLoop(
// leaves it out (ComputeSmoothGWindow) and the guard after the loop drops it. Left in the
// step, a single such frame held the rms at 0.18 and the loop was stopped as unsettled after
// 16 rounds; with it left out the same loop settles in 23. On that sweep every loop that
// stopped unsettled was one that went on to drop a frame.
std::vector<double> fitted;
for (int f = 0; f < n_frames; ++f)
if (scaled_any[f] && g[f] > 0.0) fitted.push_back(g[f]);
double g_floor = 0.0;
if (!fitted.empty()) {
std::nth_element(fitted.begin(), fitted.begin() + fitted.size() / 2, fitted.end());
g_floor = fitted[fitted.size() / 2] * MIN_CREDIBLE_SCALE_RATIO;
}
// stopped unsettled was one that went on to drop a frame. The scales are pinned to the
// typical frame by now, so the floor is the ratio itself.
const double g_floor = MIN_CREDIBLE_SCALE_RATIO;
double ss = 0.0, sw = 0.0;
int n = 0;
for (int f = 0; f < n_frames; ++f) {
@@ -4439,7 +4439,7 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector<int32_t> &cell, int
rounds_used, gauge_ramp, n_clamped, ncell);
}
void RotationScaleMerge::ComputeSmoothGWindow(const std::vector<double> &g, int window,
void RotationScaleMerge::ComputeSmoothGWindow(const std::vector<double> &g, int window, double g_floor,
std::vector<double> &g_smooth) const {
const int n = static_cast<int>(g.size());
const int half = window / 2;
@@ -4453,19 +4453,8 @@ void RotationScaleMerge::ComputeSmoothGWindow(const std::vector<double> &g, int
// stays NAN, which the caller reads as "leave this frame alone"), so what reaches
// DropCollapsedScales is the scale the fit actually produced.
//
// The floor is the guard's, read the same way: MIN_CREDIBLE_SCALE_RATIO of the MEDIAN fitted
// scale, not of the precision-weighted typical frame - see DropCollapsedScales for why the median.
std::vector<double> fitted;
fitted.reserve(n);
for (int f = 0; f < n; ++f)
if (frame_scaled_scratch[f] && std::isfinite(g[f]) && g[f] > 0.0)
fitted.push_back(g[f]);
if (fitted.empty())
return;
const size_t mid = fitted.size() / 2;
std::nth_element(fitted.begin(), fitted.begin() + mid, fitted.end());
const double g_floor = fitted[mid] * MIN_CREDIBLE_SCALE_RATIO;
// The floor is the guard's, read the same way: MIN_CREDIBLE_SCALE_RATIO of the run's typical frame,
// passed in by the loop.
// Each frame's logarithm once, not once per window it falls in.
std::vector<double> log_g(n, 0.0);
std::vector<uint8_t> credible(n, 0);
@@ -4488,7 +4477,7 @@ void RotationScaleMerge::ComputeSmoothGWindow(const std::vector<double> &g, int
}
}
// Drop the observations of any frame whose scale fitted far below the run's median, reporting the
// Drop the observations of any frame whose scale fitted far below the run's typical frame, reporting the
// per-frame factor the caller has to apply to its corr (1 where nothing changed, 0 where the frame is
// dropped). Returns false when there was nothing to drop.
//
@@ -4510,11 +4499,19 @@ bool RotationScaleMerge::DropCollapsedScales(const std::vector<uint8_t> &fitted_
apply.assign(n, 0);
ratio.assign(n, 1.0);
// Against the MEDIAN frame, not the precision-weighted typical frame the sweep ledger reads: a
// frame that collapsed toward zero and was not quite dropped has its intensities amplified by
// 1/G, the fulls built from them are re-fitted with a scale of that size, and a cubic mean is then
// theirs - measured on a sweep 60% out of the beam, the fulls' precision-weighted typical scale
// sat on such frames and the floor dropped three quarters of every live frame's fulls.
// Against the run's typical frame (TypicalFrameScale), the frame the scaling loop pins its gauge to,
// not against the median frame. On a sweep that spent most of its turn out of the beam the median is
// itself a dead frame, and it walks: the dead frames' scales slide toward zero round after round
// (an empty frame's scale has no fixed point above zero), the median and the floor slide with them,
// and the frames on the floor cross it, re-enter the smoothing window at their neighbours' scale and
// cross it again - so the loop never settled, and which frames were merged depended on the round it
// stopped on. Measured on a lysozyme sweep two thirds out of the beam, where the median frame sat at
// a sixth of the typical one: the floor fell 2.4x in one loop, and the same loop ended by two
// stopping rules merged R_meas 287% and 685%; with the floor on the typical frame it merges 81%,
// CC1/2 0.93 -> 0.97, model R-free 0.330 -> 0.327. Where most of the sweep diffracts the
// median and the typical frame are the same frame and nothing changes. The fulls keep the median
// (DropCollapsedFullScales): their scales are not gauge-pinned, and amplified fulls of a frame that
// was not quite dropped pull a cubic mean up onto themselves.
std::vector<double> fitted;
fitted.reserve(n);
for (int f = 0; f < n; ++f)
@@ -4525,7 +4522,11 @@ bool RotationScaleMerge::DropCollapsedScales(const std::vector<uint8_t> &fitted_
const size_t mid = fitted.size() / 2;
std::nth_element(fitted.begin(), fitted.begin() + mid, fitted.end());
const double g_typ = fitted[mid];
const double g_median = fitted[mid];
std::vector<int64_t> n_obs(n, 0);
for (int f = 0; f < n; ++f)
if (fitted_mask[f] && std::isfinite(g[f]) && g[f] > 0.0) n_obs[f] = frame_count[f];
const double g_typ = TypicalFrameScale(g, n_obs);
const double g_floor = g_typ * MIN_CREDIBLE_SCALE_RATIO;
int n_dropped = 0;
@@ -4536,12 +4537,12 @@ bool RotationScaleMerge::DropCollapsedScales(const std::vector<uint8_t> &fitted_
apply[f] = 1;
ratio[f] = 0.0; // zero corr = out of the combine, the merge and the error model
worst = std::max(worst, g_typ / g[f]);
g[f] = g_typ; // keep the median in g so this frame cannot drag a smoothing window
g[f] = g_median; // keep the median in g so this frame cannot drag a smoothing window
++n_dropped;
}
if (n_dropped > 0)
logger->Warning("Dropped {} frame(s) whose scale came out more than {:.0f}x below the run median "
"(worst {:.0f}x) - a scale that small says the frame holds no diffraction, and "
logger->Warning("Dropped {} frame(s) whose scale came out more than {:.0f}x below the run's typical "
"frame (worst {:.0f}x) - a scale that small says the frame holds no diffraction, and "
"using it would amplify the frame's intensities, and its sigmas with them, "
"by the same factor",
n_dropped, 1.0 / MIN_CREDIBLE_SCALE_RATIO, worst);
image_analysis/scale_merge/RotationScaleMerge.h
+3 -3
View File
@@ -883,11 +883,11 @@ private:
// The windowed geometric mean of G over frames; the scaling loop applies the resulting ratio to
// the observations' corr (on the GPU, to the resident corr in a kernel).
void ComputeSmoothGWindow(const std::vector<double> &g, int window,
void ComputeSmoothGWindow(const std::vector<double> &g, int window, double g_floor,
std::vector<double> &g_smooth) const;
// Drop the observations of any frame whose fitted per-frame scale collapsed far below the run
// median, reporting the per-frame corr factor the caller has to apply. See the .cpp for why nothing
// Drop the observations of any frame whose fitted per-frame scale collapsed far below the run's
// typical frame, reporting the per-frame corr factor the caller has to apply. See the .cpp for why nothing
// downstream can catch a collapsed scale on its own.
bool DropCollapsedScales(const std::vector<uint8_t> &fitted_mask, std::vector<double> &g,
std::vector<uint8_t> &apply, std::vector<double> &ratio) const;