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); } }