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) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_011GxZqDiFP3KqriBhNdcR56
This commit is contained in:
@@ -2,8 +2,11 @@
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#include <catch2/catch_all.hpp>
|
||||
|
||||
#include <random>
|
||||
#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<MergedReflection> &v, std::mt19937 &rng,
|
||||
double s_from, double s_to, int n, bool correlated) {
|
||||
std::normal_distribution<double> g(0.0, 1.0);
|
||||
for (int j = 0; j < n; ++j) {
|
||||
MergedReflection m;
|
||||
m.d = static_cast<float>(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<float>(a);
|
||||
m.I_half[1] = static_cast<float>(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<MergedReflection> 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<MergedReflection> 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
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user