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