From d00175c8b478ea5dba8bbd2ef76a609ddc62628b Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 14 Sep 2026 14:17:39 +0200 Subject: [PATCH] resolution: a cut the shell table refutes is not quoted One run was written to 1.089 A where I/sigma was zero and R_meas 1776 per cent, and the summary said only "Merged to 1.09 A". Two things let that through. The logistic's crossing was taken wherever the fitted curve met the target, even when that point lay past every bin the curve was fitted through - an extrapolation of a fall-off the data never showed, quoted as a measurement. The crossing now has to lie inside the fitted bins, with the one-shell extension still there to spare; outside them the number is read off the bins themselves instead. And the shipped guard asked only for the finest shell whose CC1/2 still reached the target, with no requirement that the curve get there monotonically. A noise shell that climbs back over the bar therefore became the edge of the data. It is the climbing back that disqualifies it, and the program already noticed - it printed "CC1/2 is not monotone with resolution" and then cut there anyway, because the test was developer-only and decided nothing. It decides now, in the report everyone reads, and the duplicate is gone so there is one such test rather than two. The dataset above is written to 1.527 A: CC1/2 0.209 to 0.688, I/sigma 0.40 to 1.26, R_meas 245 to 145 per cent. Forty-four of fifty-one reports are unchanged to the character; of the seven that move, five are cut coarser and every headline number of all five improves, one loses a fit that was an extrapolation without changing anything written, and one gains a quotable fit and loses a warning. This is a guard, not the cause. On four of those five the reflections doing the damage are ice, which the resolution fit now leaves out for its own reasons; the guard still earns its place, because on those four removing the ice alone does not stop the cut being quoted past what the shells support. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_011GxZqDiFP3KqriBhNdcR56 --- docs/CHANGELOG.md | 1 + .../scale_merge/ResolutionCutoff.cpp | 31 ++++++-- image_analysis/scale_merge/ResolutionCutoff.h | 7 +- rugnux/ResultReport.cpp | 70 ++++++++++--------- tests/MergeScaleTest.cpp | 59 ++++++++++++++++ tests/ResultReportTest.cpp | 15 ++++ 6 files changed, 145 insertions(+), 38 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 7738baad8..0167b7dcb 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -4,6 +4,7 @@ ### 1.0.0-rc.169 * Building Jungfraujoch no longer needs zlib or Eigen installed on the machine: the build downloads and builds both, like every other dependency. +* rugnux no longer writes a run to a resolution its own shell table refutes: a CC1/2 curve that drops below the cut's target and climbs back above it is not read as a fall-off, and the report says so in place of quoting FITTED_RESOLUTION. * rugnux drops the stretches of a rotation sweep whose removal measurably improves the merged intensities (delta-CC1/2), and reports what became of every frame - merged, downgraded or rejected, in frames and in degrees. * The viewer draws its spot markers over a black outline, so they stay visible on a light colour map and are no longer mistaken for a magenta bad pixel or the coral beam stop. * The viewer has a background (black-point) slider beside the foreground one; hold B and scroll the wheel to adjust it. diff --git a/image_analysis/scale_merge/ResolutionCutoff.cpp b/image_analysis/scale_merge/ResolutionCutoff.cpp index 2d1005734..071a32b56 100644 --- a/image_analysis/scale_merge/ResolutionCutoff.cpp +++ b/image_analysis/scale_merge/ResolutionCutoff.cpp @@ -107,9 +107,11 @@ ResolutionCutoffResult ComputeCCHalfLogisticCutoff(const std::vector the fall-off crossing), NOT the full measured range: when the detector reaches far // past where the crystal diffracts (a high-res-configured detector on a low-res crystal), the // full range is dominated by high-res noise, so s_range/10 would be a huge over-extension. - const double s_cross = s0 + std::log(1.0 / cc_target - 1.0) / k; + double s_cross = s0 + std::log(1.0 / cc_target - 1.0) / k; + + // The crossing has to be one the bins show, not one the logistic extrapolates to. On a real + // fall-off it sits between the last bin at the target and the first bin below it, with the two + // extension bins to spare. A crossing past the last fitted bin means the logistic never followed + // these data down - it settled on a shallow slope through a fall-off region that is not one - + // and every resolution read off it is one the bins refute. The bins' own crossing is the honest + // answer then: where the measured CC1/2 passes the target, between the two bins that straddle it. + bool logistic = true; + if (s_cross > s_bin.back()) { + if (first_below >= cc_bin.size()) { + result.note = "CC1/2 does not fall off within the measured range"; + return result; + } + logistic = false; + const double f = (cc_bin[first_below - 1] - cc_target) + / (cc_bin[first_below - 1] - cc_bin[first_below]); + s_cross = s_bin[first_below - 1] + f * (s_bin[first_below] - s_bin[first_below - 1]); + } + const double delta_s = (s_cross - s_lo) / SHELLS_FOR_EXTENSION; const double s_final = s_cross + delta_s; @@ -194,7 +215,9 @@ ResolutionCutoffResult ComputeCCHalfLogisticCutoff(const std::vector ApplyResolutionCutoff(std::vector &merge if (rc.d_cut) { effective_d_min = rc.d_cut; if (rc.d_fit) - logger.Info("Auto resolution cutoff: {:.2f} A (CC1/2 fit crosses {:.2f} at {:.2f} A, " + logger.Info("Auto resolution cutoff: {:.2f} A ({} crosses {:.2f} at {:.2f} A, " "written one shell past it; override with --scaling-high-resolution)", - *rc.d_cut, cc_target, *rc.d_fit); + *rc.d_cut, rc.note, cc_target, *rc.d_fit); else logger.Info("Auto resolution cutoff: {:.2f} A ({}; override with --scaling-high-resolution)", *rc.d_cut, rc.note); diff --git a/image_analysis/scale_merge/ResolutionCutoff.h b/image_analysis/scale_merge/ResolutionCutoff.h index 229ac21b0..b6913aea3 100644 --- a/image_analysis/scale_merge/ResolutionCutoff.h +++ b/image_analysis/scale_merge/ResolutionCutoff.h @@ -18,12 +18,15 @@ // Method (see docs/rugnux_resolution_cutoff_design.md): bin CC1/2 against s = 1/d^2 in fine bins, // fit a logistic CC1/2(s) = 1/(1+exp(k*(s-s0))) to the contiguous-from-low-res fall-off, take the s // where the fit crosses cc_target, then extend by one mean (10-shell) shell width in s ("one shell -// too far", generous). d_cut is nullopt when the fit is degenerate (too few bins, flat/non-monotone) +// too far", generous). Where the fitted crossing lands past the bins the fit was made over - a +// fall-off region too ragged for a logistic to follow, so the crossing is an extrapolation rather +// than something the bins show - the crossing is read off the bins themselves instead, and +// everything after it is unchanged. d_cut is nullopt when the fit is degenerate (too few bins, flat) // or CC1/2 never falls below cc_target inside the measured range - the caller then keeps the full // range. cc_target in (0,1); merged must carry finite I_half[0]/I_half[1] to contribute. struct ResolutionCutoffResult { std::optional d_cut; // high-resolution limit (A); nullopt => keep the full range - // Where the fitted CC1/2 crosses cc_target, BEFORE the deliberate one-shell extension - i.e. the + // Where CC1/2 crosses cc_target, BEFORE the deliberate one-shell extension - i.e. the // resolution the data are judged to reach, as opposed to the (coarser in s, finer in d) limit the // reflections are actually written to. This is the number to quote. Set whenever the fit produced // a crossing inside the measured range, even when no cut was applied. diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 7137f885d..92e17e35c 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -191,6 +191,25 @@ namespace { finest = sh.d_min; return finest; // 0.0 = no shell reaches the target at all } + + // Where CC1/2 climbs back above the cut's target after it has already fallen below it. Below the + // target a shell carries no signal the cut would keep, so a shell past that point which + // correlates again is noise correlating with noise: the curve is not a fall-off, and no single + // resolution can be read off it. Returns that shell's d_min, 0.0 when the curve never climbs + // back. A wobble that stays ABOVE the target is not this - a fall-off with a noisy shell in it + // is still a fall-off. + double CCHalfClimbsBackAboveTarget(const MergeStatistics &ms, double target) { + bool below = false; + for (const auto &sh : ms.shells) { + if (!std::isfinite(sh.cc_half)) + continue; + if (sh.cc_half < target) + below = true; + else if (below) + return sh.d_min; + } + return 0.0; + } } // ---------------------------------------------------------------------------- the emitter @@ -560,21 +579,29 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, // than the finest shell that still reaches the target asserts a resolution the shell // table on this very page refutes. const double finest_supported = FinestSupportedShell(ms, cc_target); + const double climbs_back = CCHalfClimbsBackAboveTarget(ms, cc_target); std::string fit_suppressed; // Only where there are shells to judge the fit against; with no table there is no // evidence either way, and suppressing on no evidence is its own false claim. - if (result.resolution_fit_A && !ms.shells.empty()) { - if (finest_supported == 0.0) + if (!ms.shells.empty()) { + if (climbs_back > 0.0) fit_suppressed = fmt::format( - "CC1/2 never reaches the cut's target of {:.2f} in any shell, so there is no" - " fall-off to fit", cc_target); - else if (*result.resolution_fit_A < finest_supported - 1e-6) - fit_suppressed = fmt::format( - "the fit reads {:.2f} A, finer than the finest shell whose CC1/2 still reaches" - " {:.2f} ({:.2f} A), so the curve it was read off is not a fall-off", - *result.resolution_fit_A, cc_target, finest_supported); - else - fitted_resolution = *result.resolution_fit_A; + "CC1/2 falls below the cut's target of {:.2f} and then climbs back above it" + " at {:.2f} A, so this curve is not a fall-off and no one resolution can be" + " read off it", cc_target, climbs_back); + else if (result.resolution_fit_A) { + if (finest_supported == 0.0) + fit_suppressed = fmt::format( + "CC1/2 never reaches the cut's target of {:.2f} in any shell, so there is" + " no fall-off to fit", cc_target); + else if (*result.resolution_fit_A < finest_supported - 1e-6) + fit_suppressed = fmt::format( + "the fit reads {:.2f} A, finer than the finest shell whose CC1/2 still" + " reaches {:.2f} ({:.2f} A), so the curve it was read off is not a" + " fall-off", *result.resolution_fit_A, cc_target, finest_supported); + else + fitted_resolution = *result.resolution_fit_A; + } } if (fitted_resolution) Add(s, KeyReal("FITTED_RESOLUTION", *fitted_resolution, "{:.2f}")); @@ -1289,27 +1316,6 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, "interpretation should allow for it; no intensity has been corrected for it here", an.delta_b)); } - - // CC1/2 that falls and then climbs again is not a fall-off. An observation about the - // resolution cut's own internals, which no user can act on: it belongs in the developer - // report, because a warning nobody can act on teaches people to skip the section. - { - // A climb only counts while the shell still carries signal: past CC1/2 ~ 0.3 the data - // are beyond any cut this run would take and the number is oscillating in noise. - constexpr double CC_HALF_RISE = 0.05, CC_HALF_ALIVE = 0.30; - const auto &sh = result.merge_statistics.shells; - size_t fell = 0, rose_after = 0; - for (size_t i = 1; i < sh.size(); ++i) { - if (sh[i].cc_half < sh[i - 1].cc_half) fell = i; - else if (fell > 0 && sh[i].cc_half > sh[i - 1].cc_half + CC_HALF_RISE - && sh[i].cc_half > CC_HALF_ALIVE) rose_after = i; - } - if (rose_after > 0) - Warn(doc, PathologyCode::RESOLUTION_FIT, fmt::format( - "CC1/2 is not monotone with resolution (it climbs again at {:.2f}-{:.2f} A) - the " - "fall-off the resolution cut is read off does not describe these data", - sh[rose_after].d_max, sh[rose_after].d_min), true); - } } doc.sections.push_back(std::move(s)); } diff --git a/tests/MergeScaleTest.cpp b/tests/MergeScaleTest.cpp index e5b90857e..5ca2294cb 100644 --- a/tests/MergeScaleTest.cpp +++ b/tests/MergeScaleTest.cpp @@ -2,8 +2,11 @@ // SPDX-License-Identifier: GPL-3.0-only #include + +#include #include "../image_analysis/scale_merge/HKLKey.h" #include "../image_analysis/scale_merge/Merge.h" +#include "../image_analysis/scale_merge/ResolutionCutoff.h" #include "gemmi/reciproc.hpp" TEST_CASE("HKLKey_NoSG_noMergeFriedel") { @@ -255,3 +258,59 @@ TEST_CASE("MergeStats_NoReferenceCellLeavesCompletenessUnmeasured") { CHECK(Completeness(stats.overall) == 0.0); } } + +// ---------------------------------------------------------------- the automatic resolution cutoff +namespace { + // Half-set pairs spread uniformly in s = 1/d^2 over [s_from, s_to), either correlated with each + // other (signal) or drawn independently (noise, CC1/2 ~ 0), so a whole CC1/2 curve can be built + // band by band. + void AddBand(std::vector &v, std::mt19937 &rng, + double s_from, double s_to, int n, bool correlated) { + std::normal_distribution g(0.0, 1.0); + for (int j = 0; j < n; ++j) { + MergedReflection m; + m.d = static_cast(1.0 / std::sqrt(s_from + (j + 0.5) * (s_to - s_from) / n)); + const double a = g(rng), b = g(rng); + m.I_half[0] = static_cast(a); + m.I_half[1] = static_cast(correlated ? a : b); + v.push_back(m); + } + } +} + +// A clean fall-off: CC1/2 crosses the target where the signal stops, and the cut is written one +// shell past it. +TEST_CASE("ResolutionCutoff_CleanFallOff") { + Logger logger("test"); + std::mt19937 rng(12345); + std::vector merged; + AddBand(merged, rng, 0.01, 0.25, 480, true); // signal to 1/sqrt(0.25) = 2.00 A + AddBand(merged, rng, 0.25, 0.51, 520, false); // noise beyond it + + const auto rc = ComputeCCHalfLogisticCutoff(merged, 0.30, logger); + REQUIRE(rc.d_fit); + CHECK(*rc.d_fit == Catch::Approx(2.0).margin(0.15)); + REQUIRE(rc.d_cut); + CHECK(*rc.d_cut < *rc.d_fit); // the deliberate one-shell extension + CHECK(*rc.d_cut == Catch::Approx(1.92).margin(0.15)); +} + +// A fall-off region a logistic cannot follow: CC1/2 drops through the target and comes straight back +// up. The fitted crossing is then an extrapolation far past the bins it was made over, and reading +// the cut off it writes the data deep into the noise; the crossing the bins themselves show is where +// the signal stopped, and that is what must be used. +TEST_CASE("ResolutionCutoff_RaggedFallOffIsReadOffTheBins") { + Logger logger("test"); + std::mt19937 rng(12345); + std::vector merged; + AddBand(merged, rng, 0.01, 0.13, 240, true); // signal to 1/sqrt(0.13) = 2.77 A + AddBand(merged, rng, 0.13, 0.17, 80, false); // a hole below the target + AddBand(merged, rng, 0.17, 0.25, 160, true); // correlated again - not a fall-off + AddBand(merged, rng, 0.25, 0.51, 520, false); + + const auto rc = ComputeCCHalfLogisticCutoff(merged, 0.30, logger); + REQUIRE(rc.d_fit); + CHECK(*rc.d_fit == Catch::Approx(2.77).margin(0.20)); + REQUIRE(rc.d_cut); + CHECK(*rc.d_cut > 2.30); // coarser than the band that correlates again +} diff --git a/tests/ResultReportTest.cpp b/tests/ResultReportTest.cpp index ec02e1b5d..9dcec514a 100644 --- a/tests/ResultReportTest.cpp +++ b/tests/ResultReportTest.cpp @@ -572,6 +572,21 @@ TEST_CASE("ResultReport_FittedResolutionSuppressed", "[Diagnostics]") { result.resolution_fit_A = 0.95; CHECK(RenderResultReport("p", "in.h5", x, result).find("\nFITTED_RESOLUTION= 0.95\n") != std::string::npos); + + // CC1/2 that drops through the target and climbs back above it is not a fall-off either, and the + // finest shell reaching the target is then noise correlating with noise - so a fit the test above + // accepts must still be suppressed, and the user told so in the default report. + ProcessResult climbs = result; + climbs.merge_statistics.shells.back().cc_half = 0.05; + MergeStatisticsShell noise = climbs.merge_statistics.shells.back(); + noise.d_max = 0.80f; + noise.d_min = 0.70f; + noise.cc_half = 0.60; + climbs.merge_statistics.shells.push_back(noise); + const auto climbs_text = RenderResultReport("p", "in.h5", x, climbs); + CHECK(climbs_text.find("FITTED_RESOLUTION=") == std::string::npos); + CHECK(climbs_text.find("climbs back above it at 0.70 A") != std::string::npos); + CHECK(climbs_text.find("RESOLUTION_FIT") != std::string::npos); } // The completeness of the range where the signal is, beside the completeness of the range that was