diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 70cc2f6a0..3b9fbdbe3 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -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 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 &cell, int rounds_used, gauge_ramp, n_clamped, ncell); } -void RotationScaleMerge::ComputeSmoothGWindow(const std::vector &g, int window, +void RotationScaleMerge::ComputeSmoothGWindow(const std::vector &g, int window, double g_floor, std::vector &g_smooth) const { const int n = static_cast(g.size()); const int half = window / 2; @@ -4453,19 +4453,8 @@ void RotationScaleMerge::ComputeSmoothGWindow(const std::vector &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 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 log_g(n, 0.0); std::vector credible(n, 0); @@ -4488,7 +4477,7 @@ void RotationScaleMerge::ComputeSmoothGWindow(const std::vector &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 &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 fitted; fitted.reserve(n); for (int f = 0; f < n; ++f) @@ -4525,7 +4522,11 @@ bool RotationScaleMerge::DropCollapsedScales(const std::vector &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 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 &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); diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index e3b4ab0d1..4ef64c4fb 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -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 &g, int window, + void ComputeSmoothGWindow(const std::vector &g, int window, double g_floor, std::vector &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 &fitted_mask, std::vector &g, std::vector &apply, std::vector &ratio) const;