mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 67 Wiki Activity

Merging: do not floor the merged sigma at the systematic term

Browse Source 
The merged sigma was floored at b*|I|, so I/sigma could never exceed the reported ISa.
On one dataset every merged reflection came out at I/sigma <= 12.96 with a 99th
percentile of 12.77 in every resolution shell alike, while the scatter of the
observations implied about 44 and XDS reported 58.

The floor is wrong in principle. `b` is fitted from the scatter BETWEEN a reflection's
symmetry equivalents, i.e. from the part that is not common to them, so it averages
down with multiplicity exactly like the counting term. 1/sqrt(sum_w) with the
b-inflated per-observation sigma already gives b*I/sqrt(n); flooring at b*|I| puts the
sqrt(n) back. That is the whole effect: 12.96 * sqrt(21.6) = 60, against XDS's 58.

It was introduced on a comparison of our MERGED I/sigma against XDS's UNMERGED
I/sigma. XDS's own merged low-resolution I/sigma exceeds its reported ISa on 30 of the
39 reference datasets here, median ratio 1.78 and up to 4.23.

Merged low-shell I/sigma now lands where XDS's does: 22.4 -> 46.2 against 46.2 on one
crystal, 26.7 -> 115.7 against 96.6 on another, 12.5 -> 45.0 against 58.0 on a third.
Over the 38-crystal battery the space groups, the merged reflection sets, R_meas and
CC1/2 are all unchanged - every one of them is sigma-independent, which is what makes
them the right control - and <I/sigma> rises on 35 crystals with none worse.

