From 3e0285c4250defd4f37ab650784029b3f33b6b8f Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 21 Oct 2026 19:52:43 +0200 Subject: [PATCH] StillsPartialityRefine: the partiality width is the run's, measured on the crystals, not the prediction window's The stills partiality width came from a per-crystal least-squares fit of de^2 against d*^2 over the integrated reflections. Those fill the prediction window evenly, so the fit returned the window's own width (W^2/3, constant in the linear Ewald distance, gamma_e ~ 0): at low resolution several times the crystal's rocking width, which mis-corrected every partial there. The width is now measured once per run from what the reflections RECORD: the intensity-weighted second moment of the excitation error, sum I de^2 / sum I, which is the width of a Gaussian partiality whatever the window, in 8 shells of d*^2 over all crystals. Each shell is corrected for the window cutting the curve's tails (a curve wider than the window is taken as wide as the window), and the shells are fitted as s^2 = gamma0^2 + (gamma_e d*)^2 + (bw lambda d*^2 / 2)^2 with every term non-negative. The per-crystal fit stays as the fallback when the run is too small to measure (< 800 reflections). The prediction window is unchanged (2x the per-frame profile radius). Measured: on the SwissFEL lysozyme stills the width comes out ~0.035 deg angular at every resolution, the ~0.03-0.045 deg that the empirical rocking profile of the observations gives; the window-moment model was 3-4x wider at low resolution. dimple R-free (deposited free set, matrix subsets/cuts, PARTIALITY_MATRIX.md), shipped -> this: SwissFEL run0237 4k 0.2692 -> 0.2454 same, -C/-S given 0.2679 -> 0.2513 MicroMAX mono 4k 0.3341 -> 0.3189 MicroMAX pink 4k 0.2415 -> 0.2442 SACLA 243 subset 0.2576 -> 0.2454 SACLA 237 subset 0.2342 -> 0.2331 LCLS 222 subset 0.2344 -> 0.2330 PAL 86 subset 0.3499 -> 0.3629 (binned CCD, stencil-limited) For comparison --bandwidth 0.003 by hand: 0.2616 / - / 0.3108 / 0.2377 / 0.2460 / 0.2270 / 0.2404 / 0.3469. Stills only: rotation never runs StillsPartialityRefine. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu --- .../scale_merge/StillsPartialityRefine.cpp | 148 +++++++++++++++++- .../scale_merge/StillsPartialityRefine.h | 19 ++- tests/MergeScaleTest.cpp | 44 +++++- 3 files changed, 199 insertions(+), 12 deletions(-) diff --git a/image_analysis/scale_merge/StillsPartialityRefine.cpp b/image_analysis/scale_merge/StillsPartialityRefine.cpp index 21785a263..98664cef7 100644 --- a/image_analysis/scale_merge/StillsPartialityRefine.cpp +++ b/image_analysis/scale_merge/StillsPartialityRefine.cpp @@ -11,9 +11,11 @@ #include #include +#include #include #include "Merge.h" +#include "../../common/Logger.h" #include "../../common/ThreadAffinity.h" namespace { @@ -66,6 +68,29 @@ namespace { return std::exp(-0.5 * de * de / sig2); } + // The prediction window cuts the rocking curve. Reflections fill a window of half-width W evenly, so + // 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) { + 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))); + }; + const double t = m / u; + double lo = 1.0, hi = 50.0; + if (!(t < g(lo))) + return 3.0 * u; + for (int i = 0; i < 60; ++i) { + const double mid = 0.5 * (lo + hi); + (g(mid) > t ? lo : hi) = mid; + } + const double a = 0.5 * (lo + hi); + 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, @@ -221,7 +246,8 @@ StillsPartialityRefine::StillsPartialityRefine(const DiffractionExperiment &x) bandwidth_sigma_(x.GetBandwidthFWHM().value_or(0.0f) / 2.3548f) {} double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, - const std::map &reference) const { + const std::map &reference, + const std::optional &run_width) const { if (outcome.reflections.empty()) return 0.0; @@ -230,7 +256,6 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, const Coord Cstar = outcome.latt.Cstar(); const Coord S0 = outcome.geom.GetScatteringVector(); const double inv_lambda = 1.0 / outcome.geom.GetWavelength_A(); - const double bw = bandwidth_sigma_; auto base_q = [&](const Reflection &r) { return Astar * static_cast(r.h) + Bstar * static_cast(r.k) @@ -283,7 +308,8 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, if (obs.size() < MIN_FIT_REFLECTIONS || n_de == 0) return 0.0; - // Solve the 2x2 normal equations for de^2 = A + B d*^2. A degenerate spread in d* (all reflections in + // Without the run's measured width (too few reflections in the run to measure it), the width comes + // from this crystal's own reflections. Solve the 2x2 normal equations for de^2 = A + B d*^2. A degenerate spread in d* (all reflections in // one shell) leaves B undetermined, so fall back to the pure angular width there; a negative fitted // component is unphysical and is clamped to zero, which reduces to the previous model. // @@ -292,8 +318,8 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, // outside any measured range), i.e. a width constant in the LINEAR Ewald distance. That is // structural, not a fluke: prediction accepts reflections on a fixed linear |dist_ewald| cutoff, so // the accepted population's de^2 is flat in d*^2 by construction and the slope is genuinely ~0. The - // truncated population cannot constrain an angular term; the resolution-independent one is what the - // data actually support. + // truncated population cannot constrain an angular term, and its width is the window's, not the + // crystal's - which is why the run's measured width is used wherever there is one. const double det = m_n * m_xx - m_x * m_x; double A = 0.0, B = 0.0; if (std::fabs(det) > 1e-30) { @@ -302,9 +328,15 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, } else { B = m_x > 0.0 ? m_y / m_x : 0.0; } - const double gamma0 = std::sqrt(std::max(0.0, A)); + double gamma0 = std::sqrt(std::max(0.0, A)); const double gamma_e_fit = std::max(std::sqrt(std::max(0.0, B)), 1e-9); - const double gamma_e = settings_.gamma_e > 0.0 ? settings_.gamma_e : gamma_e_fit; + double gamma_e = settings_.gamma_e > 0.0 ? settings_.gamma_e : gamma_e_fit; + double bw = bandwidth_sigma_; + if (run_width) { + gamma0 = run_width->gamma0; + gamma_e = std::max(run_width->gamma_e, 1e-9); + bw = run_width->bandwidth_sigma; + } double psi[2] = {0.0, 0.0}; double G = 1.0; @@ -391,6 +423,105 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, return tilt_deg; } +std::optional StillsPartialityRefine::MeasureRunWidth( + const std::vector &outcomes, const std::optional &reference_cell) const { + // The rocking width measured on the crystals rather than on the prediction window. A lattice point + // sits anywhere across the window with equal likelihood, while what it records is I_full * p(de), so + // the intensity-weighted second moment of the excitation error, sum I de^2 / sum I, is the width s^2 + // of a Gaussian p - corrected for the window cutting the curve's tails (UntruncatedVariance). The + // plain moment of de over the predicted reflections, which the per-crystal fit in RefineOne uses, is + // 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. + const auto &indexing = experiment_.GetIndexingSettings(); + struct Point { double x, de2, w; }; + std::vector pts; + for (const auto &o: outcomes) { + if (!o.image_scale_g || !(*o.image_scale_g > 0.0)) + continue; + if (reference_cell && !o.latt.GetUnitCell().is_close(*reference_cell, indexing.GetUnitCellDistTolerance(), + indexing.GetUnitCellAngleTolerance_deg())) + continue; + for (const auto &r: o.reflections) { + if (r.on_ice_ring || !AcceptReflection(r, d_min_limit_, d_max_limit_)) + continue; + const double w = r.I * r.prescaling_corr * r.qe_corr * r.flight_corr / *o.image_scale_g; + if (!std::isfinite(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}); + } + } + constexpr size_t N_SHELLS = 8; + 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"); + std::array sx{}, sy{}; + for (size_t k = 0; k < N_SHELLS; ++k) { + double sw = 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; + swde2 += pts[i].w * pts[i].de2; + sde2 += pts[i].de2; + sumx += pts[i].x; + } + if (!(sw > 0.0)) + return std::nullopt; + const double n = static_cast(hi - lo); + sx[k] = sumx / n; + sy[k] = UntruncatedVariance(swde2 / sw, sde2 / n); + logger.Debug("Stills rocking width at d = {:.2f} A: {:.2e} 1/A ({:.3f} deg)", 1.0 / std::sqrt(sx[k]), + std::sqrt(sy[k]), std::sqrt(sy[k] / sx[k]) * kRadToDeg); + } + + // Least squares for s^2 = A + B d*^2 + C d*^4 over the shells. Every term is a variance, so a negative + // one is unphysical: of the fits on each subset of the three terms, the best with none negative is kept. + std::array best{}; + double best_residual = INFINITY; + for (int mask = 1; mask < 8; ++mask) { + std::vector terms; + for (int t = 0; t < 3; ++t) + if (mask & (1 << t)) + terms.push_back(t); + const int m = static_cast(terms.size()); + Eigen::MatrixXd M = Eigen::MatrixXd::Zero(m, m); + Eigen::VectorXd v = Eigen::VectorXd::Zero(m); + for (size_t k = 0; k < N_SHELLS; ++k) + for (int i = 0; i < m; ++i) { + const double fi = std::pow(sx[k], terms[i]); + v(i) += fi * sy[k]; + for (int j = 0; j < m; ++j) + M(i, j) += fi * std::pow(sx[k], terms[j]); + } + const Eigen::VectorXd c = M.ldlt().solve(v); + if (!c.allFinite() || (c.array() < 0.0).any()) + continue; + std::array coef{}; + for (int i = 0; i < m; ++i) + coef[terms[i]] = c(i); + double residual = 0.0; + for (size_t k = 0; k < N_SHELLS; ++k) { + const double r = sy[k] - coef[0] - coef[1] * sx[k] - coef[2] * sx[k] * sx[k]; + residual += r * r; + } + if (residual < best_residual) { + best_residual = residual; + best = coef; + } + } + if (!std::isfinite(best_residual) || !(best[0] + best[1] + best[2] > 0.0)) + return std::nullopt; + // The d*^4 term is a bandwidth: its radial smear bw |q_z|, with |q_z| = lambda d*^2 / 2 on the sphere. + const RunWidth width{.gamma0 = std::sqrt(best[0]), .gamma_e = std::sqrt(best[1]), + .bandwidth_sigma = 2.0 * std::sqrt(best[2]) / experiment_.GetWavelength_A()}; + logger.Info("Stills rocking width measured: {:.2e} 1/A + {:.3f} deg angular + {:.4f} bandwidth (FWHM)", + width.gamma0, width.gamma_e * kRadToDeg, width.bandwidth_sigma * 2.3548); + return width; +} + double StillsPartialityRefine::Run(std::vector &outcomes, size_t nthreads, const std::optional &reference_cell) const { const auto &indexing = experiment_.GetIndexingSettings(); @@ -404,6 +535,7 @@ double StillsPartialityRefine::Run(std::vector &outcomes, si nthreads = std::max(1, nthreads); double last_mean_tilt = 0.0; + const auto run_width = MeasureRunWidth(outcomes, reference_cell); for (int outer = 0; outer < settings_.outer_iterations; ++outer) { // Reference full intensities from the current corrections. @@ -419,7 +551,7 @@ double StillsPartialityRefine::Run(std::vector &outcomes, si auto worker = [&]() { size_t i = next.fetch_add(1); while (i < outcomes.size()) { - const double t = merged_crystal(outcomes[i]) ? RefineOne(outcomes[i], reference) : 0.0; + const double t = merged_crystal(outcomes[i]) ? RefineOne(outcomes[i], reference, run_width) : 0.0; if (t > 0.0) { double prev = tilt_sum.load(); while (!tilt_sum.compare_exchange_weak(prev, prev + t)) {} diff --git a/image_analysis/scale_merge/StillsPartialityRefine.h b/image_analysis/scale_merge/StillsPartialityRefine.h index f74568b65..a3600f0ca 100644 --- a/image_analysis/scale_merge/StillsPartialityRefine.h +++ b/image_analysis/scale_merge/StillsPartialityRefine.h @@ -3,6 +3,7 @@ #pragma once +#include #include #include "../../common/DiffractionExperiment.h" @@ -26,7 +27,9 @@ // dist_ewald = | |q + S0| - 1/lambda |) - NO pixel re-integration. The loop is: merge -> per-crystal // refine tilt (G profiled out by the same robust IRLS ScaleOnTheFly uses) -> recompute p and // image_scale_corr -> re-merge, iterated a few times. Mutates each reflection's `partiality` and -// `image_scale_corr` in place; the existing MergeOnTheFly then consumes the improved corrections. +// `image_scale_corr` in place; the existing MergeOnTheFly then consumes the improved corrections. The +// width of the Gaussian partiality is the run's, measured once from the intensity-weighted excitation +// errors (MeasureRunWidth). class StillsPartialityRefine { public: struct Settings { @@ -37,7 +40,7 @@ public: double prior_sigma_deg = 0.02; // soft prior pulling dpsi toward 0 (0 = off); tames weak-data overfit, // inert on strong data (well-supported tilts overcome it) double robust_k = 3.0; // Cauchy loss scale (sigma units) - double gamma_e = 0.0; // angular width sigma(d*) = gamma_e*d* (0 = estimate per crystal from data) + double gamma_e = 0.0; // angular width sigma(d*) = gamma_e*d* (0 = measure from data) }; explicit StillsPartialityRefine(const DiffractionExperiment &x); @@ -48,6 +51,15 @@ public: double Run(std::vector &outcomes, size_t nthreads, const std::optional &reference_cell) const; + // The run's rocking width, sigma^2 = gamma0^2 + (gamma_e d*)^2 + (bandwidth_sigma |q_z|)^2 (1/A, rad, + // relative), measured from the intensity-weighted excitation errors; empty when the run has too + // few reflections to measure it. + struct RunWidth { + double gamma0, gamma_e, bandwidth_sigma; + }; + std::optional MeasureRunWidth(const std::vector &outcomes, + const std::optional &reference_cell) const; + private: const DiffractionExperiment experiment_; const Settings settings_{}; @@ -59,5 +71,6 @@ private: // Refine one crystal against the reference map; returns |dpsi| in degrees (0 if skipped). double RefineOne(IntegrationOutcome &outcome, - const std::map &reference) const; + const std::map &reference, + const std::optional &run_width) const; }; diff --git a/tests/MergeScaleTest.cpp b/tests/MergeScaleTest.cpp index 70f08b06b..46148bbdb 100644 --- a/tests/MergeScaleTest.cpp +++ b/tests/MergeScaleTest.cpp @@ -720,6 +720,7 @@ TEST_CASE("StillsPartialityRefine_RecoversTheTilt") { Reflection r{}; r.h = h; r.k = k; r.l = l; r.d = 1.0f / dstar; + r.dist_ewald = std::fabs((q + S0).Length() - 1.0f / WAVELENGTH); r.I = static_cast(G * p_true * full_I(h, k, l)); r.sigma = static_cast(1.0 + 0.01 * r.I); r.prescaling_corr = 1.0f; @@ -745,7 +746,7 @@ TEST_CASE("StillsPartialityRefine_RecoversTheTilt") { const double after = mean_abs_error(); CHECK(mean_tilt > 0.015); CHECK(mean_tilt < 0.06); - CHECK(after < 0.6 * before); // measured 0.54: the prior holds the larger tilts back + CHECK(after < 0.6 * before); // measured 0.41 (0.54 with the per-crystal window width): the prior holds the larger tilts back } namespace { @@ -914,3 +915,44 @@ TEST_CASE("MergeOnTheFly_RsplitKeepsEveryFrameWhole") { CHECK(sh.r_split == Catch::Approx(expected).epsilon(1e-5)); CHECK_FALSE(SerialShellTable(stats).empty()); } + +// The run's rocking width is measured on the crystals, not on the prediction window: reflections spread +// evenly across a window four times wider than the curve, each recording I_full * p(de), give back the +// curve's gamma0 and angular width, while the plain moment of de would return the window's. +TEST_CASE("StillsPartialityRefine_MeasureRunWidth") { + 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(5); + 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(200); + for (auto &o: outcomes) { + 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 = 4.0 * std::sqrt(GAMMA0 * GAMMA0 + GAMMA_E * GAMMA_E / (1.5 * 1.5)); + const double de = window * (2.0 * unit(rng) - 1.0); + 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) * std::exp(-0.5 * de * de / (s * s)) + + 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.1)); + CHECK(width->bandwidth_sigma < 1e-3); +}