From aa87fb8bdefc73139bc7b2e6087e7d2d9e9bf4cf Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 19 Sep 2026 10:36:09 +0200 Subject: [PATCH] rugnux: print ISa as undetermined when the error model's b is not resolved from zero The (a, b) error-model fit already refuses to report ISa when its top intensity bin has no leverage on b. It still printed 1/b when b had leverage but came out at zero within its own error - the strongest bins scatter no more than counting statistics say (the fitted b^2 is negative and clamps to 0), or there are too few samples for the bin medians to mean anything. ISa then reports the noise in b as an I/sigma: on one low-resolution sweep consecutive merges of the same data gave ISa 0, 28, 46, 113 and 130. The standard error of b^2 is taken from the 16 bins' own scatter about the fitted line, and when b^2 is less than two standard errors above zero the merge result carries isa_resolved = false. Only what is printed follows it: the report's ISA key reads "undetermined", the summary line says so, the asymptote is not printed and the mmCIF carries "?". The fitted ISa itself is unchanged and is still what the space-group search's present-reflection cut and the refused-point-group arbitration read, so no decision moves. Co-Authored-By: Claude Opus 5 (1M context) --- .../scale_merge/RotationScaleMerge.cpp | 29 +++++++++++++++++-- .../scale_merge/RotationScaleMerge.h | 4 +++ rugnux/ResultReport.cpp | 12 ++++++-- rugnux/Rugnux.cpp | 12 ++++++-- rugnux/Rugnux.h | 1 + rugnux/rugnux_cli.cpp | 10 +++++-- tests/ResultReportTest.cpp | 24 +++++++++++++++ 7 files changed, 81 insertions(+), 11 deletions(-) diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 7457c5c05..bfde3cd48 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -3548,6 +3548,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool const std::vector *chi2_pool = nullptr; double chi2_a = 0.0, chi2_b2 = 0.0; bool error_model_b_unmeasured = false; // b had no leverage; ISa is not reported + bool error_model_b_resolved = true; // b^2 stood clear of its own standard error std::vector &samples = em_samples; samples.clear(); constexpr int n_bins = 16; @@ -3659,10 +3660,29 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool chi2_pool = ∈ chi2_a = error_model_a; chi2_b2 = 0.0; } } else if (det > 1e-10 * Ass * AII) { - error_model_a = std::clamp((Bs * AII - BI * AsI) / det, 0.25, 100.0); - const double b2 = std::max((Ass * BI - AsI * Bs) / det, 0.0); + const double a_fit = (Bs * AII - BI * AsI) / det; + const double b2_fit = (Ass * BI - AsI * Bs) / det; + error_model_a = std::clamp(a_fit, 0.25, 100.0); + const double b2 = std::max(b2_fit, 0.0); error_model_b = std::sqrt(b2); error_model_b_unmeasured = false; + // Leverage is not enough: b^2 can have it and still come out at zero within its own error - + // the strong bins scatter no more than counting says, or there are too few samples for a + // bin median to mean anything. 1/b then reads the noise in b as an I/sigma, any value from + // ~30 to infinity between merges of the same data. The standard error of b^2 is taken from + // the bins' own scatter about the fitted line; the fit and the sigmas are left as they are, + // and so is the fitted ISa every decision reads - this only decides what is printed. + // Measured: z = 13.6 and 22 on two crystals whose b is plainly visible, 0.1-0.5 on a + // 3000-sample low-multiplicity set and -1.2 to -6 on one whose strongest bins scatter + // LESS than counting says (both printed ISa 96-130). + double chi = 0.0; + for (int bin = 0; bin < n_bins; ++bin) { + const double d2w = std::max(bd2[bin], dev2_floor); + const double r = bd2[bin] - a_fit * bs2[bin] - b2_fit * bI2[bin]; + chi += r * r / (d2w * d2w); + } + const double b2_se = std::sqrt(chi / (n_bins - 2) * Ass / det); + error_model_b_resolved = b2_fit >= 2.0 * b2_se; error_model_active = true; chi2_pool = ∈ chi2_a = error_model_a; chi2_b2 = b2; } @@ -4069,6 +4089,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool { const auto em = ToXdsErrorModel(error_model_a, error_model_b); result.isa = em.isa; + result.isa_resolved = error_model_b_resolved; result.error_model_a = em.a; result.error_model_b = em.b; } @@ -4099,6 +4120,10 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool "not reach far enough for a systematic error to be seen, not that there is " "none; a resolution range matched to the signal would measure it", em.a, error_model_chi2); + else if (!error_model_b_resolved) + logger.Info("Error model (XDS convention): a={:.3f} b={:.3e} - b is not resolved from zero " + "(under two standard errors), so ISa is undetermined (1/b = {:.1f})", + em.a, em.b, em.isa); else if (!full_stats) // Neither the asymptote nor the chi2 median was measured on this merge, so neither is // reported for it. diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 3855f1a0f..533242ad3 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -49,6 +49,10 @@ public: // 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; + // b^2 stood at least two standard errors clear of zero. When it did not, `isa` is still the + // fitted 1/b and every decision reads it as before; only what is PRINTED says "undetermined", + // since 1/b of a b that is zero within its error is the noise in b, not an I/sigma. + bool isa_resolved = true; double error_model_a = 0.0; // XDS convention: sigma^2 = a*(sigma0^2 + b*I^2) double error_model_b = 0.0; // The overall CC1/2 as it stood BEFORE the correction surfaces were folded in. That is what the diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index f2f76ab82..cbb2ae8ba 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -695,8 +695,12 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, // 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}")); - Add(s, KeyReal("ISA", result.error_model_isa, "{:.2f}")); - if (result.error_model_isa_asymptotic > 0.0) + // A b that is zero within its error gives a 1/b that is the noise in b, not an I/sigma. + 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) @@ -1763,7 +1767,9 @@ 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 > 0.0) + 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)) diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 44a16b013..4d97df765 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -4916,6 +4916,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b /*measure_cc_before_corrections=*/true); result.error_model_isa = r.isa; result.error_model_isa_asymptotic = r.isa_asymptotic; + result.error_model_isa_resolved = r.isa_resolved; if (r.resolution_fit_A) result.resolution_fit_A = r.resolution_fit_A; result.error_model_a = r.error_model_a; result.error_model_b = r.error_model_b; @@ -6135,6 +6136,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // comparison merges are not the run's answer, so it goes back if they lose. const double em_isa = result.error_model_isa; const double em_isa_asymptotic = result.error_model_isa_asymptotic; + const bool em_isa_resolved = result.error_model_isa_resolved; const double em_a = result.error_model_a; const double em_b = result.error_model_b; const auto res_fit = result.resolution_fit_A; @@ -6201,6 +6203,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b experiment_.SetSpaceGroup(sg); result.error_model_isa = em_isa; result.error_model_isa_asymptotic = em_isa_asymptotic; + result.error_model_isa_resolved = em_isa_resolved; result.error_model_a = em_a; result.error_model_b = em_b; result.resolution_fit_A = res_fit; @@ -6697,9 +6700,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b && !superseded) { 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_isa > 0 && result.error_model_isa_resolved + ? fmt::format("{:.2f}", result.error_model_isa) : "?", + 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()}; WriteReflections(sm.merged, *result.consensus_cell, experiment_, sm.statistics, @@ -6776,6 +6780,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // multiplicity, so its fall-off is not this crystal's. const double em_isa = result.error_model_isa; const double em_isa_asymptotic = result.error_model_isa_asymptotic; + const bool em_isa_resolved = result.error_model_isa_resolved; const double em_a = result.error_model_a; const double em_b = result.error_model_b; const auto res_fit = result.resolution_fit_A; @@ -6825,6 +6830,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b experiment_.SetSpaceGroup(determined_group); result.error_model_isa = em_isa; result.error_model_isa_asymptotic = em_isa_asymptotic; + result.error_model_isa_resolved = em_isa_resolved; result.error_model_a = em_a; result.error_model_b = em_b; result.resolution_fit_A = res_fit; diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 2dc87f29f..d56907983 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -239,6 +239,7 @@ struct ProcessResult { MergeStatistics merge_statistics; 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 + 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; bool has_reference = false; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index ecd6e07cc..6ffbc469d 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1727,6 +1727,7 @@ static int RunRugnux(int argc, char **argv) { MergeStatistics merged_statistics; double error_model_isa = 0.0; double error_model_isa_asymptotic = 0.0; + bool error_model_isa_resolved = true; 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; @@ -1755,6 +1756,7 @@ static int RunRugnux(int argc, char **argv) { merged_statistics = std::move(r.statistics); error_model_isa = r.isa; error_model_isa_asymptotic = r.isa_asymptotic; + error_model_isa_resolved = r.isa_resolved; error_model_a = r.error_model_a; error_model_b = r.error_model_b; resolution_fit_A = r.resolution_fit_A; @@ -1952,9 +1954,10 @@ static int RunRugnux(int argc, char **argv) { 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_isa > 0 && error_model_isa_resolved + ? fmt::format("{:.2f}", error_model_isa) : "?", + 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()}; WriteReflections(merged_reflections, *experiment.GetUnitCell(), experiment, merged_statistics, @@ -1990,6 +1993,7 @@ static int RunRugnux(int argc, char **argv) { { std::ostringstream s; s << merged_statistics; scale_result.merge_statistics_text = s.str(); } scale_result.error_model_isa = error_model_isa; scale_result.error_model_isa_asymptotic = error_model_isa_asymptotic; + 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.resolution_fit_A = resolution_fit_A; diff --git a/tests/ResultReportTest.cpp b/tests/ResultReportTest.cpp index 9a246edb8..594306c74 100644 --- a/tests/ResultReportTest.cpp +++ b/tests/ResultReportTest.cpp @@ -651,6 +651,30 @@ TEST_CASE("ResultReport_FittedResolutionSuppressed", "[Diagnostics]") { CHECK(climbs_text.find("RESOLUTION_FIT") != std::string::npos); } +// ISa of a b that is zero within its error is the noise in b; the report says undetermined, not a number. +TEST_CASE("ResultReport_IsaUndetermined", "[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.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 = 1e-4; + result.error_model_isa = 95.0; + result.error_model_isa_asymptotic = 120.0; + + const auto resolved = RenderResultReport("p", "in.h5", x, result); + CHECK(resolved.find("\nISA= 95.00\n") != std::string::npos); + + result.error_model_isa_resolved = false; + const auto text = RenderResultReport("p", "in.h5", x, result); + CHECK(text.find("\nISA= undetermined\n") != std::string::npos); + CHECK(text.find("ISA_ASYMPTOTIC") == 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