From a970e5223b3ce23d3faee03746f1977f85235fd2 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 21 Sep 2026 08:45:51 +0200 Subject: [PATCH] Rotation scaling: the smoothness restraint is released where the frames outvote it The per-frame scale is smoothed over the rotation inside both scaling loops, which is what stopped a fine-sliced sweep of a small cell from walking its own scales down fifteen decades. But the window is a prior about the dose and the crystal, and on a merge where one frame carries hundreds of observations the frame's own data outvote it. Smoothing them anyway does not remove the per-frame structure, it moves it into the error model: on a 360 degree merge at 14x multiplicity the scale the window averaged away is 2.2% rms, which is most of the 2.4% of unmodelled error the change added there, and it cost ISa 29.0 -> 23.8 and 43.3 -> 38.0 on two long high-multiplicity protein sweeps. So the loop over the fulls now ends with one free fit against the settled reference, shrunk frame by frame back toward the restrained scale by how much of the difference is real. Real and noise separate with nothing to tune, because a rocking event spans several frames and the fit's noise does not: whatever moved a frame's scale moved its neighbour's, so the lag-1 covariance of the difference is the systematic's variance and what is left of the total is the fit's. A sweep that measures its scales to a fraction of a percent keeps nearly all of what it measured; one sliced at a twentieth of its rocking curve, where a frame holds twenty-five fulls and the difference has no correlation with its neighbour at all, keeps none of it and is scaled exactly as before. Neither is told which it is. Only the fulls, and only off a settled loop. A partial's scale is the direction the collapse runs in and handing any of it back re-opens it - shrunk or not, it costs ISa 23.8 -> 14.8 on the same merge and takes a small-cell monoclinic set from 14.2 to 2.7 - while a full is a whole rocking event and cannot enter it; and a loop that walked until it was stopped has measured nothing to hand back. Measured against the current behaviour. Two long high-multiplicity protein sweeps: ISa 23.8 -> 27.0 and 38.0 -> 41.3, with R_meas and the anomalous signal moving with them (SigAno 0.929 -> 0.949, CC_anom 0.18 -> 0.22 on the first). A long-wavelength, high-multiplicity series collected at three goniometer settings - the case the restraint is right for, sliced twenty times finer than its rocking curve: ISa 13.7 -> 13.9, 14.2 -> 14.6, 12.5 -> 12.8, SigAno up on all three, and the anomalous peak heights at the deposited model's own sites unchanged to a tenth of a sigma. The fine-sliced, small-cell and powder sets the restraint was added for are unchanged to the last reported digit, as are the short-sweep tetragonal pair and every space group. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_013nW6FNRP1bBJJ8pfHiByAT --- .../scale_merge/RotationScaleMerge.cpp | 124 +++++++++++++++++- .../scale_merge/RotationScaleMerge.h | 8 ++ 2 files changed, 128 insertions(+), 4 deletions(-) diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index c2445a590..602d26fa5 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -41,6 +41,10 @@ namespace { // The per-frame scale is smoothed over at least this many rocking curves (see Run). constexpr int SMOOTH_G_MIN_ROCKING_EVENTS = 6; + // ShrinkToRestrained reads two variances off the frames of the sweep; on fewer frames than this + // neither of them means anything, and the restraint is kept whole. + constexpr size_t MIN_FRAMES_FOR_SHRINKAGE = 20; + // 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. // It matters because the scale enters as 1/G, so a G that collapses toward zero multiplies every @@ -1409,8 +1413,86 @@ void RotationScaleMerge::UpdateCorr(std::vector &obs, const std::vector &g_restrained, + const std::vector &fitted, + const std::vector &frame_obs_count, + std::vector &g) const { + // g comes in as a free fit and goes out shrunk toward g_restrained. d is what a frame's own + // observations say about its scale beyond its neighbourhood's, and it is two things added: the + // per-frame systematic the smoothing window averaged away, and the noise of the frame's own fit. + // + // They separate without a constant to choose, because the first is shared with the next frame + // and the second is not: a rocking event spans several frames, so whatever moved this frame's + // scale moved its neighbour's, while the fit's noise comes from disjoint observations. The + // lag-1 covariance of d is therefore the systematic's variance and what is left of var(d) is + // the fit's, and each frame keeps the fraction of its d that its own observations support - + // everything on a merge that measures the scale to a fraction of a percent, nothing on one + // where d is the fit talking to itself. + // + // Measured: on a 360 deg, 14x lysozyme merge the fulls' d is 2.2% rms with a lag-1 + // autocorrelation of 0.44, and keeping it is worth ISa 23.8 -> 27.0; on a 1440 deg sweep sliced + // at a twentieth of its rocking curve, where a frame holds 25 fulls, d is 7.0% rms with a lag-1 + // autocorrelation of 0.02 and is discarded. Neither sweep was told which it was. + const int n = static_cast(g.size()); + std::vector d(n, 0.0); + std::vector d2, cross, inv_n; + auto usable = [&](int f) { + return fitted[f] && frame_obs_count[f] > 0 && g_restrained[f] > 0.0 + && std::isfinite(g[f]) && g[f] > 0.0; + }; + for (int f = 0; f < n; ++f) + if (usable(f)) d[f] = std::log(g[f] / g_restrained[f]); + for (int f = 0; f < n; ++f) { + if (!usable(f)) continue; + d2.push_back(d[f] * d[f]); + inv_n.push_back(1.0 / frame_obs_count[f]); + if (f + 1 < n && usable(f + 1)) cross.push_back(d[f] * d[f + 1]); + } + // Too few frames to read either variance off: leave the restrained scale where it is. + if (d2.size() < MIN_FRAMES_FOR_SHRINKAGE || cross.size() < MIN_FRAMES_FOR_SHRINKAGE) { + for (int f = 0; f < n; ++f) + if (usable(f)) g[f] = g_restrained[f]; + return; + } + // var(d) from the median of d^2 against the median of a chi^2 with one degree of freedom, so a + // handful of frames the fit lost control of cannot set the scale of the whole sweep. + std::vector d2_copy = d2; + const double var_d = median_of(d2_copy) / 0.4549; + // The cross products are signed, and a median has no useful null for a product, so the mean of + // the middle nine tenths instead - the same purpose, the same robustness, no distributional + // assumption. + std::sort(cross.begin(), cross.end()); + const size_t lo = cross.size() / 20, hi = cross.size() - cross.size() / 20; + double sum_cross = 0.0; + for (size_t i = lo; i < hi; ++i) sum_cross += cross[i]; + const double cov_1 = hi > lo ? sum_cross / static_cast(hi - lo) : 0.0; + const double var_real = std::clamp(cov_1, 0.0, var_d); + const double var_fit = std::max(0.0, var_d - var_real); + // The fit's variance goes as 1/n_obs, so one number sets it for every frame. + double mean_inv_n = 0.0; + for (double v : inv_n) mean_inv_n += v; + mean_inv_n /= static_cast(inv_n.size()); + const double fit_var_coeff = mean_inv_n > 0.0 ? var_fit / mean_inv_n : 0.0; + double kept_sum = 0.0; + int kept_n = 0; + for (int f = 0; f < n; ++f) { + if (!usable(f)) continue; + const double var_f = fit_var_coeff / frame_obs_count[f]; + const double kept = var_real + var_f > 0.0 ? var_real / (var_real + var_f) : 0.0; + g[f] = g_restrained[f] * std::exp(kept * d[f]); + kept_sum += kept; + ++kept_n; + } + logger.Info("Per-frame scaling ({}): the scales carry {:.0f}% of what each frame's own " + "observations say beyond its neighbourhood's - that deviation is {:.2f}% rms, of " + "which {:.2f}% is shared with the next frame and the rest is the fit's own noise", + what, 100.0 * (kept_n ? kept_sum / kept_n : 0.0), 100.0 * std::sqrt(var_d), + 100.0 * std::sqrt(var_real)); +} + RotationScaleMerge::ScalingLoopOutcome RotationScaleMerge::RunScalingLoop( const char *what, std::vector &g, const std::vector &frame_obs_count, int window, + bool release_restraint, 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 @@ -1434,6 +1516,15 @@ RotationScaleMerge::ScalingLoopOutcome RotationScaleMerge::RunScalingLoop( // 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 restraint is a prior, though, not a measurement, and on a merge where a frame carries + // hundreds of observations the frame's own data outvote it. So a loop that is asked to + // (release_restraint) ends with ONE free fit against the settled reference, shrunk frame by + // frame toward the restrained scale by how much of the difference is real - see + // ShrinkToRestrained. Only the loop over the fulls is asked to: a partial's scale is the + // direction the collapse runs in, and handing any of it back re-opens it (measured: ISa 23.8 -> + // 22.1 on a 360 deg, 14x merge), while a full is a whole rocking event and cannot enter it. And + // only off a loop that settled, since a loop still walking has measured nothing to hand back. + // // 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 @@ -1499,6 +1590,27 @@ RotationScaleMerge::ScalingLoopOutcome RotationScaleMerge::RunScalingLoop( g_prev = g; } frame_scaled_scratch = scaled_any; + // Only off a settled loop: the deviation below is a measurement of the frames, and a loop that + // walked until it was stopped has not measured them. + if (release_restraint && window > 1 && out.converged) { + const std::vector g_restrained = g; + iterate(); // one free fit against the settled reference + for (int f = 0; f < n_frames; ++f) + scaled_any[f] |= frame_scaled_scratch[f]; + const std::vector g_free = g; + ShrinkToRestrained(what, g_restrained, scaled_any, frame_obs_count, g); + 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) { + if (!scaled_any[f] || !(std::isfinite(g[f]) && g[f] > 0.0)) { ratio[f] = 1.0; continue; } + ratio[f] = g_free[f]; + g[f] /= g_ref; + ratio[f] /= g[f]; + } + rescale(scaled_any, ratio); + frame_scaled_scratch = scaled_any; + } if (out.converged) logger.Info("Per-frame scaling ({}): settled after {} iterations (rms |dlogG| {:.1e})", what, out.iterations, out.step); @@ -4791,7 +4903,8 @@ 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, smooth_window, [&] { + partial_loop = RunScalingLoop("partials", g_partial, frame_count, smooth_window, + /*release_restraint=*/false, [&] { 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) @@ -4804,7 +4917,8 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st } #endif if (!scaled_on_gpu) { - partial_loop = RunScalingLoop("partials", g_partial, frame_count, smooth_window, [&] { + partial_loop = RunScalingLoop("partials", g_partial, frame_count, smooth_window, + /*release_restraint=*/false, [&] { 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); @@ -5026,7 +5140,8 @@ 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, smooth_window, [&] { + fulls_loop = RunScalingLoop("fulls", g_full, f_count, smooth_window, + /*release_restraint=*/true, [&] { gpu_->ScaleFulls(1, min_partiality); gpu_->GetG(g_dev.data(), scaled_dev.data()); for (int f = 0; f < n_frames; ++f) @@ -5086,7 +5201,8 @@ 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, smooth_window, [&] { + fulls_loop = RunScalingLoop("fulls", g_full, fulls_frame_count, smooth_window, + /*release_restraint=*/true, [&] { 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 a07017460..68fa1bf83 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -452,10 +452,18 @@ private: 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, int window, + bool release_restraint, const std::function &iterate, const std::function &, const std::vector &)> &rescale); + // Shrink a free per-frame scale back toward the restrained one, frame by frame, by how much of + // the difference between them is real (see ShrinkToRestrained in the .cpp). + void ShrinkToRestrained(const char *what, const std::vector &g_restrained, + const std::vector &fitted, + const std::vector &frame_obs_count, + std::vector &g) const; + // corr = prescaling_corr / (partiality * G[frame]); leaves corr unchanged for frames that could not be fit. void UpdateCorr(std::vector &obs, const std::vector &g, const std::vector &frame_scaled) const;