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;