From ddc785590071bb96079eee5dd399e54d0cee24aa Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 20 Sep 2026 00:04:29 +0200 Subject: [PATCH] Rotation scaling: frame guards and delta-CC1/2 read the frames by what they carry Three follow-ups to the converging scale loop, all on sweeps with a stretch the crystal barely diffracted on: - The sweep ledger's scale channel (and with it delta-CC1/2's "normal frame" test) and the space-group search's scale floor measure a frame against the precision-weighted typical frame G_ref (TypicalFrameScale) instead of the run median: on a sweep that spent most of its turn out of the beam the median frame is itself a dead one, and every stretch then reads as typical. The two collapse guards stay on the median on purpose - a frame that collapsed toward zero and was not quite dropped has its fulls re-fitted with a scale of 1/G, a cubic mean is then theirs, and the floor read against it dropped three quarters of every live frame's fulls (measured). - delta-CC1/2's sigma-tau statistic enters each reflection with the information it carries (the same per-frame factor as the CC1/2 weight), so a dead stretch, whose scaled-up noise flooded the mean error variance of every mixed reflection and read as harm up to the 25% cap, now costs about nothing and is left to the ledger. - The scale loop's step test weighs each frame by its merge weight (its observations at its scale squared): a dead frame's scale is fitted on noise and wanders by orders of magnitude every iteration, carries nothing into the merge, is dropped after the loop, and must not hold the loop open. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_013nW6FNRP1bBJJ8pfHiByAT --- .../scale_merge/RotationScaleMerge.cpp | 93 +++++++++++++------ 1 file changed, 66 insertions(+), 27 deletions(-) diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index a01ca9958..2a6409224 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -1430,15 +1430,23 @@ RotationScaleMerge::ScalingLoopOutcome RotationScaleMerge::RunScalingLoop( } rescale(scaled_any, ratio); if (it > 0) { - double ss = 0.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]); - ss += d * d; + const double w = frame_obs_count[f] * g[f] * g[f]; + ss += w * d * d; + sw += w; ++n; } - out.step = n > 0 ? std::sqrt(ss / n) : 0.0; + out.step = sw > 0.0 ? std::sqrt(ss / sw) : 0.0; logger.Debug("Per-frame scaling ({}): iteration {}, rms |dlogG| {:.2e} over {} frames, " "gauge {:.4f}", what, it + 1, out.step, n, g_ref); if (out.step < SCALING_TOLERANCE) { @@ -1959,21 +1967,30 @@ void RotationScaleMerge::MeasureSweepQuality(const std::vector &partial return; // Run gauges. Only the ratio to the run matters - the absolute scale is degenerate with the merge's. - std::vector g_fitted, cc_fitted; + // The scale's gauge is the precision-weighted typical frame (TypicalFrameScale), not the median + // frame: on a sweep that spent most of its turn out of the beam the median frame IS a dead one, + // and every stretch then reads as typical (measured: a run whose good frames sat at 1.3-1.7 and + // whose median was 0.21 reported its dead stretch at "scale 0.24"). The CC channel keeps its + // median - a correlation is not a scale, and dead frames carry no CC to drag it with. + std::vector cc_fitted; + std::vector n_obs(n_frames, 0); + int n_scaled = 0; for (int f = 0; f < n_frames; ++f) { - if (partial_scaled[f] && std::isfinite(g_partial[f]) && g_partial[f] > 0.0) - g_fitted.push_back(g_partial[f]); + if (partial_scaled[f] && std::isfinite(g_partial[f]) && g_partial[f] > 0.0) { + n_obs[f] = frame_count[f]; + ++n_scaled; + } if (cc_n[f] >= MIN_REFLECTIONS_FOR_IMAGE_CC && std::isfinite(cc[f])) cc_fitted.push_back(cc[f]); } - if (static_cast(g_fitted.size()) < SWEEP_MIN_SCALED_FRAMES) + if (n_scaled < SWEEP_MIN_SCALED_FRAMES) return; auto median_of = [](std::vector &v) { const size_t mid = v.size() / 2; std::nth_element(v.begin(), v.begin() + mid, v.end()); return v[mid]; }; - const double g_typ = median_of(g_fitted); + const double g_typ = TypicalFrameScale(g_partial, n_obs); const double cc_typ = cc_fitted.empty() ? 0.0 : median_of(cc_fitted); // The two channels. The scale is the sample side alone: DivideOutIncidentFlux took the beam out of the @@ -2177,15 +2194,24 @@ namespace { // the variance of the unique reflections' mean intensities against the error variance of a half // dataset's mean, accumulated one reflection at a time. No half-set split, so no random numbers - // the same answer every run, and more precise than the split estimate. + // + // Each reflection enters with the information it carries: the precision its observations would + // have at the run's typical frame scale over the precision they have (the same per-frame factor + // as the CC1/2 weight, CCHalfFrameFactors), so a reflection measured only on frames the crystal + // barely diffracted on, whose scaled-up noise would otherwise flood the mean error variance of + // every mixed reflection, counts for what it knows. Unweighted, a dead stretch read as harm + // (-0.15 on a sweep 60% out of the beam) and was convicted up to the cap; weighted, it costs + // about nothing and is left to the ledger, which is what it is - a stretch that adds nothing. struct SigmaTau { - double n = 0, sum_mean = 0, sum_mean2 = 0, sum_eps = 0; - void Add(double mean, double eps) { - n += 1; sum_mean += mean; sum_mean2 += mean * mean; sum_eps += eps; + double n = 0, w = 0, sum_mean = 0, sum_mean2 = 0, sum_eps = 0; + void Add(double mean, double eps, double weight) { + n += 1; w += weight; sum_mean += weight * mean; sum_mean2 += weight * mean * mean; + sum_eps += weight * eps; } [[nodiscard]] double CC() const { - if (n < 4) return NAN; - const double var_y = (sum_mean2 - sum_mean * sum_mean / n) / (n - 1); - const double eps = sum_eps / n; + if (n < 4 || !(w > 0.0)) return NAN; + const double var_y = (sum_mean2 - sum_mean * sum_mean / w) / w * n / (n - 1); + const double eps = sum_eps / w; return var_y + eps / 2 > 0.0 ? (var_y - eps / 2) / (var_y + eps / 2) : NAN; } }; @@ -2293,12 +2319,17 @@ void RotationScaleMerge::MeasureBatchDeltaCCHalf(int n_groups) { // The observations, once: scaled intensity, ASU group (-1 = not in this merge, the same test the // merge itself applies), and the frame order that makes any frame range a contiguous span. - std::vector obs_I(n_full); + std::vector obs_I(n_full), obs_info(n_full); std::vector obs_group(n_full), order(n_full), fstart(n_frames + 1, 0); + std::vector n_obs_frame(n_frames, 0); + for (int i = 0; i < n_full; ++i) + if (UsableFull(fulls[i])) ++n_obs_frame[fulls[i].frame]; + const std::vector frame_factor = CCHalfFrameFactors(frame_scale_fulls, n_obs_frame); for (int i = 0; i < n_full; ++i) { const Obs &o = fulls[i]; obs_I[i] = o.I * o.corr; obs_group[i] = UsableFull(o) ? o.group : -1; + obs_info[i] = 1.0 / frame_factor[o.frame]; // the information the observation carries, <= 1 ++fstart[o.frame + 1]; } for (int f = 0; f < n_frames; ++f) fstart[f + 1] += fstart[f]; @@ -2310,16 +2341,18 @@ void RotationScaleMerge::MeasureBatchDeltaCCHalf(int n_groups) { // Per-reflection totals over the frames still in, and the per-reflection contribution of one frame // range - which is all a leave-one-out needs, because a count, a sum and a sum of squares subtract. std::vector tot_n(n_groups, 0), d_n(n_groups, 0), touched; - std::vector tot_s1(n_groups, 0.0), tot_s2(n_groups, 0.0), - d_s1(n_groups, 0.0), d_s2(n_groups, 0.0); + std::vector tot_s1(n_groups, 0.0), tot_s2(n_groups, 0.0), tot_w(n_groups, 0.0), + d_s1(n_groups, 0.0), d_s2(n_groups, 0.0), d_w(n_groups, 0.0); auto build_totals = [&]() { std::fill(tot_n.begin(), tot_n.end(), 0); std::fill(tot_s1.begin(), tot_s1.end(), 0.0); std::fill(tot_s2.begin(), tot_s2.end(), 0.0); + std::fill(tot_w.begin(), tot_w.end(), 0.0); for (int i = 0; i < n_full; ++i) { const int g = obs_group[i]; if (g < 0 || frame_rejected[fulls[i].frame]) continue; tot_n[g] += 1; tot_s1[g] += obs_I[i]; tot_s2[g] += obs_I[i] * obs_I[i]; + tot_w[g] += obs_info[i]; } }; @@ -2340,6 +2373,7 @@ void RotationScaleMerge::MeasureBatchDeltaCCHalf(int n_groups) { if (g < 0) continue; if (d_n[g] == 0) touched.push_back(g); d_n[g] += 1; d_s1[g] += obs_I[i]; d_s2[g] += obs_I[i] * obs_I[i]; + d_w[g] += obs_info[i]; } } // Only reflections that are still measurable WITHOUT the range count, which is eq. 5's own @@ -2351,11 +2385,11 @@ void RotationScaleMerge::MeasureBatchDeltaCCHalf(int n_groups) { if (tot_n[g] >= 2 && n_without >= 2) { double mean, eps; ReflectionMoments(tot_n[g], tot_s1[g], tot_s2[g], mean, eps); - with.Add(mean, eps); + with.Add(mean, eps, tot_w[g] / tot_n[g]); ReflectionMoments(n_without, tot_s1[g] - d_s1[g], tot_s2[g] - d_s2[g], mean, eps); - without.Add(mean, eps); + without.Add(mean, eps, (tot_w[g] - d_w[g]) / n_without); } - d_n[g] = 0; d_s1[g] = 0.0; d_s2[g] = 0.0; + d_n[g] = 0; d_s1[g] = 0.0; d_s2[g] = 0.0; d_w[g] = 0.0; } const double cc_with = with.CC(), cc_without = without.CC(); if (with.n < DELTA_CC_MIN_REFLECTIONS || !std::isfinite(cc_with) || !std::isfinite(cc_without)) @@ -3180,6 +3214,11 @@ 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. std::vector fitted; fitted.reserve(n); for (int f = 0; f < n; ++f) @@ -3256,6 +3295,7 @@ bool RotationScaleMerge::DropCollapsedFullScales(bool from_staging) { && o.corr != 1.0f) g_frame[o.frame] = 1.0 / static_cast(o.corr); + // The median, not the precision-weighted typical frame - see DropCollapsedScales. std::vector fitted; fitted.reserve(g_frame.size()); for (const double gf : g_frame) @@ -4834,14 +4874,13 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st // --- 2a'. Search pass only: drop the frames the crystal barely diffracted on, whose 1/G would // amplify their intensities into a P1 merge with no multiplicity to catch them. --- if (for_search) { - std::vector fitted; + std::vector n_obs(n_frames, 0); for (int f = 0; f < n_frames; ++f) if (frame_scaled_scratch[f] && std::isfinite(g_partial[f]) && g_partial[f] > 0.0) - fitted.push_back(g_partial[f]); - if (!fitted.empty()) { - const size_t mid = fitted.size() / 2; - std::nth_element(fitted.begin(), fitted.begin() + mid, fitted.end()); - const double floor_g = fitted[mid] * SEARCH_MIN_SCALE_RATIO; + n_obs[f] = frame_count[f]; + if (std::any_of(n_obs.begin(), n_obs.end(), [](int64_t c) { return c > 0; })) { + // Against the precision-weighted typical frame - see MeasureSweepQuality. + const double floor_g = TypicalFrameScale(g_partial, n_obs) * SEARCH_MIN_SCALE_RATIO; std::vector reject(n_frames, 0); int n_rejected = 0; for (int f = 0; f < n_frames; ++f) @@ -4865,7 +4904,7 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st if (reject[partials[i].frame]) partials[i].corr = 0.0f; }); logger.Info("Space-group search: ignoring {} of {} frames whose scale came out below " - "1/{:.0f} of the run median (the crystal barely diffracted there)", + "1/{:.0f} of the run's typical frame (the crystal barely diffracted there)", n_rejected, n_frames, 1.0 / SEARCH_MIN_SCALE_RATIO); } }