From cf7676f5d5825d8e086e26e272cf121d250eba81 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 25 Oct 2026 19:45:42 +0100 Subject: [PATCH] Stills merge: outlier rejection on by default (10 sigma), the band linear in the error model's sigma The stills merge had no working outlier rejection. --reject-outliers was 0 for stills, and when set its band took the systematic part as a factor exp(+-n s), s the spread of ln(I / median) - on stills a reflection's observations are partials and per-crystal scales, many near zero, so s is wide and the band let through observations hundreds of sigma out. On a high-multiplicity broad-bandpass run with a spotty powder contaminant, single grains in a weak reflection's box were that reflection's whole merged intensity (half-set means -2 vs 2284 on 141 observations); the shells holding the powder rings and the lowest shell merged at CC1/2 0.6-0.75 beside 0.99 elsewhere, at multiplicity ~400-1000. Stills now test the plain linear band, n sqrt(a sigma^2 + (b )^2) about the median - the error model's b is the spread as it is - at n = 10 by default (rotation keeps 6 and its multiplicative band). The stills merge statistics now report the rejected count (OBSERVATIONS_REJECTED was always 0 there). Measured (REFMAC 0-cycle R-free of an isomorphous room-temperature model, depositor's free set, scale mode on stored integrations, vs the same code with --reject-outliers 0): - broad-bandpass in-house stills, 50k frames: 0.318 -> 0.265; 20k: 0.274 -> 0.244; R_split 0.29 -> 0.06, CC1/2 0.79 -> 0.98. 0.1% of the observations are rejected. - monochromatic in-house stills, 30k frames: 0.216 -> 0.212; R_split 0.126 -> 0.087. - n = 6 gains ~0.003 more there but costs the public XFEL sets 0.003-0.006 on three of four; 20 keeps most of the gain. Still battery, 12 non-float sets, against rc180: R_split better on all 12; R-free change mean -0.0005 (range -0.006..+0.003, plus pal86 +0.009 that is within that set's run-to-run spread of ~0.007 - a rerun of rc180's behaviour gave 0.192 against the baseline's 0.185); one draw becomes a win (lcls242_s), nothing else changes class. Rotation is untouched (the stills band and default are gated on the experiment not being a rotation). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu (cherry picked from commit 3a7142801056f21c720956255f1207e0bf1a8f43) --- docs/CPU_DATA_ANALYSIS_DECISIONS.md | 2 +- docs/RUGNUX_ADVANCED.md | 2 +- image_analysis/scale_merge/Merge.cpp | 13 +++++++++- image_analysis/scale_merge/Merge.h | 3 ++- rugnux/RugnuxDefaults.cpp | 9 +++++-- rugnux/rugnux_cli.cpp | 2 +- tests/MergeScaleTest.cpp | 39 ++++++++++++++++++++++++++++ 7 files changed, 63 insertions(+), 7 deletions(-) diff --git a/docs/CPU_DATA_ANALYSIS_DECISIONS.md b/docs/CPU_DATA_ANALYSIS_DECISIONS.md index 08a820e98..46457e451 100644 --- a/docs/CPU_DATA_ANALYSIS_DECISIONS.md +++ b/docs/CPU_DATA_ANALYSIS_DECISIONS.md @@ -74,7 +74,7 @@ Beside them, on rotation data with `--finalist-ledger`, the run reports **twin-i ### 13.3 Outlier rejection -Merging applies an optional per-observation median-based $N\sigma$ cut (`--reject-outliers`, default 6σ for `rot3d`, off otherwise). On the rotation path the median is weighted by each full's counting variance taken at the reflection's **expected** intensity, as the merge weights are (§10.4): an observation's own Poisson variance falls with its own intensity, so own-variance weights hand the median to whichever equivalents came out low — on a strongly absorbing crystal, where some frames read equivalents ten to a hundred times down, the median sat on the absorbed ones and the test removed correct measurements far above it. The band about the median keeps the counting part linear, $\pm N\sigma_\text{counting}$, but takes the systematic part as a **factor**, $e^{\pm N s}$ about the expected intensity: the error model's systematic term is multiplicative — absorption, illuminated volume, a scale that is off — and as likely to be a factor $1/f$ as $f$, while a band symmetric in $I$ reaches far below the median and only a little above it once that term is large. $s$ is the spread of $\ln(I/\text{median})$ measured on the merge's own fulls (their weighted median of $|\ln(I/\text{median})|$, weighted by $(\langle I\rangle/\sigma_\text{counting})^2$ so the fulls whose ratio counting noise does not blur carry it), with the error model's $b$ as the fallback where it cannot be measured; to first order in $Ns$ the band is the linear one. The same $N\sigma$ cut is fed back into the error model: after an initial $a,b$ fit the parameters are re-fit once on the reflections that survive rejection (dropping any whose squared deviation exceeds $N^2\,[a\,\sigma^2 + (b\,\langle I\rangle)^2]$), so the calibrated errors describe the reflections that actually enter the merge rather than the pre-rejection pool. +Merging applies a per-observation median-based $N\sigma$ cut (`--reject-outliers`, default 6σ on rotation, 10σ on stills). On the rotation path the median is weighted by each full's counting variance taken at the reflection's **expected** intensity, as the merge weights are (§10.4): an observation's own Poisson variance falls with its own intensity, so own-variance weights hand the median to whichever equivalents came out low — on a strongly absorbing crystal, where some frames read equivalents ten to a hundred times down, the median sat on the absorbed ones and the test removed correct measurements far above it. The band about the median keeps the counting part linear, $\pm N\sigma_\text{counting}$, but takes the systematic part as a **factor**, $e^{\pm N s}$ about the expected intensity: the error model's systematic term is multiplicative — absorption, illuminated volume, a scale that is off — and as likely to be a factor $1/f$ as $f$, while a band symmetric in $I$ reaches far below the median and only a little above it once that term is large. $s$ is the spread of $\ln(I/\text{median})$ measured on the merge's own fulls (their weighted median of $|\ln(I/\text{median})|$, weighted by $(\langle I\rangle/\sigma_\text{counting})^2$ so the fulls whose ratio counting noise does not blur carry it), with the error model's $b$ as the fallback where it cannot be measured; to first order in $Ns$ the band is the linear one. The same $N\sigma$ cut is fed back into the error model: after an initial $a,b$ fit the parameters are re-fit once on the reflections that survive rejection (dropping any whose squared deviation exceeds $N^2\,[a\,\sigma^2 + (b\,\langle I\rangle)^2]$), so the calibrated errors describe the reflections that actually enter the merge rather than the pre-rejection pool. On **stills** the band is the plain linear one, $\pm N\sqrt{a\,\sigma^2 + (b\,\langle I\rangle)^2}$ about the median: a still's observations of one reflection are partials and per-crystal scales, many of them near zero, so the spread of $\ln(I/\text{median})$ is wide and a band drawn from it let through a powder grain or a neighbour's spot in the box at hundreds of $\sigma$ — on a high-multiplicity run a single such observation was a weak reflection's whole merged intensity, and they took the low-resolution shells' $\mathrm{CC}_{1/2}$ down to 0.6. The error model's $b$ is that spread as it is. $N$ is wider than on rotation because a still's partials scatter with heavier tails than the band's normal core; what has to go sits tens to hundreds of $\sigma$ out. The median needs three observations, so a reflection measured once or twice — typically one good observation and one artefact (a hot pixel, a zinger) after Friedel merging — is never tested by it, and an artefact that looks precise (thousands of counts, a small relative $\sigma$) can outweigh mates from weak frames and become the median itself. Every observation of the written merge is therefore also judged by **Wilson statistics**. Beyond 4 Å its $E^2 = I/(\varepsilon\,\langle I/\varepsilon\rangle_\mathrm{shell})$ is compared with a bound set by a budget of $\alpha = 0.01$ expected false rejections per dataset: with $N$ observations tested, $\ln(2N/\alpha)$ for acentrics and twice that for centrics, and the observation's lower confidence limit $I - z\sigma$ (with its own $\sigma$, and $z$ from the same budget) must exceed it, so noise in weak shells does not trigger it; a shell whose $\langle I\rangle$ is not itself established at that significance is not judged. The bound is widened by the tail scale of the dataset's own intensity distribution, measured on its well-measured observations (1 for a Wilson crystal; larger under a pseudo-translation or anisotropy). An improbable singleton is dropped; an improbable observation with company only when most of the reflection's other observations are probable and it disagrees with their mean beyond the errors — several large observations confirm each other. An observation whose signal disk lost pixels to the mask or to saturation is never counted as such a witness, since its profile estimate may be low: a strong reflection is not replaced by its clipped mate. The count is `OBSERVATIONS_REJECTED_WILSON=` in the report (included in `OBSERVATIONS_REJECTED=`), and the developer report lists each observation with its image and detector position. diff --git a/docs/RUGNUX_ADVANCED.md b/docs/RUGNUX_ADVANCED.md index cc4f8fc27..1ba91b3fb 100644 --- a/docs/RUGNUX_ADVANCED.md +++ b/docs/RUGNUX_ADVANCED.md @@ -449,7 +449,7 @@ Scaling and merging: | `--min-partiality ` | Minimum partiality to accept a reflection (default: 0.02) | | `--ice-min-score ` | Ice-presence gate: the measured per-run ice score (1 = no ice) a dataset must reach before **any** ice handling is applied — the flagging and the exclusion from scaling (default: 1.5; 0 = no gate). The fixed hexagonal bands cover 16–26 % of the unique reflections whether or not the crystal has ice, so handling ice on a clean crystal only costs completeness | | `--ice-min-spot-ratio ` | The 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) | -| `--reject-outliers ` | Per-observation outlier rejection, N σ from the per-reflection median (default: 6 for `rot3d`, off otherwise) | +| `--reject-outliers ` | Per-observation outlier rejection, N σ from the per-reflection median (default: 6 for rotation, 10 for stills; 0 = off) | | `--min-image-cc ` | Per-image CC limit, percent (default: no limit) | | `--search-min-zeta ` | De-novo space-group search only: also search a merge of just the observations whose Lorentz geometry \|ζ\| 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 — as it always has for the systematic absences | | `--mosaicity ` | Diagnostic: fix the scaling mosaicity (°) instead of using the per-image seed | diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index b07eb3593..b9454e546 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -72,7 +72,8 @@ MergeOnTheFly::MergeOnTheFly(const DiffractionExperiment &x) capture_coeff(scaling_settings.GetCaptureUncertaintyCoeff()), generator(scaling_settings.GetMergeFriedel(), space_group), reject_outliers(scaling_settings.GetOutlierRejectNsigma() > 0.0), - reject_nsigma(scaling_settings.GetOutlierRejectNsigma()) { + reject_nsigma(scaling_settings.GetOutlierRejectNsigma()), + stills(!x.IsRotationIndexing()) { } MergeOnTheFly &MergeOnTheFly::ReferenceCell(const std::optional &cell) { @@ -419,6 +420,15 @@ void MergeOnTheFly::RefineErrorModel(const std::vector &outc error_model_a = a; error_model_b = std::sqrt(b2); error_model_active = true; + // A still's observations of one reflection spread by their partiality and per-crystal scale, not by a + // multiplicative factor about the median: many read near zero, so the spread of ln(I / median) is + // wide and the band drawn from it lets through even a powder grain or a neighbour's spot at hundreds + // of sigma. The error model's b is that spread measured as it is, so stills test the plain band of + // n sigma about the median. + if (stills) { + reject_up = reject_nsigma * error_model_b; + reject_down = reject_nsigma * error_model_b; + } // Achieved goodness of fit: the median of the per-observation dev2/(a*sigma^2 + (b*)^2). That // ratio is chi-square(1)-distributed (median 0.4549) when the sigmas are correct, so normalize by @@ -842,6 +852,7 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m MergeStatistics out; out.shells.resize(n_shells); + out.n_observations_rejected = reject_count; for (int s = 0; s < n_shells; ++s) { const auto &sa = acc[s]; diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 3fe75bea5..cefedba66 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -333,9 +333,10 @@ class MergeOnTheFly { 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. + // ln(I / median) measured on this merge's observations; on stills n b, the linear band. double reject_up = 0.0, reject_down = 0.0; size_t reject_count = 0; + bool stills = false; bool Mask(const IntegrationOutcome &outcome) const; diff --git a/rugnux/RugnuxDefaults.cpp b/rugnux/RugnuxDefaults.cpp index 72ec86dae..8ef13bbb4 100644 --- a/rugnux/RugnuxDefaults.cpp +++ b/rugnux/RugnuxDefaults.cpp @@ -11,6 +11,10 @@ namespace { constexpr double SMOOTH_G_DEG = 5.0; // Per-observation outlier rejection for the rotation combine, in sigma. constexpr double REJECT_OUTLIERS_NSIGMA = 6.0; + // The same for the stills merge, set wider: a still's partials scatter with heavier tails than the + // error model's normal band, and what has to go - a powder grain or a neighbour's spot in the box - sits + // tens to hundreds of sigma out. + constexpr double REJECT_OUTLIERS_NSIGMA_STILLS = 10.0; } void ApplyRugnuxExperimentDefaults(DiffractionExperiment &experiment) { @@ -63,8 +67,9 @@ ScalingSettings RugnuxDefaultScalingSettings(bool rotation) { s.CaptureUncertaintyCoeff(rotation ? 1.0 : 0.0); // A single un-rejected artifact wrecks the intensity-correlation and second-moment statistics the - // space-group and twinning tests are built on. - s.OutlierRejectNsigma(rotation ? REJECT_OUTLIERS_NSIGMA : 0.0); + // space-group and twinning tests are built on. On stills one powder grain or neighbour's spot in a + // reflection's box can be its whole merged intensity, whatever its multiplicity. + s.OutlierRejectNsigma(rotation ? REJECT_OUTLIERS_NSIGMA : REJECT_OUTLIERS_NSIGMA_STILLS); // Run the de-novo space-group search a second time on a merge of only the well-measured // observations and keep whichever found more symmetry. It cannot lose symmetry, so an imperfect cut diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 095a9db93..2c27ec893 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -234,7 +234,7 @@ void print_usage() { 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, the systematic part of the band multiplicative (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 on rotation and linear on stills (default: 6 for rotation, XDS/DIALS-style, and 10 for stills; 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 d8fffe5f0..32301970c 100644 --- a/tests/MergeScaleTest.cpp +++ b/tests/MergeScaleTest.cpp @@ -1113,6 +1113,45 @@ TEST_CASE("MergeOnTheFly_RejectsAnOutlierInTheBand") { CHECK(merge.RejectedCount() == 1); } +// Stills: the observations of a reflection spread from near zero to its full intensity (partials), so +// the spread of ln(I / median) is wide. One powder grain in a weak reflection's box, at 200 times its +// intensity, must still be rejected - the band is the error model's, n sigma about the median - while +// none of the partials is. +TEST_CASE("MergeOnTheFly_StillsRejectsAGrainAmongPartials") { + DiffractionExperiment x; + x.SetSpaceGroup(gemmi::get_spacegroup_by_name("P 1")); + ScalingSettings s = x.GetScalingSettings(); + s.OutlierRejectNsigma(6.0); + x.ImportScalingSettings(s); + REQUIRE(!x.IsRotationIndexing()); + + std::mt19937 rng(5); + std::normal_distribution gauss(0.0f, 1.0f); + std::uniform_real_distribution partial(0.1f, 1.0f); + std::vector outcomes(40); + 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 * partial(rng) + 3.0f * gauss(rng); + if (i == 7 && h == 1) + r.I = 200.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") {