From ac06b5c64f7cce2215024c16cd8058305ccdd791 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 10 Aug 2026 18:50:22 +0200 Subject: [PATCH] Merging: report the error model in XDS's convention rugnux fits sigma^2 = a*sigma0^2 + (b*)^2, so its `b` is a fraction of the intensity. XDS fits sigma^2 = a*(sigma0^2 + b*I^2) and prints ISa = 1/sqrt(a*b). The two `a` are the same number, but the two `b` are not - b_xds = b^2/a - so the pair rugnux printed could not be read against a CORRECT.LP, which is the only reason anyone looks at it. Convert at the report. The fit, the merge weights and both engines' variance expressions are untouched, so this is a re-expression and not a change: on a rotation dataset the merged intensities move strictly less between before and after than they do between two runs of the SAME binary (99.9% identical, max |dI/I| 9.1e-4 against the run-to-run control's 7.5e-3), with the same reflection set. The rotation path also printed the wrong ISa for the comparison it invites. What it calls ISa is the strong-reflection asymptote, a tier XDS has no equivalent of and which can only ever be the more optimistic of the two; XDS's ISa is the whole-range 1/sqrt(a*b), which in rugnux units is exactly 1/b. Print both, labelled. On a broadband rotation dataset that is 13.2 (whole range) and 15.6 (asymptote) against XDS's 21.18 - so the number previously compared was flattering rugnux by 2.4. A third, unrelated `b` lives in the space-group search: fitted with the sigma^2 coefficient held at 1, with gate constants calibrated in that convention, and a ratio bound does not survive the mapping (1.90 would have to become 3.61) while the absolute floor has no correct value at all, there being no `a`. It is now commented as such, since making the three consistent is the obvious wrong move. Also corrects three comments and two doc passages that still described a merged-sigma systematic floor deleted in 72efb75a8. Co-Authored-By: Claude Opus 5 (1M context) --- docs/CHANGELOG.md | 1 + docs/CPU_DATA_ANALYSIS.md | 6 ++++-- image_analysis/scale_merge/Merge.h | 21 +++++++++++++++++++ .../scale_merge/RotationScaleMerge.cpp | 13 +++++++++--- .../scale_merge/RotationScaleMerge.h | 2 +- .../scale_merge/SearchSpaceGroup.cpp | 4 ++++ rugnux/Rugnux.cpp | 11 ++++++---- rugnux/Rugnux.h | 2 +- rugnux/rugnux_cli.cpp | 11 ++++++---- 9 files changed, 56 insertions(+), 15 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 18ca7f8c..da95c5f1 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -3,6 +3,7 @@ ### 1.0.0-rc.161 This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use. +* rugnux: The error-model **a** and **b** are reported in XDS's convention, with the whole-range ISa printed next to the strong-reflection asymptote. * rugnux: The background ring can be **elongated radially per reflection** for broadband data (`--integration-stencil `, default 0 = the fixed circular ring), by `k` times the beam's radial streak; the `r1` signal box stays circular. * rugnux: New **beam-stop shadow detection**, **on by default** (`--detect-beam-stop[=N|off]`), finds the beam stop and its holder in a projection of N images (default 60) and adds them to the pixel mask as bit 9, which is cleared at the start of every run. * Viewer: the detected beam-stop shadow is drawn in coral, with a "Show beam stop" switch in the side panel. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index b80ae089..3aa16a38 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -717,7 +717,9 @@ Per-shell and overall merging statistics are computed on corrected intensities, - completeness against the enumerated reflections for the cell and symmetry, - the anomalous signal-to-noise $\mathrm{SigAno}$ (below). -The error model is refined as $\sigma_\mathrm{corr}^2 = a\,\sigma^2 + (b\,\langle I\rangle)^2$, with $a$ set by the scatter of weak (counting-limited) reflections and $b$ the intensity-proportional systematic scatter of the strong ones. On the **rotation** path, **ISa** is the asymptotic ($I\to\infty$) signal-to-noise — by definition the reproducibility limit of the strongest reflections (Diederichs, *Acta Cryst.* **D66** (2010) 733) — and is read directly from the strong symmetry equivalents as the counting-subtracted fractional scatter of well-measured reflection groups (a robust median over strong groups; the $I/\sigma$ threshold is relaxed on weak or radiation-damaged data that has few strong reflections), rather than as $1/b$ of the whole-range fit, whose $b$ is raised slightly by an intermediate-intensity excess and so understates the limit. The reported ISa and the merged-intensity systematic floor $\sigma \ge b_\mathrm{ISa}\,|I|$ both use this asymptotic value, so a high-multiplicity merged $I/\sigma$ approaches ISa; the per-observation $\sigma_\mathrm{corr}$ (the merge weights) uses the whole-range $a,b$ and is unchanged. The **stills** path has no asymptotic estimate: it reports $\mathrm{ISa}=1/b$ and floors the merged $\sigma$ with the same whole-range $b$. +The error model is refined as $\sigma_\mathrm{corr}^2 = a\,\sigma^2 + (b\,\langle I\rangle)^2$, with $a$ set by the scatter of weak (counting-limited) reflections and $b$ the intensity-proportional systematic scatter of the strong ones. On the **rotation** path, **ISa** is the asymptotic ($I\to\infty$) signal-to-noise — by definition the reproducibility limit of the strongest reflections (Diederichs, *Acta Cryst.* **D66** (2010) 733) — and is read directly from the strong symmetry equivalents as the counting-subtracted fractional scatter of well-measured reflection groups (a robust median over strong groups; the $I/\sigma$ threshold is relaxed on weak or radiation-damaged data that has few strong reflections), rather than as $1/b$ of the whole-range fit, whose $b$ is raised slightly by an intermediate-intensity excess and so understates the limit. The asymptotic value is **report-only** — nothing downstream reads it, and the merged $\sigma$ is not floored at $b|I|$ (that floor was removed). The per-observation $\sigma_\mathrm{corr}$ (the merge weights) uses the whole-range $a,b$. The **stills** path has no asymptotic estimate and reports $\mathrm{ISa}=1/b$ directly. + +$a$ and $b$ are **reported in XDS's convention**, which is $\sigma^2 = a(\sigma_0^2 + b I^2)$ with $\mathrm{ISa}=1/\sqrt{ab}$, so the printed pair can be read straight against a `CORRECT.LP`. The internal fit keeps the form above; only the report converts, as $b_\mathrm{XDS} = b^2/a$. Note that $a$ is the same in both conventions and that the two ISa expressions are the same number, $1/\sqrt{a\cdot b^2/a} = 1/b$ — so the rotation log prints **two** ISa, the whole-range $1/b$ (XDS's meaning) and the strong-reflection asymptote beside it, which can only ever be the more optimistic of the two. A third, unrelated $b$ appears in the space-group search (§11); it is fitted with the $\sigma^2$ coefficient held at 1 and its gate constants are calibrated in that convention. **Anomalous signal-to-noise (SigAno).** The strength of the anomalous signal is reported per shell and overall as $\mathrm{SigAno}=\langle|\Delta I|\rangle / \langle\sigma(\Delta I)\rangle$, where $\Delta I = I(+)-I(-)$ over acentric reflections measured in both Bijvoet hands and $\sigma(\Delta I)=\sqrt{\sigma_+^2+\sigma_-^2}$. It is computed from the **full-multiplicity** inverse-variance $I(+)/I(-)$ split (the same one written to the output), i.e. from all observations rather than a half-set. For pure noise $\mathrm{SigAno}$ approaches the half-normal value $\sqrt{2/\pi}\approx0.8$, and it rises above $1$ once a real anomalous difference is present. A half-set anomalous correlation ("$\mathrm{CC}_\mathrm{anom}$") is **not** reported. Its two half estimates $\Delta I_0,\Delta I_1$ are complementary partitions of one observation pool ($\Delta I_0+\Delta I_1=2\,\Delta I_\mathrm{full}$), and subtracting the two Bijvoet hands cancels the large common intensity that keeps $\mathrm{CC}_{1/2}$ non-negative, leaving the small anomalous signal against the per-half split noise; below an anomalous signal-to-noise of $1$ per half that correlation tends towards $-1$ rather than $0$. $\mathrm{SigAno}$ has no such floor. It is emitted only when an anomalous split was made, using the standard PDBx items `_reflns.pdbx_absDiff_over_sigma_anomalous` (overall) and `_reflns_shell.pdbx_absDiff_over_sigma_anomalous` (per shell), and appears as the `SigAno` column of the printed merge-statistics table. @@ -848,7 +850,7 @@ A **dataset-wide** Wilson $B$ is also estimated over the merged reflections — 2. **Error-model $b$** (the intensity-proportional systematic). A genuine symmetry step gains multiplicity without inflating $b$; merging a twin law's extra operator inflates it. A $\chi^2$-passing promotion is vetoed when $b$ rises past a bound relative to the confirmed subgroup. 3. **Operator disagreement**, a sigma-free statistic $H=\mathrm{median}\,|I_1-I_2|/(I_1+I_2)$, formed as the ratio of the operators a promotion *adds* to the parent's own, measured on the same reflections. Normalising against the parent divides out the systematic floor that symmetry mates carry on real data, which varies by crystal and by operator; a median is used because a twin perturbs every pair whereas a badly-measured minority perturbs only the tail. Where a candidate has several parents of the same order, it is judged against the worst of them, since a rival subgroup can itself contain the twin laws. - The operator correlations are taken on reflections above an $I/\sigma$ cut, and that cut is **capped at the merge's own $I/\sigma$ quantile** rather than applied as a fixed number. A merged $\sigma$ carries the systematic floor $b|I|$ (§10.4), so no reflection in a merge can read above $\mathrm{ISa}=1/b$; a fixed cut of 3 therefore selects nothing at all on a search merge whose ISa is below 3, leaving every operator correlation undefined and collapsing the point group to 1. The cap keeps at least the strongest quarter and is inert — the cut stays exactly 3.0 — on a healthy merge. + The operator correlations are taken on reflections above an $I/\sigma$ cut, and that cut is **capped at the merge's own $I/\sigma$ quantile** rather than applied as a fixed number. On a search merge whose ISa is below 3, a fixed cut of 3 selects nothing at all, leaving every operator correlation undefined and collapsing the point group to 1. The cap keeps at least the strongest quarter and is inert — the cut stays exactly 3.0 — on a healthy merge. Several space groups may share an absence pattern exactly. Where they do, the search scores them identically and **all of them are named** in the result rather than one being reported as the answer: some are enantiomorph pairs, which merged intensities cannot distinguish in principle, and others differ only by a screw condition that the centering condition already implies, so the screw has no observable signature at all. The representative reported first is the lowest space-group number, which is a convention and not a measurement. diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index dc250a15..8bd626f5 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -76,6 +76,27 @@ struct MergeAccum { size_t n_half[2] = {0, 0}; }; +// XDS's error-model convention (Diederichs, Acta Cryst. D66 (2010) 733) is sigma^2 = a*(sigma0^2 + +// b*I^2), reported with ISa = 1/sqrt(a*b). Jungfraujoch fits sigma^2 = a*sigma0^2 + (b*)^2 - the +// same `a`, but a `b` that is a FRACTION of the intensity - so b_xds = b^2/a, and the two ISa +// expressions are the same number: 1/sqrt(a * b^2/a) = 1/b. Only `b` needs converting, and only +// where it is reported: doing it here rather than in the fit leaves every merge weight untouched. +// +// NOTE this is NOT the `b` of SearchSpaceGroup's merge_systematic_b, which is a third, unrelated +// quantity (a fraction of I, fitted with no `a` at all) whose gate constants are calibrated in that +// convention. Do not "make them consistent". +struct XdsErrorModel { + double a = 1.0; + double b = 0.0; + double isa = 0.0; +}; + +inline XdsErrorModel ToXdsErrorModel(double a, double b) { + if (!(a > 0.0) || !(b > 0.0)) + return {a, 0.0, 0.0}; + return {a, b * b / a, 1.0 / b}; +} + class MergeOnTheFly { mutable std::mutex merged_mutex; diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 83696826..364a65b9 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -2066,9 +2066,16 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // as unmeasured rather than emit a spurious extreme. constexpr double MIN_ASYMPTOTIC_B = 0.01; // ISa cap 100 result.isa = error_model_b_asymptotic >= MIN_ASYMPTOTIC_B ? 1.0 / error_model_b_asymptotic : 0.0; - if (error_model_active) - logger.Info("Error model: a={:.3f} b={:.3f} ISa={:.1f} chi2={:.2f}", error_model_a, error_model_b, - result.isa, error_model_chi2); + if (error_model_active) { + // Reported in XDS's convention so the numbers can be read against a CORRECT.LP directly. + // Two ISa are printed and they are different quantities: the whole-range 1/sqrt(a*b), which + // is what XDS's ISa means, and the strong-reflection asymptote above, which XDS has no + // equivalent of and which can only ever be the more optimistic of the two. + const auto em = ToXdsErrorModel(error_model_a, error_model_b); + logger.Info("Error model (XDS convention): a={:.3f} b={:.3e} ISa={:.1f} " + "(strong-reflection asymptote {:.1f}) chi2={:.2f}", + em.a, em.b, em.isa, result.isa, error_model_chi2); + } AssignRfreeFlags(result.merged, x.GetSpaceGroupNumber().value_or(1), rfree_fraction); // French-Wilson (F, and F(+)/F(-) from the anomalous split) is deferred until after the anomalous diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 6464d664..a124dde6 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -42,7 +42,7 @@ public: struct Result { std::vector merged; MergeStatistics statistics; - double isa = 0.0; // 1/b of the fitted error model (0 if the model stayed at identity) + double isa = 0.0; // 1/b_asymptotic: the strong-reflection tier, NOT the whole-range 1/b (0 if the model stayed at identity) }; // experiment: read live (its space group is changed by the caller between Run() calls). diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index 0c0d1fd0..8f264ca3 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -463,6 +463,10 @@ SearchSpaceGroupResult SearchSpaceGroup( // the systematic part of the scatter, which the fixed-sigma chi^2 ratio cannot: a genuine but // imperfectly-scaled high-symmetry merge and a twin can share a chi^2 ratio (~2) yet differ // sharply here (a genuine cubic step b x1.04 vs a merohedral twin b x1.6). + // This `b` is a fraction of I fitted with the sigma^2 coefficient held at 1 - it is NOT the merge + // error model's b, and NOT XDS's b either. The gate constants below are calibrated in this + // convention; converting them to any other silently squares the ratios and makes the absolute + // floor a-dependent, on a quantity that has no a. Leave it alone. auto merge_systematic_b = [&](const Orbits& orb) -> double { const auto& grp = orb.grp; const auto& rep = orb.rep; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index dc9cb468..22dd6dcb 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1436,10 +1436,13 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b && experiment_.GetScalingSettings().GetMinCCForImage() > 0.0); merge_engine.RefineErrorModel(merge_input); if (merge_engine.ErrorModelActive()) - logger.Info("Error model: a={:.3f} b={:.3f} ISa={:.1f} chi2={:.2f}", merge_engine.ErrorModelA(), - merge_engine.ErrorModelB(), - merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0, - merge_engine.ErrorModelChi2()); + { + // XDS convention, so the pair reads directly against a CORRECT.LP; the stills ISa is + // already the whole-range 1/b, there being no asymptote tier on this path. + const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB()); + logger.Info("Error model (XDS convention): a={:.3f} b={:.3e} ISa={:.1f} chi2={:.2f}", + em.a, em.b, em.isa, merge_engine.ErrorModelChi2()); + } for (size_t i = 0; i < merge_input.size(); ++i) merge_engine.AddImage(merge_input[i], static_cast(i)); diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index dcf5ff64..00cf9dd1 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -123,7 +123,7 @@ struct ProcessResult { std::string merge_statistics_text; // populated when scaling/merging ran // Structured merge statistics (per-shell + overall), populated when scaling/merging ran. ISa is - // the error-model asymptotic I/sigma (1/b); has_reference is true when a reference MTZ drove CCref. + // the error-model asymptotic I/sigma (the strong-reflection tier; the whole-range XDS-equivalent is 1/b); has_reference is true when a reference MTZ drove CCref. bool has_merge_statistics = false; MergeStatistics merge_statistics; double error_model_isa = 0.0; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 96d2b060..15e32116 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1396,10 +1396,13 @@ static int RunRugnux(int argc, char **argv) { // model and produce much worse stills intensities (no (b*I)^2 systematic term at all). merge_engine.RefineErrorModel(reflections); if (merge_engine.ErrorModelActive()) - logger.Info("Error model: a={:.3f} b={:.3f} ISa={:.1f} chi2={:.2f}", merge_engine.ErrorModelA(), - merge_engine.ErrorModelB(), - merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0, - merge_engine.ErrorModelChi2()); + { + // XDS convention, so the pair reads directly against a CORRECT.LP; the stills ISa is + // already the whole-range 1/b, there being no asymptote tier on this path. + const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB()); + logger.Info("Error model (XDS convention): a={:.3f} b={:.3e} ISa={:.1f} chi2={:.2f}", + em.a, em.b, em.isa, merge_engine.ErrorModelChi2()); + } for (size_t i = 0; i < reflections.size(); ++i) merge_engine.AddImage(reflections[i], static_cast(i)); merged_reflections = merge_engine.ExportReflections();