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"