rugnux: weight every CC1/2 by the information a reflection carries
The resolution cutoff, the shell table and the overall CC1/2 counted every unique reflection equally. The merge itself is inverse-variance weighted, so an observation from a frame the crystal barely diffracted on enters it at 1/G^2 of a good one - honestly, with its sigma - but a reflection measured only on such frames is scaled-up noise that then counts as much as a well-measured pair in every Pearson CC1/2 read off the merge. On a sweep where half the frames are weak the curve collapses at every resolution: the cut lands at 3.7 A on a P1 crystal whose good frames reach 1.5 A, and the UNUSABLE verdict fires (CC1/2 0.40 beside I/sigma 10). Each merged reflection now carries cc_weight: the precision its half-sets would have had with every observation at the run's typical frame scale, over the precision they have. G_ref = sum G^3 / sum G^2 over the usable observations is the precision-weighted typical scale, which the dead frames cannot drag down however many there are; the factor per observation is max(1, (G_ref/G)^2), with G the frame's total scale (partial scale, flux and the fulls' own G) taken before the correction surfaces and before collapsed frames are dropped, so nothing intensity- or resolution-dependent enters it. On a sweep without a weak stretch every weight is 1 and the CC1/2 is the plain Pearson it was. The cutoff fit, the shell table and the overall CC1/2 (and so the UNUSABLE verdict and the report's shell checks) all read the same weighted statistic. The two extra per-group sums are accumulated on both merge paths, the host loop and MergeAccumKernel, from one per-frame factor array; a host recompute of the device sums agrees to 1e-15 relative on every merge after the corrected corr is uploaded, and a host-merge run gives the same cut, space group, CC1/2 and ISa on five sets. Weighting by the half-set error variance alone (1/(v0+v1)) is not this: v grows with the intensity, so it weights the weak end of the intensity distribution and biases homogeneous data coarser. Measured (written resolution, together with the weighted outlier median): a P1 sweep with a long weak stretch 3.73 -> 1.40 A against a 1.63 A XDS reference, UNUSABLE withdrawn (overall CC1/2 0.40 -> 0.97); a second 3.61 -> 3.00 A; one with most of the sweep out of beam 7.35 -> 5.20 A. Lysozyme, thaumatin and two insulin sets unchanged to 0.01 A (one lysozyme sweep with a weak wedge 1.13 -> 1.16 A, from the median), same space groups throughout. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This commit is contained in:
@@ -306,14 +306,17 @@ namespace {
|
||||
// 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);
|
||||
double s_from, double s_to, int n, bool correlated, double sigma = 1.0) {
|
||||
std::normal_distribution<double> g(0.0, sigma);
|
||||
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);
|
||||
// The pair's weight in a CC1/2: the precision it would have had at the typical frame scale
|
||||
// over the precision it has - 1/sigma^2 for a band scaled up by sigma from dead frames.
|
||||
m.cc_weight = static_cast<float>(1.0 / (sigma * sigma));
|
||||
v.push_back(m);
|
||||
}
|
||||
}
|
||||
@@ -336,6 +339,54 @@ TEST_CASE("ResolutionCutoff_CleanFallOff") {
|
||||
CHECK(*rc.d_cut == Catch::Approx(1.92).margin(0.15));
|
||||
}
|
||||
|
||||
// A sweep with a long stretch where the crystal barely diffracted: the reflections measured only there
|
||||
// are noise scaled up by 1/G, sigmas with them, and there are many of them at every resolution. Counted
|
||||
// by the information they carry they must not hide where the well-measured reflections stop.
|
||||
TEST_CASE("ResolutionCutoff_ScaledUpNoiseDoesNotHideTheFallOff") {
|
||||
Logger logger("test");
|
||||
std::mt19937 rng(12345);
|
||||
std::vector<MergedReflection> merged;
|
||||
AddBand(merged, rng, 0.01, 0.25, 480, true);
|
||||
AddBand(merged, rng, 0.25, 0.51, 520, false);
|
||||
AddBand(merged, rng, 0.01, 0.51, 300, false, 50.0);
|
||||
|
||||
const auto rc = ComputeCCHalfLogisticCutoff(merged, 0.30, logger);
|
||||
REQUIRE(rc.d_fit);
|
||||
CHECK(*rc.d_fit == Catch::Approx(2.0).margin(0.15));
|
||||
|
||||
// Counted as equals, the same reflections hide it: the curve reads noise from the first bin on.
|
||||
for (auto &m : merged) m.cc_weight = 1.0f;
|
||||
const auto unweighted = ComputeCCHalfLogisticCutoff(merged, 0.30, logger);
|
||||
CHECK((!unweighted.d_fit || *unweighted.d_fit > 3.0));
|
||||
}
|
||||
|
||||
// The frame factor of the CC1/2 weight is 1 on a sweep whose frames all sit at one scale, however that
|
||||
// scale is spread over the frames, and 1 for any frame brighter than typical.
|
||||
TEST_CASE("CCHalfFrameFactors_UniformScaleIsOne") {
|
||||
const auto f = CCHalfFrameFactors({1.7, 1.7, 1.7, 1.7}, {100, 3, 250, 0});
|
||||
for (double x : f) CHECK(x == 1.0);
|
||||
const auto g = CCHalfFrameFactors({1.0, 2.0}, {100, 100});
|
||||
CHECK(g[1] == 1.0);
|
||||
CHECK(g[0] == Catch::Approx(9.0 / 5.0 * 9.0 / 5.0)); // G_ref = (1 + 8) / (1 + 4)
|
||||
}
|
||||
|
||||
// Most of the sweep at 2% of the good frames' scale: the typical scale is still the good frames', not
|
||||
// the run median, and an observation from a dead frame carries 2500x the variance of a good one.
|
||||
TEST_CASE("CCHalfFrameFactors_DeadStretchDoesNotSetTheTypicalScale") {
|
||||
std::vector<double> scale;
|
||||
std::vector<int64_t> n;
|
||||
for (int f = 0; f < 100; ++f) {
|
||||
scale.push_back(f < 60 ? 0.02 : 1.0);
|
||||
n.push_back(50);
|
||||
}
|
||||
scale.push_back(NAN); // a frame without a scale counts as G = 1
|
||||
n.push_back(50);
|
||||
const auto factor = CCHalfFrameFactors(scale, n);
|
||||
CHECK(factor[0] == Catch::Approx(2500.0).epsilon(0.01));
|
||||
CHECK(factor[99] == Catch::Approx(1.0).epsilon(1e-3));
|
||||
CHECK(factor[100] == Catch::Approx(1.0).epsilon(1e-3));
|
||||
}
|
||||
|
||||
// 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
|
||||
|
||||
Reference in New Issue
Block a user