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();