From adf87e86758982c8962635e8426e3ca1408ddeae Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 10 Aug 2026 18:55:46 +0200 Subject: [PATCH] Merging: export the XDS-comparable ISa under jfjoch_diffrn_ISa The mmCIF's _reflns.jfjoch_diffrn_ISa carried the strong-reflection asymptote, a tier XDS has no equivalent of, while the name invites comparison with XDS's ISa - which is the whole-range 1/sqrt(a*b). rugnux_vs_xds.py reads that item for the battery's ISa column, so the comparison that column exists to make was between two different quantities, flattering rugnux by the difference between the tiers. Write the whole-range value there, move the asymptote to _reflns.jfjoch_diffrn_ISa_asymptotic, and add _reflns.jfjoch_error_model_a and _b in XDS's convention so the number can be re-derived from the file rather than taken on trust. On a broadband rotation dataset the battery column now reads 13.25 against XDS's 21.18 where it read 15.6 before, and the two error models can be compared term by term for the first time: a 1.538 vs 1.249 and b 3.71e-03 vs 1.78e-03, so the gap is in BOTH the counting and the systematic term (1.23x and 2.08x, and sqrt(1.23*2.08) = 1.60 = 21.18/13.25). This is a deliberate redefinition of an exported item, not an addition: a file written by an earlier version carries the asymptote under the old name and there is no version marker to tell them apart. Noted in the changelog and in docs/CPU_DATA_ANALYSIS.md. Nothing reads the item back into the pipeline - it is written and never parsed by rugnux itself - so no stored file is reinterpreted in a way that changes a result. Co-Authored-By: Claude Opus 5 (1M context) --- docs/CHANGELOG.md | 2 +- docs/CPU_DATA_ANALYSIS.md | 2 +- image_analysis/WriteReflections.cpp | 20 ++++++++++++--- image_analysis/WriteReflections.h | 15 +++++++++-- .../scale_merge/RotationScaleMerge.cpp | 10 ++++++-- .../scale_merge/RotationScaleMerge.h | 9 ++++++- rugnux/Rugnux.cpp | 19 +++++++++++--- rugnux/Rugnux.h | 5 +++- rugnux/rugnux_cli.cpp | 25 ++++++++++++++++--- 9 files changed, 88 insertions(+), 19 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index da95c5f1..be240e67 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -3,7 +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 error-model **a** and **b** are reported in XDS's convention, and `_reflns.jfjoch_diffrn_ISa` now carries the whole-range `1/sqrt(a*b)` that XDS's ISa denotes; the strong-reflection asymptote moves to `_reflns.jfjoch_diffrn_ISa_asymptotic`. **A file written by an earlier version carries the asymptote under the old name.** * 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 3aa16a38..2b792ac9 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -719,7 +719,7 @@ Per-shell and overall merging statistics are computed on corrected intensities, 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. +$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. The mmCIF follows the same split: `_reflns.jfjoch_diffrn_ISa` is the whole-range value, directly comparable with a `CORRECT.LP`, and the asymptote is written separately as `_reflns.jfjoch_diffrn_ISa_asymptotic`, with `_reflns.jfjoch_error_model_a` and `_b` alongside so the number can be re-derived. Note that a file written before this change carries the *asymptote* under the plain `ISa` name. 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. diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index 334445c5..73ddedce 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -148,7 +148,7 @@ void WriteMmcifReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, const MergeStatistics &statistics, - const std::string &isa, + const ErrorModelReport &error_model, const TwinningAnalysisResult &twinning, const std::string &filename) { @@ -237,7 +237,19 @@ void WriteMmcifReflections(const std::vector &reflections, if (has_anom) out << "_reflns.pdbx_absDiff_over_sigma_anomalous " << Fmt(ov.abs_diff_over_sigma_anomalous, 3) << " # SigAno = <|dano|>/\n"; - out << "_reflns.jfjoch_diffrn_ISa " << CifStr(isa) << " # asymptotic I/sigma (Diederichs)\n"; + // ISa in XDS's sense: the whole-range 1/sqrt(a*b) of the error model below, so this item can + // be read straight against a CORRECT.LP. The strong-reflection asymptote - a tier XDS does not + // have, and always the more optimistic of the two - is written separately rather than here. + out << "_reflns.jfjoch_diffrn_ISa " << CifStr(error_model.isa) + << " # 1/sqrt(a*b), the XDS convention\n"; + if (!error_model.isa_asymptotic.empty()) + out << "_reflns.jfjoch_diffrn_ISa_asymptotic " << CifStr(error_model.isa_asymptotic) + << " # strong-reflection asymptote (Diederichs); rotation path only\n"; + if (!error_model.a.empty()) + out << "_reflns.jfjoch_error_model_a " << CifStr(error_model.a) + << " # sigma^2 = a*(sigma0^2 + b*I^2), XDS convention\n"; + if (!error_model.b.empty()) + out << "_reflns.jfjoch_error_model_b " << CifStr(error_model.b) << "\n"; // Dataset-wide isotropic Wilson B-factor estimate (standard PDBx item), analogous to XDS's // "WILSON LINE ... B=". Emitted only when the log-linear fit succeeded. if (std::isfinite(statistics.wilson_b) && statistics.wilson_b > 0.0) @@ -468,13 +480,13 @@ void WriteReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, const MergeStatistics &statistics, - const std::string &isa, + const ErrorModelReport &error_model, const TwinningAnalysisResult &twinning, const std::string &filename) { // Write an MTZ, an mmCIF and a SHELX HKLF-4 .hkl - each has its uses downstream (MTZ for the CCP4 / // phenix reflection tools, mmCIF for deposition and as the self-describing native format, HKLF-4 as // the SHELXC / ANODE substructure-solution input). WriteMtzReflections(reflections, unitCell, experiment, filename + ".mtz"); - WriteMmcifReflections(reflections, unitCell, experiment, statistics, isa, twinning, filename + ".cif"); + WriteMmcifReflections(reflections, unitCell, experiment, statistics, error_model, twinning, filename + ".cif"); WriteShelxHklReflections(reflections, experiment, filename + ".hkl"); } diff --git a/image_analysis/WriteReflections.h b/image_analysis/WriteReflections.h index c6537e81..3cf42c64 100644 --- a/image_analysis/WriteReflections.h +++ b/image_analysis/WriteReflections.h @@ -13,11 +13,22 @@ struct MergeStatistics; struct TwinningAnalysisResult; +// The error model as it is reported, already formatted. `isa` is the whole-range 1/sqrt(a*b), the +// same quantity XDS's ISa denotes, so a file written here is directly comparable with a CORRECT.LP; +// `isa_asymptotic` is the strong-reflection tier, which only the rotation path has. `a` and `b` are +// in XDS's convention, sigma^2 = a*(sigma0^2 + b*I^2). Empty strings are written as unknown. +struct ErrorModelReport { + std::string isa; + std::string isa_asymptotic; + std::string a; + std::string b; +}; + void WriteMmcifReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, const MergeStatistics &statistics, - const std::string &isa, + const ErrorModelReport &error_model, const TwinningAnalysisResult &twinning, const std::string &filename); @@ -35,6 +46,6 @@ void WriteReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, const MergeStatistics &statistics, - const std::string &isa, + const ErrorModelReport &error_model, const TwinningAnalysisResult &twinning, const std::string &filename); \ No newline at end of file diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 364a65b9..ad7516a0 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -2065,7 +2065,13 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // I/sigma. Real macromolecular data does not exceed ISa ~50; past a generous cap report the asymptote // 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; + result.isa_asymptotic = error_model_b_asymptotic >= MIN_ASYMPTOTIC_B ? 1.0 / error_model_b_asymptotic : 0.0; + { + const auto em = ToXdsErrorModel(error_model_a, error_model_b); + result.isa = em.isa; + result.error_model_a = em.a; + result.error_model_b = em.b; + } 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 @@ -2074,7 +2080,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool 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); + em.a, em.b, em.isa, result.isa_asymptotic, error_model_chi2); } AssignRfreeFlags(result.merged, x.GetSpaceGroupNumber().value_or(1), rfree_fraction); diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index a124dde6..86f64461 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -42,7 +42,14 @@ public: struct Result { std::vector merged; MergeStatistics statistics; - double isa = 0.0; // 1/b_asymptotic: the strong-reflection tier, NOT the whole-range 1/b (0 if the model stayed at identity) + // Two tiers, and they are different quantities. `isa` is the whole-range 1/sqrt(a*b) - which + // in this parameterisation is 1/b - and is what XDS's ISa means, so it is the one exported. + // `isa_asymptotic` is the strong-reflection tier, which XDS has no equivalent of and which can + // only ever be the more optimistic of the two. Both 0 if the model stayed at identity. + double isa = 0.0; + double isa_asymptotic = 0.0; + double error_model_a = 0.0; // XDS convention: sigma^2 = a*(sigma0^2 + b*I^2) + double error_model_b = 0.0; }; // experiment: read live (its space group is changed by the caller between Run() calls). diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 22dd6dcb..c24ec72c 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1385,6 +1385,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b phase("Scale/combine/merge (" + label + ")"); auto r = rsm->Run(for_search); result.error_model_isa = r.isa; + result.error_model_isa_asymptotic = r.isa_asymptotic; + result.error_model_a = r.error_model_a; + result.error_model_b = r.error_model_b; return ScaleMergeResult{std::move(r.merged), std::move(r.statistics)}; } // Stills (rotation goes through RotationScaleMerge above): self-scale each image against the @@ -1461,7 +1464,12 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b phase("Computing statistics"); out.statistics = merge_engine.MergeStats(out.merged, merge_input, config_.reference_data, effective_d_min); - result.error_model_isa = merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0; + { // stills: no asymptote tier, and its ISa already IS the whole-range 1/b + const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB()); + result.error_model_isa = em.isa; + result.error_model_a = em.a; + result.error_model_b = em.b; + } logger.Info("Merge complete ({} unique reflections, {})", out.merged.size(), label); return out; }; @@ -2020,9 +2028,14 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b if (result.consensus_cell && write_files && config_.write_merged) { phase("Writing reflections"); + const ErrorModelReport em_report{ + result.error_model_isa > 0 ? fmt::format("{:.2f}", result.error_model_isa) : "?", + result.error_model_isa_asymptotic > 0 ? fmt::format("{:.2f}", result.error_model_isa_asymptotic) + : std::string(), + result.error_model_a > 0 ? fmt::format("{:.3f}", result.error_model_a) : std::string(), + result.error_model_b > 0 ? fmt::format("{:.4e}", result.error_model_b) : std::string()}; WriteReflections(sm.merged, *result.consensus_cell, experiment_, sm.statistics, - result.error_model_isa > 0 ? fmt::format("{:.2f}", result.error_model_isa) : "?", - result.twinning, config_.output_prefix); + em_report, result.twinning, config_.output_prefix); // Per-image scaling table (G, B-factor, mosaicity, wedge, CC) for inspection / XDS // comparison. The offline self-scaling result is otherwise not exposed (process.h5's // per-image arrays are only filled on the online per-image path). Sourced from the diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 00cf9dd1..cb5291ec 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -126,7 +126,10 @@ struct ProcessResult { // 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; + double error_model_isa = 0.0; // whole-range 1/sqrt(a*b), XDS's meaning + double error_model_isa_asymptotic = 0.0; // strong-reflection tier; rotation path only, 0 on stills + double error_model_a = 0.0; // XDS convention: sigma^2 = a*(sigma0^2 + b*I^2) + double error_model_b = 0.0; bool has_reference = false; // Per-reflection (I, sigma) of the final merged reflections, for the ISa diagnostic: I/sigma diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 15e32116..8d801f88 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1352,6 +1352,9 @@ static int RunRugnux(int argc, char **argv) { std::vector merged_reflections; MergeStatistics merged_statistics; double error_model_isa = 0.0; + double error_model_isa_asymptotic = 0.0; + double error_model_a = 0.0; + double error_model_b = 0.0; // Rotation (rot3d): the dedicated RotationScaleMerge does the whole self-scale -> 3D combine -> // merge, including the default-on decay + absorption correction surfaces. It does not support @@ -1372,6 +1375,9 @@ static int RunRugnux(int argc, char **argv) { merged_reflections = std::move(r.merged); merged_statistics = std::move(r.statistics); error_model_isa = r.isa; + error_model_isa_asymptotic = r.isa_asymptotic; + error_model_a = r.error_model_a; + error_model_b = r.error_model_b; } else { // Scaling self-references: the reference MTZ (if any) fixes the cell/space group, reports // CCref and provides the R-free test set, but is NOT a scale anchor - scaling each image @@ -1418,7 +1424,12 @@ static int RunRugnux(int argc, char **argv) { merged_statistics = merge_engine.MergeStats(merged_reflections, reflections, reference_data, effective_d_min); - error_model_isa = merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0; + { // stills: no asymptote tier, and its ISa already IS the whole-range 1/b + const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB()); + error_model_isa = em.isa; + error_model_a = em.a; + error_model_b = em.b; + } } logger.Info("Scale + merge completed in {:.2f} s ({} unique reflections)", @@ -1451,10 +1462,16 @@ static int RunRugnux(int argc, char **argv) { const auto twinning = AnalyzeTwinning(merged_reflections, twin_sg); std::cout << std::endl << TwinningAnalysisToText(twinning) << std::endl; - if (!output_prefix.empty()) + if (!output_prefix.empty()) { + const ErrorModelReport em_report{ + error_model_isa > 0 ? fmt::format("{:.2f}", error_model_isa) : "?", + error_model_isa_asymptotic > 0 ? fmt::format("{:.2f}", error_model_isa_asymptotic) + : std::string(), + error_model_a > 0 ? fmt::format("{:.3f}", error_model_a) : std::string(), + error_model_b > 0 ? fmt::format("{:.4e}", error_model_b) : std::string()}; WriteReflections(merged_reflections, *experiment.GetUnitCell(), experiment, merged_statistics, - error_model_isa > 0 ? fmt::format("{:.2f}", error_model_isa) : "?", - twinning, output_prefix); + em_report, twinning, output_prefix); + } if (!output_prefix.empty() && !model_pdb.empty()) { const auto data_sg = experiment.GetSpaceGroupNumber();