From 789e064e6abc08497f9502d5e7c8c8881f769060 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 24 Oct 2026 15:51:06 +0200 Subject: [PATCH] Stills scaling: a model error in the per-crystal scale fit, and the rocking width measured per window group Two changes to the stills scale model; rotation is untouched (p.hkl md5 identical on two rotation sets). The per-crystal scale fits (ScaleOnTheFly and the partiality post-refinement's G and tilt) measured residuals in counting sigmas. The strong reflections of a still scatter about the model by tens of percent (robust sd of ln(I/) 0.45-0.8 by shell on one run), so with a Cauchy loss at 3 counting sigma every strong reflection was an outlier and the scale was set by the weak ones: re-fitting each frame's scale against the final merge moved it by a factor 0.45-1.5 (5-95%). The fit's variance is now the expected one at the model, background + counting + (0.3 * model)^2. The coefficient is not delicate: 0.1-1.0 fit alike in an offline replay; what matters is that the term is there. The run's rocking width came from the intensity-weighted moment of the excitation error over all frames pooled. Frames predict through windows of their own that differ by an order of magnitude, and the widest windows dominated both moments: on a run whose windows are mostly narrower than the curve at high resolution the measurement read the window at every resolution, 2.8e-4 1/A with no angular term, several times too wide at low resolution (low-partiality observations there merged at 0.1-0.3 of their reflection's mean). The moment is now taken per group of frames of similar window (4 per shell), each against its own window, and the shell takes the median of the groups whose curve is narrower than their window. On that run the width becomes 6.7e-5 1/A + 0.018 deg angular; where the pooled measurement was already right (a SASE-bandwidth run), it stays where it was. 0-cycle REFMAC R-free against the deposited models, 9 still battery sets (6 used while developing, 3 held out), change from rc179: -0.0100, -0.0095, -0.0057, -0.0055, -0.0051, -0.0045, -0.0030, -0.0024, +0.0009 (held out: -0.0055, -0.0051, +0.0009). The width change alone accounts for most of the largest gain; the scale-fit variance for 0.000-0.003 per set. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu (cherry picked from commit ccc0e01f56827bbe0648ed7c5b8049912d80650f) --- image_analysis/scale_merge/ScaleOnTheFly.cpp | 72 +++++---- .../scale_merge/StillsPartialityRefine.cpp | 142 +++++++++++------- tests/MergeScaleTest.cpp | 43 ++++++ 3 files changed, 176 insertions(+), 81 deletions(-) diff --git a/image_analysis/scale_merge/ScaleOnTheFly.cpp b/image_analysis/scale_merge/ScaleOnTheFly.cpp index 823b42706..5cce4b17d 100644 --- a/image_analysis/scale_merge/ScaleOnTheFly.cpp +++ b/image_analysis/scale_merge/ScaleOnTheFly.cpp @@ -31,45 +31,57 @@ namespace { return 1.0 / x; } + // The scatter of a still about its scale model is not counting noise: the strong reflections of a frame + // scatter by tens of percent (partiality, orientation, the crystal itself), far beyond their counting + // sigma. Measured in counting sigmas alone, every strong reflection is an outlier to the robust fit and + // the scale is decided by the weak ones. So the fit's variance carries a model error proportional to + // the predicted intensity. Its size barely matters (0.1-1.0 all fit alike); its presence does. + constexpr double SCALE_MODEL_ERROR = 0.3; + // One reflection reduced to the 1-D scale fit: predicted intensity is G * coeff (coeff is constant - // while B is fixed), measured is Iobs, weighted by 1/sigma. + // while B is fixed), measured is Iobs; sigma2 its measured variance, var_bkg the background part of it. struct ScaleObs { double coeff; double Iobs; - double weight; + double sigma2; + double var_bkg; }; - // Robust per-image scale: minimise sum_i Cauchy_k( weight_i (G*coeff_i - Iobs_i) ) over G >= 0. The - // model is linear in G, so this M-estimate is a few reweighted-least-squares steps (each a closed-form - // weighted ratio) - the same objective the Ceres path solves, without a per-image problem/autodiff/ - // trust-region. Seeded from the plain weighted-LS solution; Cauchy weight is 1/(1 + (res/k)^2). - double SolveScaleIRLS(const std::vector &obs, double robust_k) { - auto weighted_scale = [&obs](auto robust_weight) { - double num = 0.0, den = 0.0; - for (const auto &o: obs) { - const double rw = robust_weight(o); - const double w2 = o.weight * o.weight; - num += rw * w2 * o.coeff * o.Iobs; - den += rw * w2 * o.coeff * o.coeff; - } - return den > 0.0 ? num / den : NAN; - }; + // The variance of an observation at the model m = G * coeff: background, the counting noise of m and + // the model error. + double ModelVariance(const ScaleObs &o, double G) { + const double m = G * o.coeff; + return o.var_bkg + std::max(0.0, m) + SCALE_MODEL_ERROR * SCALE_MODEL_ERROR * m * m; + } - double G = weighted_scale([](const ScaleObs &) { return 1.0; }); - if (!std::isfinite(G)) + // Robust per-image scale: minimise sum_i Cauchy_k( (G*coeff_i - Iobs_i) / sqrt(var_i(G)) ) over G >= 0. + // The model is linear in G, so this M-estimate is a few reweighted-least-squares steps (each a + // closed-form weighted ratio), the variances re-evaluated at each step's G. Seeded from the plain + // least-squares solution weighted by the measured variances; Cauchy weight is 1/(1 + (res/k)^2). + double SolveScaleIRLS(const std::vector &obs, double robust_k) { + double num = 0.0, den = 0.0; + for (const auto &o: obs) { + num += o.coeff * o.Iobs / o.sigma2; + den += o.coeff * o.coeff / o.sigma2; + } + if (!(den > 0.0) || !std::isfinite(num / den)) return 1.0; - G = std::max(0.0, G); + double G = std::max(0.0, num / den); const double k2 = robust_k * robust_k; for (int iter = 0; iter < 30; ++iter) { - const double G_prev = G; - const double G_next = weighted_scale([&](const ScaleObs &o) { - const double res = o.weight * (G * o.coeff - o.Iobs); - return 1.0 / (1.0 + res * res / k2); - }); - if (!std::isfinite(G_next)) + num = den = 0.0; + for (const auto &o: obs) { + const double var = ModelVariance(o, G); + const double res = G * o.coeff - o.Iobs; + const double w = 1.0 / (var + res * res / k2); // Cauchy weight / var + num += w * o.coeff * o.Iobs; + den += w * o.coeff * o.coeff; + } + if (!(den > 0.0) || !std::isfinite(num / den)) break; - G = std::max(0.0, G_next); + const double G_prev = G; + G = std::max(0.0, num / den); if (std::abs(G - G_prev) <= 1e-7 * std::max(G, 1.0)) break; } @@ -118,7 +130,11 @@ void ScaleOnTheFly::Scale(IntegrationOutcome &integration_outcome) const { if (it == reference_data.end()) continue; const double coeff = r.partiality * SafeInv(r.prescaling_corr * r.qe_corr * r.flight_corr, 1.0) * it->second; - obs.push_back({coeff, static_cast(r.I), SafeInv(r.sigma, 1.0)}); + const double sigma2 = static_cast(r.sigma) * r.sigma; + if (!(sigma2 > 0.0) || !std::isfinite(sigma2)) + continue; + const double var_bkg = std::isfinite(r.var_bkg) && r.var_bkg > 0.0f ? r.var_bkg : sigma2; + obs.push_back({coeff, static_cast(r.I), sigma2, var_bkg}); } if (obs.size() < MIN_REFLECTIONS) { diff --git a/image_analysis/scale_merge/StillsPartialityRefine.cpp b/image_analysis/scale_merge/StillsPartialityRefine.cpp index 401c8882d..7919ce6be 100644 --- a/image_analysis/scale_merge/StillsPartialityRefine.cpp +++ b/image_analysis/scale_merge/StillsPartialityRefine.cpp @@ -22,6 +22,9 @@ namespace { constexpr size_t MIN_FIT_REFLECTIONS = 20; constexpr double kRadToDeg = 180.0 / 3.14159265358979323846; + // Model error of the per-crystal scale fit, as a fraction of the predicted intensity - ScaleOnTheFly's. + constexpr double SCALE_MODEL_ERROR = 0.3; + double SafeInv(double x, double fallback) { if (!std::isfinite(x) || x == 0.0) return fallback; @@ -30,16 +33,24 @@ namespace { // One accepted reflection reduced to the physical partiality fit: the base reciprocal vector q (crystal // frame at the stored per-image orientation), the reference full intensity, the measured intensity, the - // reciprocal of the prescaling correction (the factor that turns a true intensity back into counts) - // and the weight. dist_ewald / partiality are recomputed from q under the refined tilt. + // reciprocal of the prescaling correction (the factor that turns a true intensity back into counts), + // the measured variance and its background part, and the weight the tilt is refined with. + // dist_ewald / partiality are recomputed from q under the refined tilt. struct FitObs { double qx, qy, qz; double Iref; double Iobs; double inv_corr; // 1 / (prescaling_corr * qe_corr * flight_corr) - double weight; // 1 / sigma + double sigma2; + double var_bkg; + double weight; // 1 / sqrt(variance at the current model), set with each G }; + // The variance of an observation at the model m: background, the counting noise of m and the model error. + double ModelVariance(const FitObs &o, double m) { + return o.var_bkg + std::max(0.0, m) + SCALE_MODEL_ERROR * SCALE_MODEL_ERROR * m * m; + } + double VecLen(double x, double y, double z) { return std::sqrt(x * x + y * y + z * z); } // Analytic partiality for a reflection whose base reciprocal vector is q, tilted by (psi_x, psi_y). @@ -72,9 +83,8 @@ namespace { // their plain moment of de is u = W^2 / 3, and the intensity-weighted moment m of a Gaussian of width // s cut at +-W is u * g(a), a = W / s, with g(a) = 3/a^2 (1 - 2 a phi(a) / erf(a / sqrt 2)): 1 for a // curve much wider than the window, s^2 / u for a much narrower one. Solves g(a) = m / u for s^2 by - // bisection (g falls with a). A curve wider than the window cannot be measured through it, so it is - // taken as wide as the window (a = 1). - double UntruncatedVariance(double m, double u) { + // bisection (g falls with a). A curve wider than the window cannot be measured through it: empty. + std::optional UntruncatedVariance(double m, double u) { auto g = [](double a) { const double phi = std::exp(-0.5 * a * a) / std::sqrt(2.0 * 3.14159265358979323846); return 3.0 / (a * a) * (1.0 - 2.0 * a * phi / std::erf(a / std::sqrt(2.0))); @@ -82,7 +92,7 @@ namespace { const double t = m / u; double lo = 1.0, hi = 50.0; if (!(t < g(lo))) - return 3.0 * u; + return std::nullopt; for (int i = 0; i < 60; ++i) { const double mid = 0.5 * (lo + hi); (g(mid) > t ? lo : hi) = mid; @@ -91,36 +101,32 @@ namespace { return 3.0 * u / (a * a); } - // Robust per-crystal scale G (linear in G given the model coefficients), identical objective to - // ScaleOnTheFly::SolveScaleIRLS: minimise sum Cauchy_k( w (G*coeff - Iobs) ) over G >= 0. - double SolveScaleIRLS(const std::vector &coeff, const std::vector &Iobs, - const std::vector &weight, double robust_k) { - auto weighted_scale = [&](auto robust_weight) { - double num = 0.0, den = 0.0; - for (size_t i = 0; i < coeff.size(); ++i) { - const double rw = robust_weight(i); - const double w2 = weight[i] * weight[i]; - num += rw * w2 * coeff[i] * Iobs[i]; - den += rw * w2 * coeff[i] * coeff[i]; - } - return den > 0.0 ? num / den : NAN; - }; - - double G = weighted_scale([](size_t) { return 1.0; }); - if (!std::isfinite(G)) + // Robust per-crystal scale G (linear in G given the model coefficients), the objective of + // ScaleOnTheFly::SolveScaleIRLS: minimise sum Cauchy_k( (G*coeff - Iobs) / sqrt(var(G*coeff)) ) over G >= 0, + // the variances re-evaluated at each step's G. + double SolveScaleIRLS(const std::vector &obs, const std::vector &coeff, double robust_k) { + double num = 0.0, den = 0.0; + for (size_t i = 0; i < obs.size(); ++i) { + num += coeff[i] * obs[i].Iobs / obs[i].sigma2; + den += coeff[i] * coeff[i] / obs[i].sigma2; + } + if (!(den > 0.0) || !std::isfinite(num / den)) return 1.0; - G = std::max(0.0, G); + double G = std::max(0.0, num / den); const double k2 = robust_k * robust_k; for (int iter = 0; iter < 30; ++iter) { - const double G_prev = G; - const double G_next = weighted_scale([&](size_t i) { - const double res = weight[i] * (G * coeff[i] - Iobs[i]); - return 1.0 / (1.0 + res * res / k2); - }); - if (!std::isfinite(G_next)) + num = den = 0.0; + for (size_t i = 0; i < obs.size(); ++i) { + const double res = G * coeff[i] - obs[i].Iobs; + const double w = 1.0 / (ModelVariance(obs[i], G * coeff[i]) + res * res / k2); + num += w * coeff[i] * obs[i].Iobs; + den += w * coeff[i] * coeff[i]; + } + if (!(den > 0.0) || !std::isfinite(num / den)) break; - G = std::max(0.0, G_next); + const double G_prev = G; + G = std::max(0.0, num / den); if (std::abs(G - G_prev) <= 1e-7 * std::max(G, 1.0)) break; } @@ -284,6 +290,8 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, .Iref = it->second, .Iobs = static_cast(r.I), .inv_corr = SafeInv(r.prescaling_corr * r.qe_corr * r.flight_corr, 1.0), + .sigma2 = static_cast(r.sigma) * r.sigma, + .var_bkg = std::isfinite(r.var_bkg) && r.var_bkg > 0.0f ? r.var_bkg : static_cast(r.sigma) * r.sigma, .weight = SafeInv(r.sigma, 1.0), }); @@ -347,18 +355,18 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, const int inner = refine_tilt ? settings_.inner_iterations : 1; for (int it = 0; it < inner; ++it) { - // (1) Solve G given the current partialities. - std::vector coeff(obs.size()), Iobs(obs.size()), weight(obs.size()); + // (1) Solve G given the current partialities, and weight the tilt fit by the variances at that model. + std::vector coeff(obs.size()); for (size_t j = 0; j < obs.size(); ++j) { const double p = ComputeP(obs[j].qx, obs[j].qy, obs[j].qz, psi[0], psi[1], S0.x, S0.y, S0.z, inv_lambda, gamma0, gamma_e, bw); coeff[j] = p * obs[j].inv_corr * obs[j].Iref; - Iobs[j] = obs[j].Iobs; - weight[j] = obs[j].weight; } - G = SolveScaleIRLS(coeff, Iobs, weight, settings_.robust_k); + G = SolveScaleIRLS(obs, coeff, settings_.robust_k); if (!(G > 0.0) || !std::isfinite(G)) return 0.0; + for (size_t j = 0; j < obs.size(); ++j) + obs[j].weight = 1.0 / std::sqrt(ModelVariance(obs[j], G * coeff[j])); if (!refine_tilt) break; @@ -433,8 +441,17 @@ std::optional StillsPartialityRefine::MeasureR // instead the window's own (W^2 / 3). Measured in shells of d*^2 over the whole run, where the // intensities of thousands of crystals average out, and fitted as s^2 = gamma0^2 + gamma_e^2 d*^2 + // (bandwidth_sigma lambda d*^2 / 2)^2. + // + // Each frame predicts through a window of its own (twice its profile radius), and these differ by an + // order of magnitude within a run. Pooled, the frames with the widest windows decide both moments, and + // the correction then reads their window: on a run whose windows are mostly narrower than the curve at + // high resolution, the pooled measurement came out as the window at every resolution, a width constant + // in d* that is several times too wide at low resolution. So within a shell the frames are taken in + // WINDOW_GROUPS groups of similar window, each measured against its own, and the shell's width is the + // median of the groups whose curve is narrower than their window; a shell where no group can see the + // curve is left out. const auto &indexing = experiment_.GetIndexingSettings(); - struct Point { double x, de2, w, var_w; }; + struct Point { double x, de2, w, var_w, window; }; std::vector pts; for (const auto &o: outcomes) { if (!o.image_scale_g || !(*o.image_scale_g > 0.0)) @@ -442,6 +459,10 @@ std::optional StillsPartialityRefine::MeasureR if (reference_cell && !o.latt.GetUnitCell().is_close(*reference_cell, indexing.GetUnitCellDistTolerance(), indexing.GetUnitCellAngleTolerance_deg())) continue; + float window = 0.0f; // the largest excitation error the frame predicted + for (const auto &r: o.reflections) + if (std::isfinite(r.dist_ewald)) + window = std::max(window, std::fabs(r.dist_ewald)); for (const auto &r: o.reflections) { if (r.on_ice_ring || !AcceptReflection(r, d_min_limit_, d_max_limit_)) continue; @@ -450,38 +471,53 @@ std::optional StillsPartialityRefine::MeasureR if (!std::isfinite(w) || !std::isfinite(sigma_w) || !std::isfinite(r.dist_ewald) || !(r.d > 0.0f)) continue; pts.push_back({1.0 / (static_cast(r.d) * r.d), - static_cast(r.dist_ewald) * r.dist_ewald, w, sigma_w * sigma_w}); + static_cast(r.dist_ewald) * r.dist_ewald, w, sigma_w * sigma_w, window}); } } constexpr size_t N_SHELLS = 8; + constexpr size_t WINDOW_GROUPS = 4; if (pts.size() < 100 * N_SHELLS) return std::nullopt; std::sort(pts.begin(), pts.end(), [](const Point &a, const Point &b) { return a.x < b.x; }); Logger logger("StillsPartialityRefine"); - // A shell's moment is no better than the summed intensity it is weighted by, so a shell is used only - // where that sum is measured to 10 % (sum I / sigma(sum I) >= 10): the shells beyond the diffraction - // limit, where the sum is noise and can even be negative, say nothing about the width. + // A group's moment is no better than the summed intensity it is weighted by, so a group is used only + // where that sum is measured to 10 % (sum I / sigma(sum I) >= 10): beyond the diffraction limit the + // sum is noise and can even be negative, and says nothing about the width. constexpr double MIN_SHELL_SIGNAL = 10.0; std::vector sx, sy; for (size_t k = 0; k < N_SHELLS; ++k) { - double sw = 0.0, svar = 0.0, swde2 = 0.0, sde2 = 0.0, sumx = 0.0; const size_t lo = k * pts.size() / N_SHELLS, hi = (k + 1) * pts.size() / N_SHELLS; - for (size_t i = lo; i < hi; ++i) { - sw += pts[i].w; - svar += pts[i].var_w; - swde2 += pts[i].w * pts[i].de2; - sde2 += pts[i].de2; + double sumx = 0.0; + for (size_t i = lo; i < hi; ++i) sumx += pts[i].x; + std::sort(pts.begin() + lo, pts.begin() + hi, + [](const Point &a, const Point &b) { return a.window < b.window; }); + std::vector widths; + for (size_t g = 0; g < WINDOW_GROUPS; ++g) { + const size_t glo = lo + g * (hi - lo) / WINDOW_GROUPS, ghi = lo + (g + 1) * (hi - lo) / WINDOW_GROUPS; + double sw = 0.0, svar = 0.0, swde2 = 0.0, sde2 = 0.0; + for (size_t i = glo; i < ghi; ++i) { + sw += pts[i].w; + svar += pts[i].var_w; + swde2 += pts[i].w * pts[i].de2; + sde2 += pts[i].de2; + } + if (!(sw >= MIN_SHELL_SIGNAL * std::sqrt(svar)) || !(swde2 > 0.0)) + continue; + if (const auto v = UntruncatedVariance(swde2 / sw, sde2 / static_cast(ghi - glo))) + widths.push_back(*v); } - const double n = static_cast(hi - lo); - if (!(sw >= MIN_SHELL_SIGNAL * std::sqrt(svar)) || !(swde2 > 0.0)) + if (widths.empty()) continue; - const double x = sumx / n, y = UntruncatedVariance(swde2 / sw, sde2 / n); + std::sort(widths.begin(), widths.end()); + const size_t m = widths.size(); + const double x = sumx / static_cast(hi - lo); + const double y = m % 2 ? widths[m / 2] : 0.5 * (widths[m / 2 - 1] + widths[m / 2]); sx.push_back(x); sy.push_back(y); - logger.Debug("Stills rocking width at d = {:.2f} A: {:.2e} 1/A ({:.3f} deg), I/sigma {:.0f}", 1.0 / std::sqrt(x), - std::sqrt(y), std::sqrt(y / x) * kRadToDeg, sw / std::sqrt(svar)); + logger.Debug("Stills rocking width at d = {:.2f} A: {:.2e} 1/A ({:.3f} deg), {} of {} window groups", + 1.0 / std::sqrt(x), std::sqrt(y), std::sqrt(y / x) * kRadToDeg, m, WINDOW_GROUPS); } if (sx.size() < 3) return std::nullopt; diff --git a/tests/MergeScaleTest.cpp b/tests/MergeScaleTest.cpp index 04ba3bb63..0c84aea52 100644 --- a/tests/MergeScaleTest.cpp +++ b/tests/MergeScaleTest.cpp @@ -957,6 +957,49 @@ TEST_CASE("StillsPartialityRefine_MeasureRunWidth") { CHECK(width->bandwidth_sigma < 1e-3); } +// Frames predicted through a much wider window whose intensities do not follow the curve (a wrong lattice +// on the frame) must not decide the measured width: pooled with the others they hold most of the de^2, and +// the run's width then read as their window. Measured per group of similar window, they are a group whose +// curve looks wider than its window, which is left out. +TEST_CASE("StillsPartialityRefine_MeasureRunWidth_WideWindowFramesIgnored") { + constexpr double GAMMA0 = 1.0e-4, GAMMA_E = 0.035 * M_PI / 180.0; + DiffractionExperiment x; + x.SetSpaceGroup(gemmi::get_spacegroup_by_name("P 1")); + + std::mt19937 rng(7); + std::uniform_real_distribution unit(0.0, 1.0); + std::exponential_distribution wilson(1.0); + std::normal_distribution gauss(0.0, 1.0); + std::vector outcomes(400); + for (size_t f = 0; f < outcomes.size(); ++f) { + auto &o = outcomes[f]; + const bool wide = f % 4 == 3; + o.image_scale_g = 0.5 + unit(rng); + for (int j = 0; j < 200; ++j) { + const double d = 1.5 + 8.0 * unit(rng); + const double dstar = 1.0 / d; + const double s = std::sqrt(GAMMA0 * GAMMA0 + GAMMA_E * GAMMA_E * dstar * dstar); + const double window = (wide ? 20.0 : 4.0) * std::sqrt(GAMMA0 * GAMMA0 + GAMMA_E * GAMMA_E / (1.5 * 1.5)); + const double de = window * (2.0 * unit(rng) - 1.0); + const double p = wide ? 1.0 : std::exp(-0.5 * de * de / (s * s)); + Reflection r{}; + r.h = j; r.k = 1; r.l = 1; + r.d = static_cast(d); + r.dist_ewald = static_cast(std::fabs(de)); + r.I = static_cast(*o.image_scale_g * 1000.0 * wilson(rng) * p + 5.0 * gauss(rng)); + r.sigma = 5.0f; + r.prescaling_corr = 1.0f; + r.partiality = 1.0f; + r.image_scale_corr = static_cast(1.0 / *o.image_scale_g); + o.reflections.push_back(r); + } + } + const auto width = StillsPartialityRefine(x).MeasureRunWidth(outcomes, std::nullopt); + REQUIRE(width.has_value()); + CHECK(width->gamma0 == Catch::Approx(GAMMA0).epsilon(0.3)); + CHECK(width->gamma_e == Catch::Approx(GAMMA_E).epsilon(0.15)); +} + // Outlier rejection on the stills merge: one observation at ten times its reflection's intensity among // nineteen consistent ones is dropped, and nothing else is. TEST_CASE("MergeOnTheFly_RejectsAnOutlierInTheBand") {