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