The asymptotic estimator that fed the floor stays, for the reported ISa only, and is
repaired in the process: it subtracts a*sigma^2 rather than the raw sigma^2 (at a < 1
the difference is the same size as the b^2 being measured, which is what made it
flip between 10.9 and 62.7 on consecutive passes of the same data), it rescales each
group's variance median-unbiased before subtracting an unbiased counting term, its
I/sigma gate uses the same convention, and it is bounded by the whole-range b - an
asymptote exists to refine 1/b upward, not to report 0.3 because "strong" was selected
on a sigma scale the fit itself rejects.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-08-09 20:58:35 +02:00
co-authored by Claude Opus 5
parent f60768d49c
commit 72efb75a8c
8 changed files with 82 additions and 87 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/RotationScaleMerge.cpp
+47 -33
View File
@@ -1801,10 +1801,8 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
// intensity raises its systematic term `b`, so 1/b understates that limit. Read it instead directly
// from the strong equivalents: for each well-measured reflection group the counting-subtracted
// fractional scatter of its symmetry mates estimates the systematic term, and the robust median over
// strong groups is the asymptote. The I/sigma threshold that selects "strong" is relaxed on weak /
// radiation-damaged data that has few strong reflections (rather than fall back to the higher
// whole-range b); it falls back only when even the relaxed set is too small. Host-side over the
// merged fulls, so CPU and GPU agree.
// strong groups is the asymptote. Report-only: it is the reported ISa, nothing downstream uses it.
// Host-side over the merged fulls, so CPU and GPU agree.
double error_model_b_asymptotic = error_model_b;
if (error_model_active) {
struct GroupScatter { double sum = 0, sum_sq = 0, sum_var = 0; int n = 0; };
@@ -1816,47 +1814,60 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
auto &g = gs[o.group];
g.sum += I_corr; g.sum_sq += I_corr * I_corr; g.sum_var += sigma_corr * sigma_corr; ++g.n;
}
// The counting term to subtract. The estimator below has always used the RAW reported
// sigma^2, but the (a, b) fit immediately above concluded the counting variance is
// a*sigma^2. Where a < 1 the raw subtraction overshoots: the per-group systematic is driven
// to zero, the median lands on the boundary, and 1/b_asy reports an impossible I/sigma
// (measured: one dataset reporting 64.6 where every other statistic supports ~16, and
// flipping between 10.9 and 64.6 on consecutive merges of the same data). Opt-in for now
// (env JFJOCH_ISA_ASY_A) so the default path stays bit-identical while the two are
// batteried against each other -- an earlier attempt at this subtraction was rejected for
// over-claiming on a > 1 data, which this change would reintroduce.
const double asy_counting_scale =
std::getenv("JFJOCH_ISA_ASY_A") != nullptr ? error_model_a : 1.0;
// Per-group counting-subtracted fractional systematic variance, paired with the group's I/sigma.
// Two conventions decide whether this measures anything:
// * the counting term is the error model's OWN, a*sigma^2. Subtracting the raw sigma^2 while the
// fit has just concluded the counting variance is a*sigma^2 leaves a bias of (1-a)*sigma^2,
// which at the I/sigma admitted here is the same size as the systematic being measured - so at
// a < 1 the per-group value collapses onto zero and 1/b_asy reports an impossible I/sigma.
// The I/sigma that selects "strong" uses the same counting term, so the gate and the
// subtraction agree about what a strong reflection is.
// * the sample variance of n observations is chi^2_(n-1)-distributed, and its MEDIAN lies below
// its mean (16% at n = 5, 5% at n = 13). Taking a median across groups of a variance minus an
// unbiased counting term therefore subtracts more than it should, again by an amount
// comparable to the systematic. Rescale each group's variance to be median-unbiased first
// (Wilson-Hilferty median(chi^2_k) = k*(1-2/9k)^3, exact to 0.4% for k >= 4) - the same
// median-to-mean conversion the (a, b) fit does with CHI2_1_MEDIAN.
// The clamp at zero is gone with them: median(max(x,0)) = max(median(x),0), so it never moved a
// positive median, and a negative median is the informative answer "not measurable here".
std::vector<std::pair<double, double>> group_scatter; // (systematic b^2, I/sigma)
for (const auto &g : gs) {
if (g.n < 5) continue;
const double mean = g.sum / g.n;
const double counting = asy_counting_scale * (g.sum_var / g.n);
const double counting = error_model_a * (g.sum_var / g.n);
if (mean <= 0.0 || counting <= 0.0) continue;
const double variance = (g.sum_sq - g.sum * g.sum / g.n) / (g.n - 1);
group_scatter.push_back({std::max(variance - counting, 0.0) / (mean * mean),
const double k = g.n - 1;
const double median_of_chi2 = std::pow(1.0 - 2.0 / (9.0 * k), 3);
const double variance = (g.sum_sq - g.sum * g.sum / g.n) / k / median_of_chi2;
group_scatter.push_back({(variance - counting) / (mean * mean),
mean / std::sqrt(counting)});
}
auto asymptote_above = [&](double snr_min, size_t min_groups) -> double {
// The threshold is relaxed on weak / radiation-damaged data that has too few strong reflections
// for the tight one. A tier that HAS enough groups gives its answer and is not retried lower
// because that answer came out small - the retry is what made the report flip between 1/b and an
// absurd value on consecutive merges of statistically identical data.
auto asymptote_above = [&](double snr_min, size_t min_groups) -> std::optional<double> {
std::vector<double> b2;
for (const auto &[b2_value, snr] : group_scatter)
if (snr >= snr_min) b2.push_back(b2_value);
return b2.size() >= min_groups ? std::sqrt(median_of(b2)) : 0.0;
return b2.size() >= min_groups ? std::optional(median_of(b2)) : std::nullopt;
};
double b_asy = asymptote_above(20.0, 100); // tight threshold on data that supports it
if (b_asy <= 0.0) b_asy = asymptote_above(10.0, 50); // relaxed for weak / damaged data
if (b_asy > 0.0) error_model_b_asymptotic = b_asy;
auto b2_asy = asymptote_above(20.0, 100); // tight threshold on data that supports it
if (!b2_asy) b2_asy = asymptote_above(10.0, 50); // relaxed for weak / damaged data
// A non-positive median means the strong equivalents reproduce each other to within counting
// statistics: the asymptote is below what this data can resolve, so report the whole-range b
// rather than an extreme extrapolated from noise. The asymptote can also only ever REFINE 1/b
// upwards - that is the whole reason it is measured - so a strong-group estimate that comes out
// WORSE than the fit's own b has not measured an asymptote at all: it means "strong" was
// selected on a sigma scale that the fit itself rejects, which is what happens on data too weak
// to have strong reflections. Inert on healthy data (b_asy sits 0-40% below b there).
if (b2_asy && *b2_asy > 0.0)
error_model_b_asymptotic = std::min(std::sqrt(*b2_asy), error_model_b);
}
// Guard a degenerate low-multiplicity fit: with too few symmetry equivalents both the (a, b) fit and
// the per-group scatter collapse toward zero, and 1/b then reports an impossibly high asymptotic
// I/sigma. Real macromolecular data does not exceed ISa ~50; past a generous cap report the asymptote
// as unmeasured rather than emit a spurious extreme. Only the REPORT is dropped: the sigma floor
// below keeps using the fitted value, because a degenerate fit is precisely where the raw
// 1/sqrt(sum_w) is least trustworthy and leaving the merged sigma uncapped would be the opposite of
// what this guard is for.
// as unmeasured rather than emit a spurious extreme.
constexpr double MIN_ASYMPTOTIC_B = 0.01; // ISa cap 100
const double isa_reported = error_model_b_asymptotic >= MIN_ASYMPTOTIC_B
? 1.0 / error_model_b_asymptotic : 0.0;
@@ -1929,16 +1940,19 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
MergedReflection mr{};
mr.h = group_h[g]; mr.k = group_k[g]; mr.l = group_l[g];
mr.I = static_cast<float>(a.swI / a.sw);
// The systematic floor caps a high-multiplicity merged sigma at the reproducibility that the
// strongest reflections actually reach, i.e. the asymptotic term, so merged I/sigma approaches ISa.
mr.sigma = SigmaWithSystematicFloor(1.0 / std::sqrt(a.sw), mr.I, error_model_b_asymptotic);
// Plain inverse-variance merged sigma. The error model's systematic term (b*I)^2 is measured
// from the scatter BETWEEN a reflection's symmetry equivalents, i.e. from exactly the part of
// the error that is NOT common to them, so it averages down over the multiplicity like the
// counting part and the merge must not hold it back. XDS behaves the same way: its merged
// I/sigma runs far above its own reported ISa in the strong low-resolution shells.
mr.sigma = static_cast<float>(1.0 / std::sqrt(a.sw));
mr.I_half[0] = mr.I_half[1] = NAN;
mr.sigma_half[0] = mr.sigma_half[1] = NAN;
mr.d = a.d;
if (a.nh[0] + a.nh[1] > 0 && a.swh[0] > 0.0 && a.swh[1] > 0.0) {
for (int i = 0; i < 2; ++i) {
mr.I_half[i] = static_cast<float>(a.swIh[i] / a.swh[i]);
mr.sigma_half[i] = SigmaWithSystematicFloor(1.0 / std::sqrt(a.swh[i]), mr.I_half[i], error_model_b_asymptotic);
mr.sigma_half[i] = static_cast<float>(1.0 / std::sqrt(a.swh[i]));
}
}
if (!std::isfinite(a.d) || a.d <= 0.0f) continue;