From 419db6ba0e0d92b5cb6e43fa08a48cf1597da15c Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 21 Oct 2026 19:52:43 +0200 Subject: [PATCH] MergeOnTheFly: capture term and multiplicative outlier band for the stills merge (both opt-in) --capture-uncertainty now reaches the stills merge: each partial carries (c (1 - p) )^2 in its variance, the rotation combine's capture term per observation, in the merge weights and in the error-model samples alike. --reject-outliers on stills now uses the rotation merge's band (OutlierBand.h): exp(+-n s) about the reflection mean for the systematic part, s the weighted median |ln(I / median)| of the merge's own observations, instead of a symmetric +-n sigma band that reaches below zero when b is large (stills b 1.4-1.7). A median needs three observations. Neither is switched on for stills: measured with the run-measured partiality width before its window-truncation correction (dimple R-free, SwissFEL 4k / SACLA 243 / MicroMAX mono / MicroMAX pink): no options 0.2567 / 0.2497 / 0.3279 / 0.2477, --capture-uncertainty 1 0.2510 / 0.2470 / 0.3261 / 0.2527, --reject-outliers 6 0.2619 / 0.2487 / 0.3317 / 0.2630 - mixed, and the photon-counting pink set loses with both. Raising --min-partiality to 0.1 or 0.2 lost 0.008-0.019 on SwissFEL and mono. Defaults unchanged, so stills output is unchanged; rotation does not use MergeOnTheFly. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu --- image_analysis/scale_merge/Merge.cpp | 65 +++++++++++++++++++------ image_analysis/scale_merge/Merge.h | 8 +++- rugnux/rugnux_cli.cpp | 4 +- tests/MergeScaleTest.cpp | 72 ++++++++++++++++++++++++++++ 4 files changed, 131 insertions(+), 18 deletions(-) diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index e72daef80..16bfa0fcd 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -19,6 +19,7 @@ #include "HKLKey.h" #include "RfreeFlags.h" #include "FrenchWilson.h" +#include "OutlierBand.h" namespace { // Which half-set an observation prefers: a splitmix64 bit-mix of the image's stable index, so the @@ -69,6 +70,7 @@ MergeOnTheFly::MergeOnTheFly(const DiffractionExperiment &x) ? std::optional(scaling_settings.GetMinCCForImage()) : std::nullopt), min_partiality(scaling_settings.GetMinPartiality()), + capture_coeff(scaling_settings.GetCaptureUncertaintyCoeff()), generator(scaling_settings.GetMergeFriedel(), space_group), reject_outliers(scaling_settings.GetOutlierRejectNsigma() > 0.0), reject_nsigma(scaling_settings.GetOutlierRejectNsigma()) { @@ -100,18 +102,22 @@ size_t MergeOnTheFly::Contributions(const IntegrationOutcome &outcome, std::vect continue; auto hkl = generator(r); auto hkl_key = hkl.pack(); - sigma_corr = CorrectedSigma(I_corr, sigma_corr, r.image_scale_corr, r.var_bkg, hkl_key); + sigma_corr = CorrectedSigma(I_corr, sigma_corr, r.image_scale_corr, r.var_bkg, r.partiality, hkl_key); - // Robust outlier rejection: drop this observation if it sits more than - // reject_nsigma error-model sigmas from the reflection's median. Needs the active - // error model so sigma_corr reflects the real scatter (else the threshold is the - // bare counting sigma and would cull good partials). + // Robust outlier rejection: drop this observation if it sits outside the band about the + // reflection's median that the error model allows (OutlierBand.h). Needs the active error model + // so sigma_corr reflects the real scatter (else the threshold is the bare counting sigma and + // would cull good partials). if (reject_outliers && error_model_active) { const auto mit = reject_median_I.find(hkl_key); - if (mit != reject_median_I.end() && - std::fabs(I_corr - mit->second) > reject_nsigma * sigma_corr) { - ++rejected; - continue; + if (mit != reject_median_I.end()) { + const double I_b = error_model_mean_I.at(hkl_key); + if (OutsideOutlierBand(I_corr, mit->second, reject_nsigma, + static_cast(sigma_corr) * sigma_corr, I_b, + error_model_b * I_b, reject_up, reject_down)) { + ++rejected; + continue; + } } } @@ -188,7 +194,7 @@ void MergeOnTheFly::AddImages(const std::vector &outcomes) { } float MergeOnTheFly::CorrectedSigma(float I_corr, float sigma_corr, float image_scale_corr, - float var_bkg, + float var_bkg, float partiality, uint64_t hkl_key) const { if (!error_model_active) return sigma_corr; @@ -217,6 +223,8 @@ float MergeOnTheFly::CorrectedSigma(float I_corr, float sigma_corr, float image_ if (base > 0.0) a_var = base; } + const double capture = capture_coeff * (1.0 - std::min(1.0f, partiality)) * std::max(0.0, I_for_b); + a_var += capture * capture; const double v = error_model_a * a_var + (error_model_b * I_for_b) * (error_model_b * I_for_b); @@ -239,11 +247,13 @@ void MergeOnTheFly::RefineErrorModel(const std::vector &outc error_model_chi2 = 0.0; error_model_mean_I.clear(); reject_median_I.clear(); + reject_up = 0.0; + reject_down = 0.0; reject_count = 0; // --- 1. Collect accepted, scaled observations grouped by symmetry-equivalent hkl, // applying exactly the filters AddImage uses. --- - struct Obs { float I, sigma; }; + struct Obs { float I, sigma, p; }; struct KeyedObs { uint64_t key; Obs obs; }; std::unordered_map> groups; @@ -267,7 +277,7 @@ void MergeOnTheFly::RefineErrorModel(const std::vector &outc const float sigma_corr = r.sigma * r.image_scale_corr; if (!std::isfinite(I_corr) || !std::isfinite(sigma_corr) || sigma_corr <= 0.0f) continue; - out.push_back({generator(r).pack(), {I_corr, sigma_corr}}); + out.push_back({generator(r).pack(), {I_corr, sigma_corr, r.partiality}}); } }, [&](size_t, const std::vector &out) { @@ -281,6 +291,7 @@ void MergeOnTheFly::RefineErrorModel(const std::vector &outc // (b*I)^2 term uses the reflection mean, so the mean (not I_i) is the abscissa. --- struct Sample { double s2, I2, dev2; }; std::vector samples; + std::vector> log_spread; // |ln(I / median)|, weight for (const auto &[key, obs]: groups) { if (obs.size() < 2) @@ -296,8 +307,9 @@ void MergeOnTheFly::RefineErrorModel(const std::vector &outc const double mean = sum_wI / sum_w; error_model_mean_I[key] = static_cast(mean); // Robust centre for outlier rejection: the median intensity (resists the very - // outliers the inverse-variance mean is being protected from). Only when active. - if (reject_outliers) { + // outliers the inverse-variance mean is being protected from). Only when active, and only + // where there are three observations for a median to stand on. + if (reject_outliers && obs.size() >= 3) { std::vector iv; iv.reserve(obs.size()); for (const auto &o: obs) @@ -312,10 +324,33 @@ void MergeOnTheFly::RefineErrorModel(const std::vector &outc if (factor < 0.05) continue; const double resid = static_cast(o.I) - mean; - samples.push_back({static_cast(o.sigma) * o.sigma, I2, resid * resid / factor}); + const double capture = capture_coeff * (1.0 - std::min(1.0f, o.p)) * std::max(0.0, mean); + const double s2 = static_cast(o.sigma) * o.sigma + capture * capture; + samples.push_back({s2, I2, resid * resid / factor}); + const auto mit = reject_median_I.find(key); + if (mit != reject_median_I.end() && o.I > 0.0f && mit->second > 0.0f && mean > 0.0) + log_spread.emplace_back(std::fabs(std::log(o.I / mit->second)), I2 / s2); } } + // The spread of ln(I / median) the band's multiplicative part is drawn from: its weighted median + // absolute value, weighted by (/sigma)^2 so the observations whose ratio counting noise does not + // blur carry it - as the rotation merge measures it. + if (!log_spread.empty()) { + std::sort(log_spread.begin(), log_spread.end()); + double total = 0.0; + for (const auto &[v, w]: log_spread) total += w; + double run = 0.0; + size_t j = 0; + for (; j + 1 < log_spread.size(); ++j) { + run += log_spread[j].second; + if (run >= 0.5 * total) break; + } + const double s = 1.4826 * log_spread[j].first; + reject_up = std::exp(reject_nsigma * s) - 1.0; + reject_down = 1.0 - std::exp(-reject_nsigma * s); + } + // --- 3. Fit global dev2 = a*sigma^2 + b^2*^2. Bin by intensity (the per-observation // dev2 is chi-square-1 noisy) and take medians; weight the bins by 1/dev2^2 so it // is a *relative* fit - otherwise the strong bins (which fix b) swamp the weak diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 6d73040f0..6c8ca78a7 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -294,6 +294,9 @@ class MergeOnTheFly { // merge, the error model and MergeStats can never disagree about which images are in. bool filter_by_image_cc = false; double min_partiality = 0.02; + // The share of a partial the model could not see, (1 - p), carried in its variance as + // (capture_coeff * (1 - p) * )^2 - the rotation combine's capture term, per observation. + double capture_coeff = 0.0; // When set, ice-ring-flagged reflections are left out of this merge. Used for the P1 pass whose // merged intensities feed the space-group search and the error model - those model fits must not @@ -319,7 +322,7 @@ class MergeOnTheFly { // using the per-observation I_i instead would over-weight down-fluctuated points. std::unordered_map error_model_mean_I; [[nodiscard]] float CorrectedSigma(float I_corr, float sigma_corr, float image_scale_corr, - float var_bkg, + float var_bkg, float partiality, uint64_t hkl_key) const; // Optional per-observation outlier rejection: drop observations whose corrected @@ -330,6 +333,9 @@ class MergeOnTheFly { bool reject_outliers = false; double reject_nsigma = 6.0; std::unordered_map reject_median_I; + // The asymmetric band of OutlierBand.h: exp(+-n s) - 1 about the reflection mean, s the spread of + // ln(I / median) measured on this merge's observations. + double reject_up = 0.0, reject_down = 0.0; size_t reject_count = 0; bool Mask(const IntegrationOutcome &outcome) const; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index b2e641876..d37bf3deb 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -229,11 +229,11 @@ void print_usage() { std::cout << " --ice-min-score Ice-presence gate: measured ice score (1 = no ice) a run must reach before ANY ice handling is applied - the flagging and the exclusion from scaling (default: 1.5). The eleven hexagonal bands cover 16-26% of the unique reflections whether or not the crystal has ice, so handling ice on a clean crystal is a pure loss. 0 = no gate (always handle ice)" << std::endl; std::cout << " --ice-min-spot-ratio Second ice-presence channel: found spots on the hexagonal rings over the same q width of ice-free flanks beside them (1 = spots spread evenly). Ice in large crystallites diffracts as discrete spots and leaves the radial profile flat, so --ice-min-score alone is blind to it. Default 2.0; 0 disables this channel" << std::endl; std::cout << " --min-partiality Minimum partiality to accept reflection (default: 0.02)" << std::endl; - std::cout << " --capture-uncertainty rot3d: systematic sigma ~num*(1-captured_fraction)*I on under-captured fulls (default: 1.0 for rot3d, 0 otherwise)" << std::endl; + std::cout << " --capture-uncertainty Systematic sigma num*(1-captured fraction)*: on rot3d fulls, the share of the rocking curve not captured; on stills, the share of each partial its partiality does not see (default: 1.0 for rot3d, 0 otherwise)" << std::endl; std::cout << " --min-captured-fraction rot3d: drop a combined full whose rocking curve was captured below this fraction (edge-of-sweep truncated fulls) (default: 0.7 for rotation, 0 otherwise; 0 = off)" << std::endl; std::cout << " --mosaicity Diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed" << std::endl; std::cout << " --prediction-mosaicity Diagnostic: fix the rocking width the PREDICTION window opens to (deg), leaving the partiality on the per-image sigma_M. Prediction and partiality are one number by default, so a sigma_M that moves takes the integrated reflection population with it" << std::endl; - std::cout << " --reject-outliers Per-observation merge outlier rejection, N sigma from the per-reflection median (default: 6 for rot3d, XDS/DIALS-style; 0 = off)" << std::endl; + std::cout << " --reject-outliers Per-observation merge outlier rejection, N sigma from the per-reflection median, the systematic part of the band multiplicative (default: 6 for rot3d, XDS/DIALS-style; 0 = off)" << std::endl; std::cout << " --search-min-zeta De-novo space-group search only: also search a merge of just the observations whose Lorentz geometry |zeta| reaches this, and report both answers (default: 0.85 for rotation, 0 = single search). Reflections crossing the Ewald sphere near-tangentially are measured worst and can make a real symmetry operator look like a twin law; where the two searches disagree, the merge of ALL the observations decides" << std::endl; std::cout << " --min-image-cc Per-image CC limit in percent (default: no limit)" << std::endl; std::cout << " --scaling-iterations Cap on the per-frame scaling iterations; the loop stops earlier once the scales settle (default: 100)" << std::endl; diff --git a/tests/MergeScaleTest.cpp b/tests/MergeScaleTest.cpp index 46148bbdb..04ba3bb63 100644 --- a/tests/MergeScaleTest.cpp +++ b/tests/MergeScaleTest.cpp @@ -956,3 +956,75 @@ TEST_CASE("StillsPartialityRefine_MeasureRunWidth") { CHECK(width->gamma_e == Catch::Approx(GAMMA_E).epsilon(0.1)); CHECK(width->bandwidth_sigma < 1e-3); } + +// 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") { + DiffractionExperiment x; + x.SetSpaceGroup(gemmi::get_spacegroup_by_name("P 1")); + ScalingSettings s = x.GetScalingSettings(); + s.OutlierRejectNsigma(6.0); + x.ImportScalingSettings(s); + + std::mt19937 rng(3); + std::normal_distribution gauss(0.0f, 1.0f); + std::vector outcomes(20); + for (size_t i = 0; i < outcomes.size(); ++i) { + for (int h = 1; h <= 50; ++h) { + Reflection r{}; + r.h = h; r.k = 2; r.l = 3; + r.d = 2.0f; + const float I_true = 100.0f * static_cast(h); + r.I = I_true * (1.0f + 0.1f * gauss(rng)) + 3.0f * gauss(rng); + if (i == 7 && h == 25) + r.I = 10.0f * I_true; + r.sigma = 3.0f; + r.prescaling_corr = 1.0f; + r.partiality = 1.0f; + r.image_scale_corr = 1.0f; + outcomes[i].reflections.push_back(r); + } + } + MergeOnTheFly merge(x); + merge.RefineErrorModel(outcomes); + REQUIRE(merge.ErrorModelActive()); + merge.AddImages(outcomes); + CHECK(merge.RejectedCount() == 1); +} + +// The capture term widens the sigma of a partial by (c (1 - p) ): a partial merges at a lower weight +// than a full of the same intensity, and not at all differently with c = 0. +TEST_CASE("MergeOnTheFly_CaptureTermDownweightsPartials") { + auto merged_I = [](double c) { + DiffractionExperiment x; + x.SetSpaceGroup(gemmi::get_spacegroup_by_name("P 1")); + ScalingSettings s = x.GetScalingSettings(); + s.CaptureUncertaintyCoeff(c); + x.ImportScalingSettings(s); + std::mt19937 rng(9); + std::normal_distribution gauss(0.0f, 1.0f); + std::vector outcomes(40); + for (size_t i = 0; i < outcomes.size(); ++i) + for (int h = 1; h <= 40; ++h) { + Reflection r{}; + r.h = h; r.k = 1; r.l = 1; + r.d = 2.0f; + // Every other image records the reflection as a 30 % partial that reads 20 % high. + const bool partial = i % 2 == 1; + r.partiality = partial ? 0.3f : 1.0f; + r.image_scale_corr = 1.0f / r.partiality; + r.I = r.partiality * 100.0f * static_cast(h) * (partial ? 1.2f : 1.0f) * (1.0f + 0.05f * gauss(rng)); + r.sigma = 2.0f; + r.prescaling_corr = 1.0f; + outcomes[i].reflections.push_back(r); + } + MergeOnTheFly merge(x); + merge.RefineErrorModel(outcomes); + merge.AddImages(outcomes); + const auto m = merge.ExportReflections(); + return m.at(10).I / (100.0f * static_cast(m.at(10).h)); + }; + const float without = merged_I(0.0), with = merged_I(1.0); + CHECK(with < without); + CHECK(with < 1.08f); +}