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