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") {