From 3caa11caa4a7b8f5818f037619b06eeaa7edb2c2 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 1 Nov 2026 10:34:55 +0100 Subject: [PATCH] Stills: the error model is reported as its terms a, b, c, not as an ISa A stills merge no longer reports ISA. On stills 1/b is the I/sigma ceiling of a single still observation under the stills model - the same frames processed as a rotation sweep read several times higher - so it describes the processing model more than the instrument and the crystal, and read on CORRECT.LP's scale it calls good serial data useless. The report gives the error model as ERROR_MODEL_A, ERROR_MODEL_B (XDS convention, as before) and a new ERROR_MODEL_C, the partiality term, in the same convention: sigma^2 = a*(sigma0^2 + b*I^2 + c*((1-p)/p)^2*I^2). ERROR_MODEL_B= undetermined where b is zero within its error (ISA= said it before). The same on every other place that showed the stills ISa: the summary's Signal row, the mmCIF (_reflns.jfjoch_error_model_c instead of _reflns.jfjoch_diffrn_ISa), the log line, the live status the broker serves, and the viewer's merge-statistics window (no ISa card and no ISa diagnostic for a stills merge). MergeStatistics::stills marks a merge made by the stills path. The 1/b the space-group search reads (merge_isa) is unchanged. Rotation output is unchanged. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- docs/CPU_DATA_ANALYSIS_INTEGRATION.md | 4 +- docs/RUGNUX_INTEGRATION.md | 3 +- docs/RUGNUX_REPORT.md | 14 ++++- image_analysis/WriteReflections.cpp | 8 ++- image_analysis/WriteReflections.h | 5 +- image_analysis/scale_merge/Merge.cpp | 1 + image_analysis/scale_merge/Merge.h | 10 ++++ receiver/JFJochReceiverRugnux.cpp | 2 +- rugnux/ResultReport.cpp | 51 ++++++++++++++----- rugnux/Rugnux.h | 6 ++- rugnux/RugnuxScaleMerge.cpp | 13 +++-- rugnux/rugnux_cli.cpp | 15 ++++-- tests/ResultReportTest.cpp | 26 ++++++++++ tools/battery/still/rugnux_stills_ab.py | 4 +- viewer/windows/JFJochMergeStatsWindow.cpp | 5 +- viewer/windows/JFJochMergeStatsWindow.h | 2 +- viewer/windows/JFJochProcessingJobsWindow.cpp | 5 +- 17 files changed, 135 insertions(+), 39 deletions(-) diff --git a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md index 61fac7826..afb55d2be 100644 --- a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md +++ b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md @@ -388,9 +388,9 @@ Per-shell and overall merging statistics are computed on corrected intensities, itself count as missing instead of leaving the denominator along with the data, - the anomalous signal-to-noise $\mathrm{SigAno}$ and the half-set anomalous correlation $\mathrm{CC}_\mathrm{anom}$ (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 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. +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 adds a partiality term, $\sigma_\mathrm{corr}^2 = a\,\sigma^2 + (b\,\langle I\rangle)^2 + c^2\left(\frac{1-p}{p}\right)^2\langle I\rangle^2$, and reports **no ISa**: there $1/b$ is the $I/\sigma$ ceiling of one fully recorded still under the stills model - the same frames processed as a rotation sweep read several times higher - which describes the processing model more than the instrument and the crystal. It reports $a$, $b$ and $c$ instead (`ERROR_MODEL_A/B/C`, `_reflns.jfjoch_error_model_a/b/c`); $1/b$ is still used inside, by the space-group search. -$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 (§13.1); 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$ (and on stills $c_\mathrm{XDS} = c^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 (§13.1); 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 reported beside it: $\Delta I$ is formed once per half-set and the two are correlated over the acentric pairs where **both** hands split into two non-empty halves, per shell and overall as one correlation rather than a mean of shells. Unlike SigAno it is not a ratio against the error model, so an optimistic $\sigma$ cannot inflate it. It has no floor either: subtracting the two Bijvoet hands cancels the large common intensity that keeps $\mathrm{CC}_{1/2}$ non-negative, so on data with little anomalous signal and about two observations per mate it goes strongly negative. That is a property of the statistic and is reported as measured. It agrees with AIMLESS's `CCanom` and `phenix.merging_statistics`' `cc_anom`; XDS's `Anomal Corr` is a **different quantity and is not comparable with it** — on the same observations it reads two to three times higher in the low shells. $\mathrm{CC}_\mathrm{anom}$ is a rotation-path statistic (the stills merge forms no half-set anomalous difference), and it is reported in the `CCanom` column of the printed merge-statistics table and as `CC_ANOM=` in the results report — not in the mmCIF. SigAno 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 same table. Where either could not be measured at all — a Friedel-merged run that split no Bijvoet pair, a shell too thin to split one in both hands — the table prints `-` and the results report writes no key, which is not the same statement as a value measured to be zero. diff --git a/docs/RUGNUX_INTEGRATION.md b/docs/RUGNUX_INTEGRATION.md index b1c788a5d..c7cf964c8 100644 --- a/docs/RUGNUX_INTEGRATION.md +++ b/docs/RUGNUX_INTEGRATION.md @@ -24,9 +24,10 @@ ignores them, and one that does can find them without guessing. | item | meaning | |---|---| -| `_reflns.jfjoch_diffrn_ISa` | Asymptotic I/σ in **XDS's sense**: the whole-range `1/√(a·b)` of the error model, so it can be read directly against a `CORRECT.LP` | +| `_reflns.jfjoch_diffrn_ISa` | Asymptotic I/σ in **XDS's sense**: the whole-range `1/√(a·b)` of the error model, so it can be read directly against a `CORRECT.LP`. Rotation only: a stills merge writes the error-model terms instead | | `_reflns.jfjoch_diffrn_ISa_asymptotic` | The **strong-reflection** tier — the counting-subtracted scatter of well-measured groups. XDS has no equivalent, and it can only ever be the more optimistic of the two. Rotation path only | | `_reflns.jfjoch_error_model_a`, `_b` | The error model in XDS's convention, `σ² = a(σ₀² + b·I²)`, so the ISa above is re-derivable from the file rather than taken on trust | +| `_reflns.jfjoch_error_model_c` | Stills only: the partiality term, `σ² = a(σ₀² + b·I² + c·((1−p)/p)²·I²)` with p the observation's partiality | | `_reflns.jfjoch_second_moment_I` | Twinning second moment ⟨I²⟩/⟨I⟩² — 2.00 untwinned, 1.50 for a perfect twin | | `_reflns.jfjoch_L_test_mean_abs_L`, `_L_test_mean_L_squared` | Padilla–Yeates L-test. ⟨\|L\|⟩ is 0.500 untwinned / 0.375 for a perfect twin; ⟨L²⟩ is 0.333 / 0.200. Written only when the test found pairs | | `_reflns.jfjoch_radiation_damage_relative_B` | Relative *B* from the first to the last rotation batch (Ų); positive is the usual direction, high-resolution intensity fading with dose | diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 9994522b6..b389b81fb 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -231,10 +231,22 @@ The keys most scripts read, by section; the rest of this page explains the ones | `UNIQUE_REFLECTIONS`, `TOTAL_OBSERVATIONS`, `COMPLETENESS`, `MULTIPLICITY` | Over the written range; completeness in percent | | `I_OVER_SIGMA`, `R_MEAS`, `CC_HALF`, `WILSON_B` | The overall merging statistics; the shell table beside them is in the comment lines | | `SIGANO`, `CC_ANOM`, `FRIEDELS_LAW` | The anomalous signal, and whether Friedel pairs were merged (`FALSE` with `-A`) | -| `ERROR_MODEL_A`, `ERROR_MODEL_B`, `ISA` | The error model in XDS's convention and ISa = 1/sqrt(a·b); `ISA= undetermined` where b is zero within its error, so 1/b would be noise | +| `ERROR_MODEL_A`, `ERROR_MODEL_B`, `ISA` | The error model in XDS's convention and ISa = 1/sqrt(a·b); `ISA= undetermined` where b is zero within its error, so 1/b would be noise. Rotation only for `ISA` | +| `ERROR_MODEL_C` | Stills only: the partiality term of the stills error model, σ² = a(σ₀² + b·I² + c·((1−p)/p)²·I²) in the same convention. A stills report writes no `ISA` ([below](#the-stills-error-model)), and `ERROR_MODEL_B= undetermined` where b is zero within its error | | `R_WORK`, `R_FREE` | Section 5: the model as placed against these data ([below](#model-validation)); the refined model's are `REFINE_R_WORK`/`REFINE_R_FREE` | | `R_FREE_FLAGS` | Section 5: the free set both are read on - `DRAWN` by this run, inherited from the `-z` `REFERENCE`, or from the `--free-flags` `FILE` | +### The stills error model + +A stills merge reports its error model as its terms, `ERROR_MODEL_A`, `ERROR_MODEL_B` and +`ERROR_MODEL_C`, with σ² = a(σ₀² + b·I² + c·((1−p)/p)²·I²), p being the partiality of the observation: +b is the systematic error of a fully recorded still and c the error a partial one adds. It writes no +`ISA`. On stills 1/sqrt(a·b) is the I/σ ceiling of a single still observation under the stills model - +the same frames processed as a rotation sweep read several times higher - so it describes the +processing model more than the instrument and the crystal, and read on `CORRECT.LP`'s scale it would +call good serial data useless. The serial figures of merit are `R_SPLIT=`, `CC_HALF=` and `CC_STAR=`. +The mmCIF follows: `_reflns.jfjoch_error_model_c` instead of `_reflns.jfjoch_diffrn_ISa`. + ## The summary, and what the run decided The file opens with a **`SUMMARY`** section, above everything it summarises. It exists because the diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index 8bb1b8789..b30d77a2c 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -428,8 +428,9 @@ void WriteMmcifReflections(const std::vector &reflections, // 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 (!statistics.stills) + 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"; @@ -438,6 +439,9 @@ void WriteMmcifReflections(const std::vector &reflections, << " # 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"; + if (!error_model.c.empty()) + out << "_reflns.jfjoch_error_model_c " << CifStr(error_model.c) + << " # stills partiality term: + c*((1-p)/p)^2*I^2\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) diff --git a/image_analysis/WriteReflections.h b/image_analysis/WriteReflections.h index 14b5be3c6..badd3d60f 100644 --- a/image_analysis/WriteReflections.h +++ b/image_analysis/WriteReflections.h @@ -21,12 +21,15 @@ struct CellEsd; // geom_refinement/CellEsd.h // 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. +// in XDS's convention, sigma^2 = a*(sigma0^2 + b*I^2), and `c` is the stills partiality term in the same +// convention, sigma^2 = a*(sigma0^2 + b*I^2 + c*((1-p)/p)^2*I^2). Empty strings are written as unknown. +// A stills merge writes no ISa (MergeStatistics::stills). struct ErrorModelReport { std::string isa; std::string isa_asymptotic; std::string a; std::string b; + std::string c; }; // What the reflection files record about where they came from, beyond the experiment itself and the diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 12ac2add6..4c594d4f3 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -949,6 +949,7 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m }; MergeStatistics out; + out.stills = true; out.shells.resize(n_shells); out.n_observations_rejected = reject_count; diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 8f9ee1fa0..56ea09bca 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -182,6 +182,10 @@ struct MergeStatistics { std::vector shells; MergeStatisticsShell overall; AnomalousMergeStatistics per_hand; + // Made by the stills merge (MergeOnTheFly). Its error model is reported as its terms, not as an + // ISa: on stills 1/b is the ceiling of one still observation under the stills model, which is not + // the instrument and crystal error CORRECT.LP's ISa measures. + bool stills = false; // Dataset-wide isotropic Wilson B-factor estimate (A^2) from the log-linear fit of the shell-mean // merged intensity against 1/d^2 (CalcGlobalWilsonB) - the analogue of XDS's "WILSON LINE ... B=". @@ -280,6 +284,12 @@ inline XdsErrorModel ToXdsErrorModel(double a, double b) { return {a, b * b / a, 1.0 / b}; } +// The stills partiality term c (MergeOnTheFly::ErrorModelC) in the same convention: +// sigma^2 = a*(sigma0^2 + b*I^2 + c*((1-p)/p)^2*I^2). +inline double ToXdsPartialityTerm(double a, double c) { + return a > 0.0 ? c * c / a : 0.0; +} + // Why the stills merge took or left one observation: a bit for every test it fails, 0 = merged // (MergeOnTheFly::Verdicts; a diagnostic, --export-observations). enum ObservationFlag : uint32_t { diff --git a/receiver/JFJochReceiverRugnux.cpp b/receiver/JFJochReceiverRugnux.cpp index 08cab131f..6afcbd87a 100644 --- a/receiver/JFJochReceiverRugnux.cpp +++ b/receiver/JFJochReceiverRugnux.cpp @@ -426,7 +426,7 @@ void JFJochReceiverRugnux::RugnuxThread() { std::optional d_min, isa; if (result.has_merge_statistics) SetIfFinite(d_min, result.merge_statistics.overall.d_min); - if (result.error_model_isa_resolved && result.error_model_isa > 0.0) + if (!result.merge_statistics.stills && result.error_model_isa_resolved && result.error_model_isa > 0.0) isa = result.error_model_isa; live_status.Headline(result.space_group ? result.space_group->xhm() : "", d_min, isa); // The placed model and a copy of its map coefficients (the end message takes the originals), and the diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 132ece884..7b4cf9ee7 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -47,6 +47,8 @@ namespace { // the XDS GEOMETRY KEYWORDS section; OVERLOADED_* / OVERLAPPED_* and the OVERLOADS / OVERLAPS flags. // 17: INDEXING_AMBIGUITY_RESOLVED and the AMBIGUITY_* keys in section 2 - serial stills, the // indexing ambiguity resolved from the data alone (or measured against the reference's answer). + // 18: stills report their error model as ERROR_MODEL_A / _B / _C (C the partiality term) and no + // ISA (nor ISa in the summary's Signal row). // 8: every line that is not `KEY= value` and not blank starts with `#`, warnings included // (`# WARNING:`); the header gains RUGNUX_DOWNLOAD and BUILD_CXX_FLAGS, and RUGNUX_GIT is // the full build-time hash. @@ -1543,7 +1545,14 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, // The error model in XDS's convention, so the numbers are directly comparable with a // CORRECT.LP. Kept in the default report because beamline monitoring scripts read them. Add(s, KeyReal("ERROR_MODEL_A", result.error_model_a, "{:.4f}")); - Add(s, KeyReal("ERROR_MODEL_B", result.error_model_b, "{:.4e}")); + // On stills the terms are the whole report of the error model, so a b that is zero within + // its error says so here; on rotation ISA= says it. + if (ms.stills && !result.error_model_isa_resolved) + Add(s, KeyText("ERROR_MODEL_B", "undetermined")); + else + Add(s, KeyReal("ERROR_MODEL_B", result.error_model_b, "{:.4e}")); + if (ms.stills) + Add(s, KeyReal("ERROR_MODEL_C", result.error_model_c, "{:.4e}")); // A b that is zero within its error gives a 1/b that is the noise in b, not an I/sigma. Add(s, KeyText("SCALING_ITERATIONS", fmt::format("{} partials, {} fulls{}", result.scaling_iterations_partials, @@ -1558,12 +1567,17 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, ? "the space-group search merge and the final merge" : !result.scaling_converged ? "the final merge" : "the space-group search merge")); - if (result.error_model_isa_resolved) - Add(s, KeyReal("ISA", result.error_model_isa, "{:.2f}")); - else - Add(s, KeyText("ISA", "undetermined")); - if (result.error_model_isa_resolved && result.error_model_isa_asymptotic > 0.0) - Add(s, KeyReal("ISA_ASYMPTOTIC", result.error_model_isa_asymptotic, "{:.2f}", true)); + // No ISA on stills: there 1/b is the ceiling of one still observation under the stills + // model, a property of the processing more than of the instrument, and read on CORRECT.LP's + // scale it calls good serial data useless. + if (!ms.stills) { + if (result.error_model_isa_resolved) + Add(s, KeyReal("ISA", result.error_model_isa, "{:.2f}")); + else + Add(s, KeyText("ISA", "undetermined")); + if (result.error_model_isa_resolved && result.error_model_isa_asymptotic > 0.0) + Add(s, KeyReal("ISA_ASYMPTOTIC", result.error_model_isa_asymptotic, "{:.2f}", true)); + } // Only when TRUE: FALSE is the ordinary case and says nothing. if (result.has_reference) { Add(s, KeyBool("REFERENCE_DATA_USED", true)); @@ -1626,9 +1640,16 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, " both hands the quantity does not exist: the key is then absent here and the column is\n" " a dash in the table below, which is not the same claim as a signal measured to be\n" " zero.\n\n" - " ERROR_MODEL_A / ERROR_MODEL_B are in XDS's convention, sigma^2 = a*(sigma0^2 + b*I^2),\n" - " so ISA = 1/sqrt(a*b) means what CORRECT.LP's ISa means. ISA_ASYMPTOTIC, where present,\n" - " is the strong-reflection tier only.", true)); + + std::string(ms.stills + ? " ERROR_MODEL_A / ERROR_MODEL_B / ERROR_MODEL_C are the stills error model in XDS's\n" + " convention, sigma^2 = a*(sigma0^2 + b*I^2 + c*((1-p)/p)^2*I^2), with p the partiality\n" + " of the observation: b is the systematic error of a fully recorded still, c the added\n" + " error of a partial one. No ISa is given: 1/sqrt(a*b) here is the I/sigma ceiling of\n" + " a single still observation under the stills model, not what CORRECT.LP's ISa\n" + " measures; R_SPLIT, CC_HALF and CC_STAR are the numbers to quote." + : " ERROR_MODEL_A / ERROR_MODEL_B are in XDS's convention, sigma^2 = a*(sigma0^2 + b*I^2),\n" + " so ISA = 1/sqrt(a*b) means what CORRECT.LP's ISa means. ISA_ASYMPTOTIC, where present,\n" + " is the strong-reflection tier only."), true)); { std::ostringstream shells; @@ -3340,10 +3361,12 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, sig += fmt::format(" R_meas {:.1f} %", 100.0 * o.r_meas); if (std::isfinite(o.cc_half)) sig += fmt::format(" CC1/2 {:.3f}", o.cc_half); - if (!result.error_model_isa_resolved) - sig += " ISa undetermined"; - else if (result.error_model_isa > 0.0) - sig += fmt::format(" ISa {:.1f}", result.error_model_isa); + if (!result.merge_statistics.stills) { // no ISa on stills (section 3) + if (!result.error_model_isa_resolved) + sig += " ISa undetermined"; + else if (result.error_model_isa > 0.0) + sig += fmt::format(" ISa {:.1f}", result.error_model_isa); + } row("Signal", sig); if (std::isfinite(o.cc_anom)) row("Anomalous", o.cc_anom > CC_ANOM_SIGNAL diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 8f0e88824..f2baea81c 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -564,11 +564,15 @@ 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; // whole-range 1/sqrt(a*b), XDS's meaning + double error_model_isa = 0.0; // whole-range 1/sqrt(a*b), XDS's meaning; on stills 1/b of + // the stills model, used inside (space-group search), not reported double error_model_isa_asymptotic = 0.0; // strong-reflection tier; rotation path only, 0 on stills bool error_model_isa_resolved = true; // false: b is zero within its error, ISa printed "undetermined" double error_model_a = 0.0; // XDS convention: sigma^2 = a*(sigma0^2 + b*I^2) double error_model_b = 0.0; + // Stills only: the partiality term, in the same convention, + // sigma^2 = a*(sigma0^2 + b*I^2 + c*((1-p)/p)^2*I^2). + double error_model_c = 0.0; bool has_reference = false; // How the reference MTZ was matched: the operator the merged indices were put through to reach // its frame (h,k,l where none), the intensity CC and matched fraction in that frame, and whether diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 24f557e03..f5f03a14b 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -1788,11 +1788,12 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { merge_engine.RefineErrorModel(merge_input); if (merge_engine.ErrorModelActive()) { - // 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. + // XDS convention, sigma^2 = a*(sigma0^2 + b*I^2 + c*((1-p)/p)^2*I^2). No ISa: on stills 1/b is + // the ceiling of one still observation, not what CORRECT.LP's ISa measures. const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB()); - logger.Info("Error model (XDS convention): a={:.3f} b={:.3e} ISa={:.1f} partiality c={:.3f} chi2={:.2f}", - em.a, em.b, em.isa, merge_engine.ErrorModelC(), merge_engine.ErrorModelChi2()); + logger.Info("Error model (XDS convention): a={:.3f} b={:.3e} c={:.3e} chi2={:.2f}", + em.a, em.b, ToXdsPartialityTerm(merge_engine.ErrorModelA(), merge_engine.ErrorModelC()), + merge_engine.ErrorModelChi2()); } merge_engine.AddImages(merge_input); if (config_.export_observations && !for_search && write_output) @@ -1828,6 +1829,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { result.error_model_isa_resolved = merge_engine.ErrorModelBResolved(); result.error_model_a = em.a; result.error_model_b = em.b; + result.error_model_c = ToXdsPartialityTerm(merge_engine.ErrorModelA(), merge_engine.ErrorModelC()); } out.cc_half_before_corrections = out.statistics.overall.cc_half; // no surfaces on this path logger.Info("Merge complete ({} unique reflections, {})", out.merged.size(), label); @@ -5602,7 +5604,8 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { result.error_model_isa_asymptotic > 0 && result.error_model_isa_resolved ? 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()}; + result.error_model_b > 0 ? fmt::format("{:.4e}", result.error_model_b) : std::string(), + result.error_model_c > 0 ? fmt::format("{:.4e}", result.error_model_c) : std::string()}; MtzTable merged_mtz = MergedMtzTable(sm.merged, *result.consensus_cell, experiment_, config_.provenance); // The unmerged MTZ reads the fulls until it is made. if (unmerged_mtz.valid()) diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index d551391bf..563252346 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -2503,6 +2503,7 @@ int RugnuxInputRun::PrepareRugnuxInput() { std::optional resolution_fit_A; // CC1/2 crossing, without the one-shell extension double error_model_a = 0.0; double error_model_b = 0.0; + double error_model_c = 0.0; // stills only std::optional reference_statistics; // --report-resolution double reference_isa = 0.0; @@ -2586,11 +2587,12 @@ int RugnuxInputRun::PrepareRugnuxInput() { merge_engine.RefineErrorModel(reflections); if (merge_engine.ErrorModelActive()) { - // 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. + // XDS convention, sigma^2 = a*(sigma0^2 + b*I^2 + c*((1-p)/p)^2*I^2). No ISa: on stills 1/b is + // the ceiling of one still observation, not what CORRECT.LP's ISa measures. const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB()); - logger.Info("Error model (XDS convention): a={:.3f} b={:.3e} ISa={:.1f} partiality c={:.3f} chi2={:.2f}", - em.a, em.b, em.isa, merge_engine.ErrorModelC(), merge_engine.ErrorModelChi2()); + logger.Info("Error model (XDS convention): a={:.3f} b={:.3e} c={:.3e} chi2={:.2f}", + em.a, em.b, ToXdsPartialityTerm(merge_engine.ErrorModelA(), merge_engine.ErrorModelC()), + merge_engine.ErrorModelChi2()); } merge_engine.AddImages(reflections); merged_reflections = merge_engine.ExportReflections(); @@ -2615,6 +2617,7 @@ int RugnuxInputRun::PrepareRugnuxInput() { error_model_isa_resolved = merge_engine.ErrorModelBResolved(); error_model_a = em.a; error_model_b = em.b; + error_model_c = ToXdsPartialityTerm(merge_engine.ErrorModelA(), merge_engine.ErrorModelC()); } } @@ -2804,7 +2807,8 @@ int RugnuxInputRun::PrepareRugnuxInput() { error_model_isa_asymptotic > 0 && error_model_isa_resolved ? 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()}; + error_model_b > 0 ? fmt::format("{:.4e}", error_model_b) : std::string(), + error_model_c > 0 ? fmt::format("{:.4e}", error_model_c) : std::string()}; WriteReflections(merged_reflections, *experiment.GetUnitCell(), experiment, merged_statistics, em_report, twinning, output_provenance, output_prefix, static_cast(nthreads), &scaled_fulls); @@ -2842,6 +2846,7 @@ int RugnuxInputRun::PrepareRugnuxInput() { scale_result.error_model_isa_resolved = error_model_isa_resolved; scale_result.error_model_a = error_model_a; scale_result.error_model_b = error_model_b; + scale_result.error_model_c = error_model_c; scale_result.resolution_fit_A = resolution_fit_A; scale_result.reference_merge_statistics = reference_statistics; scale_result.reference_isa = reference_isa; diff --git a/tests/ResultReportTest.cpp b/tests/ResultReportTest.cpp index d6ec42621..7646e9c0e 100644 --- a/tests/ResultReportTest.cpp +++ b/tests/ResultReportTest.cpp @@ -832,6 +832,32 @@ TEST_CASE("ResultReport_IsaUndetermined", "[Diagnostics]") { CHECK(text.find("ISA_ASYMPTOTIC") == std::string::npos); } +// A stills merge reports its error model as its three terms and no ISa anywhere. +TEST_CASE("ResultReport_StillsErrorModelTerms", "[Diagnostics]") { + DiffractionExperiment x(DetJF(1)); + ProcessResult result; + result.consensus_cell = UnitCell{.a = 79.0f, .b = 79.0f, .c = 38.0f, + .alpha = 90.0f, .beta = 90.0f, .gamma = 90.0f}; + result.has_merge_statistics = true; + result.merge_statistics.stills = true; + result.merge_statistics.overall.d_max = 50.0f; + result.merge_statistics.overall.d_min = 1.50f; + result.merge_statistics.overall.cc_half = 0.99; + result.error_model_a = 1.1; + result.error_model_b = 0.1; + result.error_model_c = 0.02; + result.error_model_isa = 3.0; + + const auto text = RenderResultReport("p", "in.h5", x, result); + CHECK(text.find("\nERROR_MODEL_B= 1.0000e-01\n") != std::string::npos); + CHECK(text.find("\nERROR_MODEL_C= 2.0000e-02\n") != std::string::npos); + CHECK(text.find("ISA=") == std::string::npos); + CHECK(text.find("ISa 3.0") == std::string::npos); + + result.error_model_isa_resolved = false; + CHECK(RenderResultReport("p", "in.h5", x, result).find("\nERROR_MODEL_B= undetermined\n") != std::string::npos); +} + // The completeness of the range where the signal is, beside the completeness of the range that was // written. Users compare the single overall figure against another program's and conclude the run // lost data, when it is the deliberately generous cut diluting the denominator - so the report says diff --git a/tools/battery/still/rugnux_stills_ab.py b/tools/battery/still/rugnux_stills_ab.py index 15976377b..816353853 100755 --- a/tools/battery/still/rugnux_stills_ab.py +++ b/tools/battery/still/rugnux_stills_ab.py @@ -201,7 +201,7 @@ def parse_cif(path): "isig": _f(val("_reflns.pdbx_netI_over_sigmaI")), "rmeas": _f(val("_reflns.pdbx_Rrim_I_all"), 100.0), "cc": _f(val("_reflns.pdbx_CC_half"), 100.0), - "isa": _f(val("_reflns.jfjoch_diffrn_ISa")), + "em_b": _f(val("_reflns.jfjoch_error_model_b")), # the stills error model writes no ISa "dmin": _f(val("_reflns.d_resolution_high")), "dmax": _f(val("_reflns.d_resolution_low")), } @@ -535,7 +535,7 @@ SUMMARY_COLS = [ ("", "isig", 8, 2), ("Rmeas%", "rmeas", 7, 2), ("CC1/2%", "cc", 7, 2), - ("ISa", "isa", 6, 2), + ("EM b", "em_b", 7, 4), ] diff --git a/viewer/windows/JFJochMergeStatsWindow.cpp b/viewer/windows/JFJochMergeStatsWindow.cpp index a49f009ac..24281bca5 100644 --- a/viewer/windows/JFJochMergeStatsWindow.cpp +++ b/viewer/windows/JFJochMergeStatsWindow.cpp @@ -149,7 +149,8 @@ JFJochMergeStatsWindow::JFJochMergeStatsWindow(const QString &title, const Merge // by the CC1/2 logistic cutoff, so it varies per run - show it first. if (o.d_min < o.d_max) hero->addWidget(MakeCard(QString::asprintf("%.2f", o.d_min), "Resolution [Å]", this)); - hero->addWidget(MakeCard(num(isa, 1), "ISa", this)); + if (!std::isnan(isa)) // NaN: a stills merge, which reports no ISa + hero->addWidget(MakeCard(num(isa, 1), "ISa", this)); hero->addWidget(MakeCard(pct(o.cc_half * 100.0), "CC½", this)); hero->addWidget(MakeCard(pct(o.r_meas * 100.0), "R-meas", this)); hero->addWidget(MakeCard(pct(completeness), "Completeness", this)); @@ -324,7 +325,7 @@ JFJochMergeStatsWindow::JFJochMergeStatsWindow(const QString &title, const Merge table_->verticalHeader()->setVisible(false); stack_->addWidget(chart_); // page 0 stack_->addWidget(table_); // page 1 - isa_page_ = buildIsaDiagnostic(merged_i_sigma, isa); + isa_page_ = std::isnan(isa) ? nullptr : buildIsaDiagnostic(merged_i_sigma, isa); if (isa_page_) stack_->addWidget(isa_page_); // page 2 else diff --git a/viewer/windows/JFJochMergeStatsWindow.h b/viewer/windows/JFJochMergeStatsWindow.h index 7c06c2c41..7cebb5536 100644 --- a/viewer/windows/JFJochMergeStatsWindow.h +++ b/viewer/windows/JFJochMergeStatsWindow.h @@ -31,7 +31,7 @@ class JFJochMergeStatsWindow : public QWidget { public: // `report` is the rugnux result report text (_report.txt), shown as its own view. JFJochMergeStatsWindow(const QString &title, const MergeStatistics &stats, - double isa, bool has_reference, + double isa, bool has_reference, // isa NaN: none (stills) const TwinningAnalysisResult &twinning, const std::optional &space_group, const std::optional &space_group_search, diff --git a/viewer/windows/JFJochProcessingJobsWindow.cpp b/viewer/windows/JFJochProcessingJobsWindow.cpp index 1a1a1262c..224d1dee1 100644 --- a/viewer/windows/JFJochProcessingJobsWindow.cpp +++ b/viewer/windows/JFJochProcessingJobsWindow.cpp @@ -9,6 +9,8 @@ #include "../../writer/RugnuxModelFile.h" #include "../../common/ThreadAffinity.h" +#include + #include #include #include @@ -797,7 +799,8 @@ void JFJochProcessingJobsWindow::onFinished(ProcessResult result, QString report jobs_[row].report = report; jobs_[row].checks = std::move(checks); jobs_[row].merge_stats = result.merge_statistics; - jobs_[row].isa = result.error_model_isa; + // A stills merge reports its error model as terms, not as an ISa (the report says why): NaN, no card. + jobs_[row].isa = result.merge_statistics.stills ? NAN : result.error_model_isa; jobs_[row].merge_has_reference = result.has_reference; jobs_[row].twinning = result.twinning; jobs_[row].space_group = result.space_group;