From c4d74dd021a338e6934ec653d82d5da269ba56f4 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 20 Sep 2026 10:57:35 +0200 Subject: [PATCH] Rotation scaling: the per-frame scale is restrained to be smooth over the rotation inside the loop The alternating scale fit has a second degenerate direction beyond the gauge, and that one is not harmless. A frame whose scale drifts below its neighbours' loses its vote in the references it is fitted against - an observation's weight there goes as G^2 - which moves the references away from it, which drives its scale further down. Where every frame is measured against hundreds of reflections the references barely notice one frame and the loop settles; on a fine-sliced sweep of a small cell, where a frame holds two or three reflections and a rocking curve spans a dozen or more frames, the frames at the tails of the curves lose against the frames at their peaks and the loop walks the bulk of the sweep down fifteen decades. The convergence test then said the loop had settled, because its own weight (frame observations at G^2) vanishes for exactly the frames that are collapsing. So the scale is now smoothed over the rotation INSIDE the loop, at the window the merge already smoothed G over afterwards (the post-loop smoothing is the same operation and is not repeated): every frame in a window shares the window's geometric mean, no frame can lose its vote alone, and what a frame's own observations cannot determine its neighbourhood's do. The window is floored at six rocking curves - within one curve a change of scale and an error of the partiality model are the same thing, and a window holding a curve or two fits the model's error as scale. The step test is weighted by each frame's observation count alone, and a loop whose step has not fallen below nine tenths of its smallest value for five rounds stops and says so rather than walking further. Measured against the previous behaviour: a fine-sliced sweep that collapsed to P1 with CC1/2 0.000 merges in its true orthorhombic group at CC1/2 0.95; a monoclinic small-cell set goes from an unusable 2.13 A merge to 0.81 A at R_meas 11.5%; another from R_meas 33% / ISa 2.2 to 6.1% / 14.2; two powder-ring sets and a short-sweep tetragonal one improve. The space group is unchanged on every protein control measured, and the short-sweep and de-novo orthorhombic gains of the convergent loop are kept. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_013nW6FNRP1bBJJ8pfHiByAT --- .../scale_merge/RotationScaleMerge.cpp | 163 ++++++++++-------- .../scale_merge/RotationScaleMerge.h | 17 +- rugnux/ResultReport.cpp | 6 +- 3 files changed, 101 insertions(+), 85 deletions(-) diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 8bf050f1c..1a8fa8039 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -38,6 +38,8 @@ namespace { // The per-frame scaling loop has settled when the scales moved less than this between two // iterations, as the rms of |log(G_new/G_old)| over the frames fitted in both. constexpr double SCALING_TOLERANCE = 1e-3; + // The per-frame scale is smoothed over at least this many rocking curves (see Run). + constexpr int SMOOTH_G_MIN_ROCKING_EVENTS = 6; // 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. @@ -1408,7 +1410,7 @@ void RotationScaleMerge::UpdateCorr(std::vector &obs, const std::vector &g, const std::vector &frame_obs_count, + const char *what, std::vector &g, const std::vector &frame_obs_count, int window, const std::function &iterate, const std::function &, const std::vector &)> &rescale) { // Alternating least squares on a bilinear model: the reference intensities given the scales, the @@ -1417,41 +1419,66 @@ RotationScaleMerge::ScalingLoopOutcome RotationScaleMerge::RunScalingLoop( // every round, so the scales never stop moving however many rounds they are given. The gauge is // therefore pinned after every iteration: G is divided by the run's typical frame scale (the // precision-weighted G_ref every other frame measure is read against) and the observations' corr - // multiplied by it, which changes no relative scale and no merged intensity. What is left to - // converge is the relative scales, and the loop stops when they stop moving. + // multiplied by it, which changes no relative scale and no merged intensity. + // + // It has a second degenerate direction, and that one is not harmless: a frame whose scale drifts + // below its neighbours' loses its vote in the references it is fitted against (an observation's + // weight there goes as G^2), which moves the references away from it, which drives its scale + // further down. On a merge where every frame is measured against hundreds of reflections the + // references barely notice one frame and the loop settles; on a fine-sliced small-molecule sweep, + // where a frame holds two or three reflections and a rocking curve spans a dozen frames, the + // frames at the tails of the rocking curves lost against the frames at their peaks and the loop + // walked the bulk of the sweep down fifteen decades while reporting itself settled. So the scale + // is restrained to be smooth over the rotation INSIDE the loop, at the window the merge already + // smooths G over afterwards (XDS's DELPHI): every frame in a window shares the window's geometric + // mean, so no frame can lose its vote alone, and what a frame's own observations cannot determine + // its neighbourhood's do. The post-loop smoothing is then this same operation and is not repeated. + // + // The loop stops when the relative scales stop moving - the rms of |log(G_new/G_old)| over the + // frames, each weighted by its observation count and nothing else (a weight that vanishes with G + // is what let a collapsing sweep report itself settled) - or when they stop settling: a step that + // has not fallen below nine tenths of its smallest value for five rounds is a loop on a flat or + // rising valley, and more rounds of it can only walk further. Both leave the frames where they are + // and say so. ScalingLoopOutcome out; std::vector scaled_any(n_frames, 0); - std::vector g_prev(n_frames, 0.0); + std::vector g_prev(n_frames, 0.0), g_smooth; std::vector n_obs(n_frames, 0); std::vector ratio(n_frames, 1.0); + double best_step = INFINITY; + int rounds_since_best = 0; for (int it = 0; it < scaling_iter; ++it) { iterate(); out.iterations = it + 1; - for (int f = 0; f < n_frames; ++f) { + for (int f = 0; f < n_frames; ++f) scaled_any[f] |= frame_scaled_scratch[f]; - n_obs[f] = scaled_any[f] && std::isfinite(g[f]) && g[f] > 0.0 ? frame_obs_count[f] : 0; + // 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; + if (window > 1) { + frame_scaled_scratch = scaled_any; + ComputeSmoothGWindow(g, window, 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]; } + 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; const double g_ref = TypicalFrameScale(g, n_obs); for (int f = 0; f < n_frames; ++f) { - ratio[f] = scaled_any[f] ? g_ref : 1.0; - if (scaled_any[f]) g[f] /= g_ref; + if (!scaled_any[f] || !(std::isfinite(g[f]) && g[f] > 0.0)) { ratio[f] = 1.0; continue; } + g[f] /= g_ref; + ratio[f] /= g[f]; } rescale(scaled_any, ratio); if (it > 0) { - // The step, weighted by what each frame carries into the merge - its observations at its - // scale squared, which is how the inverse-variance weight of a scaled observation goes. - // A frame the crystal did not diffract on has its scale fitted on noise and wanders by - // orders of magnitude every iteration (measured on a sweep 60% out of the beam: scales - // at 1e-7 of the typical frame, still moving at the cap); it carries nothing into the - // merge and is dropped from it after the loop, so it does not hold the loop open either. double ss = 0.0, sw = 0.0; int n = 0; for (int f = 0; f < n_frames; ++f) if (scaled_any[f] && g_prev[f] > 0.0 && g[f] > 0.0) { const double d = std::log(g[f] / g_prev[f]); - const double w = frame_obs_count[f] * g[f] * g[f]; - ss += w * d * d; - sw += w; + ss += frame_obs_count[f] * d * d; + sw += frame_obs_count[f]; ++n; } out.step = sw > 0.0 ? std::sqrt(ss / sw) : 0.0; @@ -1461,6 +1488,13 @@ RotationScaleMerge::ScalingLoopOutcome RotationScaleMerge::RunScalingLoop( out.converged = true; break; } + if (out.step < 0.9 * best_step) { + best_step = out.step; + rounds_since_best = 0; + } else if (++rounds_since_best >= 5) { + out.stalled = true; + break; + } } g_prev = g; } @@ -1468,6 +1502,11 @@ RotationScaleMerge::ScalingLoopOutcome RotationScaleMerge::RunScalingLoop( if (out.converged) logger.Info("Per-frame scaling ({}): settled after {} iterations (rms |dlogG| {:.1e})", what, out.iterations, out.step); + else if (out.stalled) + logger.Warning("Per-frame scaling ({}): stopped after {} iterations - the scales stopped settling " + "(rms |dlogG| {:.1e}, not below {:.1e} for five rounds); more rounds would only walk " + "them, so everything read off this merge is read off an unsettled state", what, + out.iterations, out.step, best_step); else logger.Warning("Per-frame scaling ({}): NOT settled after the cap of {} iterations - the scales " "were still moving by rms |dlogG| {:.1e} (tolerance {:.0e}); everything read off " @@ -3182,11 +3221,14 @@ void RotationScaleMerge::ComputeSmoothGWindow(const std::vector &g, int const int half = window / 2; g_smooth.assign(n, NAN); - // No collapsed-scale handling here: DropCollapsedScales has already run on g, so every fitted - // entry is credible. That matters because the window average is a GEOMETRIC mean - a scale - // collapsing toward zero would not just corrupt its own frame, its logarithm would drag the whole - // window down, and where the only fitted frames in a window ARE the collapsed ones the mean would - // average the fault with itself and make it permanent. + // No collapsed-scale handling here, and the caller is the scaling loop, which runs BEFORE + // DropCollapsedScales - so a frame the crystal did not diffract on is in the window with the rest. + // The window average is a GEOMETRIC mean, so such a frame drags its whole window down by its + // logarithm, and is itself lifted toward its neighbours. That is the point on a fine-sliced sweep + // (a frame with two reflections cannot be told from a dead one by its own scale) and it is the + // price on a sweep with a genuinely dead stretch: the stretch is no longer far enough below the + // median for DropCollapsedScales to see it, and what removes it is the sweep ledger's + // delta-CC1/2, which asks what the frames cost rather than how small their scale came out. for (int o = 0; o < n; ++o) { double sum_log = 0.0; int count = 0; @@ -3260,22 +3302,6 @@ bool RotationScaleMerge::DropCollapsedScales(const std::vector &fitted_ return n_dropped > 0; } -void RotationScaleMerge::SmoothG(std::vector &obs, std::vector &g, int window) const { - const int n = static_cast(g.size()); - std::vector g_smooth; - ComputeSmoothGWindow(g, window, g_smooth); - for (auto &o : obs) { - const int f = o.frame; - if (!frame_scaled_scratch[f] || !std::isfinite(g[f]) || g[f] <= 0.0 || !std::isfinite(g_smooth[f])) - continue; - if (std::isfinite(o.corr)) - o.corr = static_cast(o.corr * (g[f] / g_smooth[f])); - } - for (int f = 0; f < n; ++f) - if (frame_scaled_scratch[f] && std::isfinite(g[f]) && g[f] > 0.0 && std::isfinite(g_smooth[f])) - g[f] = g_smooth[f]; -} - // The same guard on the scale-fulls pass, which refits G from scratch and so can collapse a frame the // partials guard passed. Detection and remedy are both DropCollapsedScales': the frame's fulls leave the // merge. Merging them unscaled - the old behaviour - asserts G = 1 for a frame whose scale is @@ -4717,6 +4743,27 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st std::vector partial_mean; bool scaled_on_gpu = false; + // The smoothing window (--smooth-g, XDS DELPHI-like), applied inside both scaling loops - and + // never narrower than SMOOTH_G_MIN_ROCKING_EVENTS rocking curves. Within one rocking curve a + // change of scale and an error of the partiality model are the same thing, so a window that + // holds only a curve or two fits the model's error as scale: measured on a fine-sliced + // small-molecule sweep whose curves span 17 frames, the 5 deg window (a curve and a half) let + // the loop walk to a 140x spread of scales across the sweep, three curves to 25x, six curves to + // a settled 6x. On protein data at ordinary mosaicity 5 deg already holds ten or more curves and + // the floor changes nothing. + const auto s = x.GetScalingSettings(); + const double smooth_g_deg = s.GetSmoothGDegrees(); + const auto gonio = x.GetGoniometer(); + const double osc_deg = gonio ? std::fabs(gonio->GetIncrement_deg()) : 0.0; + int smooth_window = 1; + if (smooth_g_deg > 0.0 && osc_deg > 1e-6) { + smooth_window = std::max(static_cast(std::lround(smooth_g_deg / osc_deg)), + SMOOTH_G_MIN_ROCKING_EVENTS * RockingEventFrames()); + if (smooth_window % 2 == 0) ++smooth_window; + logger.Info("Per-frame scale smoothing window: {} frames ({:.1f} deg; {} deg asked, rocking " + "curves span {} frames)", smooth_window, smooth_window * osc_deg, smooth_g_deg, + RockingEventFrames()); + } ScalingLoopOutcome partial_loop; #ifdef JFJOCH_USE_CUDA if (gpu_active_) { @@ -4727,7 +4774,7 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st // fit this iteration keeps the G it had (its corr was left where it was). std::vector g_dev(n_frames); std::vector scaled_dev(n_frames); - partial_loop = RunScalingLoop("partials", g_partial, frame_count, [&] { + partial_loop = RunScalingLoop("partials", g_partial, frame_count, smooth_window, [&] { gpu_->ScalePartials(1, min_partiality, d_min_limit.has_value()); gpu_->GetG(g_dev.data(), scaled_dev.data()); for (int f = 0; f < n_frames; ++f) @@ -4740,7 +4787,7 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st } #endif if (!scaled_on_gpu) { - partial_loop = RunScalingLoop("partials", g_partial, frame_count, [&] { + partial_loop = RunScalingLoop("partials", g_partial, frame_count, smooth_window, [&] { ReduceGroupMeans(partials, n_groups, partial_mean); FitPerFrameG(partials, frame_start, frame_count, partial_mean, /*unity=*/false, g_partial); UpdateCorr(partials, g_partial, frame_scaled_scratch); @@ -4790,37 +4837,7 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st } } - // --- 2. Smooth G across frames (XDS DELPHI-like) before the combine. --- - const auto s = x.GetScalingSettings(); - const double smooth_g_deg = s.GetSmoothGDegrees(); - const auto gonio = x.GetGoniometer(); - const double osc_deg = gonio ? std::fabs(gonio->GetIncrement_deg()) : 0.0; - if (smooth_g_deg > 0.0 && osc_deg > 1e-6) { - int window = std::max(1, static_cast(std::lround(smooth_g_deg / osc_deg))); - if (window % 2 == 0) ++window; - bool smoothed_on_gpu = false; -#ifdef JFJOCH_USE_CUDA - if (gpu_active_) { - // Apply smooth-G to the resident corr: the host builds the per-frame ratio g/g_smooth (+ an - // apply flag), the kernel multiplies corr in place. Same guard / arithmetic as host SmoothG. - std::vector g_smooth; - ComputeSmoothGWindow(g_partial, window, g_smooth); - std::vector apply(n_frames, 0); - std::vector ratio(n_frames, 1.0); - for (int f = 0; f < n_frames; ++f) - if (frame_scaled_scratch[f] && std::isfinite(g_partial[f]) && g_partial[f] > 0.0 - && std::isfinite(g_smooth[f])) { - apply[f] = 1; - ratio[f] = g_partial[f] / g_smooth[f]; - } - gpu_->SmoothCorr(apply.data(), ratio.data()); - for (int f = 0; f < n_frames; ++f) if (apply[f]) g_partial[f] = g_smooth[f]; - smoothed_on_gpu = true; - } -#endif - if (!smoothed_on_gpu) - SmoothG(partials, g_partial, window); - } + // --- 2. G is smoothed across frames (XDS DELPHI-like) inside the scaling loop above. --- // Per-frame CC + write G/CC/mosaicity back onto the partials (once). On the GPU the group means + // per-frame CC run on the resident (already smoothed) corr; only the tiny per-frame cc/cc_n come back. @@ -4992,7 +5009,7 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st if (static_cast(g_full.size()) != n_frames) g_full.assign(n_frames, 1.0); std::vector g_dev(n_frames); std::vector scaled_dev(n_frames); - fulls_loop = RunScalingLoop("fulls", g_full, f_count, [&] { + fulls_loop = RunScalingLoop("fulls", g_full, f_count, smooth_window, [&] { gpu_->ScaleFulls(1, min_partiality); gpu_->GetG(g_dev.data(), scaled_dev.data()); for (int f = 0; f < n_frames; ++f) @@ -5052,7 +5069,7 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st // --- 4. Scale the fulls (XDS order, Unity model). --- if (scale_fulls && !scaled_fulls_on_gpu) { std::vector full_mean; - fulls_loop = RunScalingLoop("fulls", g_full, fulls_frame_count, [&] { + fulls_loop = RunScalingLoop("fulls", g_full, fulls_frame_count, smooth_window, [&] { ReduceGroupMeans(fulls, n_groups, full_mean); FitPerFrameG(fulls, fulls_frame_start, fulls_frame_count, full_mean, /*unity=*/true, g_full); UpdateCorr(fulls, g_full, frame_scaled_scratch); diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 6cc953c67..a07017460 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -445,12 +445,13 @@ private: // The alternating per-frame scale fit run to convergence: `iterate` fits every frame once against // the current reference (leaving G in `g` and the frames it fitted in frame_scaled_scratch), - // `rescale` multiplies the corr of the flagged frames' observations by a per-frame ratio. Pins the - // gauge after every iteration and stops when the scales settle or the cap is reached. On return - // frame_scaled_scratch flags every frame fitted in any iteration. - struct ScalingLoopOutcome { int iterations = 0; bool converged = false; double step = 0.0; }; + // `rescale` multiplies the corr of the flagged frames' observations by a per-frame ratio. Smooths G + // over `window` frames and pins the gauge after every iteration; stops when the scales settle, stop + // settling, or the cap is reached. On return frame_scaled_scratch flags every frame fitted in any + // iteration. + struct ScalingLoopOutcome { int iterations = 0; bool converged = false; bool stalled = false; double step = 0.0; }; ScalingLoopOutcome RunScalingLoop(const char *what, std::vector &g, - const std::vector &frame_obs_count, + const std::vector &frame_obs_count, int window, const std::function &iterate, const std::function &, const std::vector &)> &rescale); @@ -459,10 +460,8 @@ private: void UpdateCorr(std::vector &obs, const std::vector &g, const std::vector &frame_scaled) const; - void SmoothG(std::vector &obs, std::vector &g, int window) const; - - // The windowed geometric mean of G over frames (the shared first half of SmoothG); the GPU path - // applies the resulting ratio to the resident corr in a kernel instead of the host obs loop. + // 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, std::vector &g_smooth) const; diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 0243ff06e..ee1d83a79 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -754,9 +754,9 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, result.scaling_converged ? "" : " (NOT settled)"), true)); if (!result.scaling_converged || !result.search_scaling_converged) Warn(doc, PathologyCode::SCALING_NOT_CONVERGED, fmt::format( - "the per-frame scales of {} were still moving when the iteration cap stopped " - "them, so everything read off that merge is read off an unsettled state " - "(--scaling-iterations raises the cap)", + "the per-frame scales of {} were still moving when the scaling loop stopped " + "(at its cap, or because they had stopped settling - the log says which), so " + "everything read off that merge is read off an unsettled state", !result.scaling_converged && !result.search_scaling_converged ? "the space-group search merge and the final merge" : !result.scaling_converged ? "the final merge"