From 72efb75a8cf790dacec15f528fec397e2ec5f5aa Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 9 Aug 2026 20:58:35 +0200 Subject: [PATCH 1/2] Merging: do not floor the merged sigma at the systematic term 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 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) --- image_analysis/scale_merge/Merge.cpp | 6 +- image_analysis/scale_merge/Merge.h | 12 --- .../scale_merge/RotationScaleMerge.cpp | 80 +++++++++++-------- .../scale_merge/SearchSpaceGroup.cpp | 30 +++---- image_analysis/scale_merge/SearchSpaceGroup.h | 4 +- rugnux/rugnux_cli.cpp | 2 +- tests/SearchSpaceGroupTwinTest.cpp | 8 +- tests/SyntheticMergedReflections.h | 27 +++---- 8 files changed, 82 insertions(+), 87 deletions(-) diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 0db26709..7ad95166 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -358,8 +358,7 @@ std::vector MergeOnTheFly::ExportReflections() { .k = accum.k, .l = accum.l, .I = static_cast(accum.sum_wI / accum.sum_w), - .sigma = SigmaWithSystematicFloor(1.0 / std::sqrt(accum.sum_w), - static_cast(accum.sum_wI / accum.sum_w), error_model_b), + .sigma = static_cast(1.0 / std::sqrt(accum.sum_w)), .I_half = {NAN, NAN}, .sigma_half = {NAN, NAN}, .d = accum.d @@ -368,8 +367,7 @@ std::vector MergeOnTheFly::ExportReflections() { if (accum.n_half[0] + accum.n_half[1] > 0 && accum.sum_w_half[0] > 0.0 && accum.sum_w_half[1] > 0.0) { for (int i = 0; i < 2; ++i) { mr.I_half[i] = static_cast(accum.sum_wI_half[i] / accum.sum_w_half[i]); - mr.sigma_half[i] = SigmaWithSystematicFloor(1.0 / std::sqrt(accum.sum_w_half[i]), - mr.I_half[i], error_model_b); + mr.sigma_half[i] = static_cast(1.0 / std::sqrt(accum.sum_w_half[i])); } } diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 464ac152..dc250a15 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -16,18 +16,6 @@ #include "HKLKey.h" -// The error model splits a reflection's variance into a statistical part (a*sigma^2, which averages -// down with multiplicity) and a systematic part ((b*I)^2 - absorption, beam flicker, partiality, -// detector non-uniformity - correlated across a reflection's repeats). Inverse-variance merging -// (sigma = 1/sqrt(sum_w)) wrongly divides BOTH by the multiplicity, so high-multiplicity reflections -// get an unphysically small merged sigma (merged I/sigma far above ISa). Floor the merged sigma at -// b*|I| so the systematic term survives the merge and ISa = 1/b stays the asymptotic I/sigma ceiling. -// error_model_b <= 0 (no active error model) leaves the sigma unchanged. -inline float SigmaWithSystematicFloor(double inv_variance_sigma, float merged_I, double error_model_b) { - const auto floor = static_cast(error_model_b * std::abs(static_cast(merged_I))); - return std::max(static_cast(inv_variance_sigma), floor); -} - struct MergeStatisticsShell { float d_min = 0.0f; float d_max = 0.0f; diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index d73ace1f..9e5d055d 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -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> 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 { std::vector 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(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(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(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(1.0 / std::sqrt(a.swh[i])); } } if (!std::isfinite(a.d) || a.d <= 0.0f) continue; diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index af484c75..69a3d661 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -221,14 +221,11 @@ SearchSpaceGroupResult SearchSpaceGroup( pass_absence[i] = in_range; } - // The merge floors the merged sigma at b|I| (Merge.h SigmaWithSystematicFloor) so that ISa = 1/b - // stays the asymptotic I/sigma ceiling: no reflection in the merge can read above it. A fixed - // I/sigma cut is therefore not a per-reflection test but a switch on the error model - every - // reflection at the floor reads 1/b exactly, however strong it is, and on a merge whose ISa falls - // below the cut NOTHING passes, so every operator is left with no pairs and the point group - // collapses to 1. Measured over the rotation battery: max I/sigma equals 1/b on every merge, and - // the four crystals that lose symmetry are exactly the four whose search merge sits at ISa 3.5-4.0 - // - just above this cut - and drops below 3 when the integration background changes. + // A fixed I/sigma cut is a statement about the error model as much as about a reflection: on a weak + // merge, whose error-model sigmas are large for every reflection, NOTHING passes, every operator is + // left with no pairs and the point group collapses to 1. Measured over the rotation battery, the four + // crystals that lose symmetry are exactly the four whose search merge is weakest - just above this + // cut - and that drop below it when the integration background changes. // So cap the cut at the merge's own I/sigma quantile: the correlation stage always keeps at least // its strongest quarter. This is a no-op on any merge where the fixed cut already keeps that many. constexpr double MIN_PRESENT_FRACTION = 0.25; @@ -754,14 +751,13 @@ SearchSpaceGroupResult SearchSpaceGroup( // screw-axis violation; it only relaxes "present", so it cannot over-call a screw whose // predicted-absent class carries real intensity. // - // present_cut, not the fixed cut: merged sigma is floored at b*|I| (SigmaWithSystematicFloor), - // so no reflection can read I/sigma above ISa = 1/b. On a merge whose ISa sits at or below the - // fixed cut NOTHING is ever present - screw_violations is identically zero, so every screw axis - // passes unchallenged, and present_strong is zero, so the centering rescue below switches - // itself off on exactly the weak data it exists for. Stage A already caps the cut at the - // merge's own 75th percentile; reusing it here keeps the two stages on one definition. Where - // the fixed cut is already the smaller of the two - any merge with ISa comfortably above it - - // present_cut EQUALS it and this is a no-op. + // present_cut, not the fixed cut: on a merge weak enough that nothing clears the fixed cut, + // screw_violations is identically zero, so every screw axis passes unchallenged, and + // present_strong is zero, so the centering rescue below switches itself off on exactly the + // weak data it exists for. Stage A already caps the cut at the merge's own 75th percentile; + // reusing it here keeps the two stages on one definition. Where the fixed cut is already the + // smaller of the two - any merge with a healthy I/sigma - present_cut EQUALS it and this is + // a no-op. const bool present = IoverSigma[i] > present_cut && (opt.present_e_squared <= 0.0 || Esq[i] > opt.present_e_squared); @@ -971,7 +967,7 @@ std::string SearchSpaceGroupResultToText(const SearchSpaceGroupResult& result, os << " absent/viol = reflections the group predicts absent, and how many are nonetheless present.\n" " E2 screw / E2 row = median I/(shell) of the SCREW-absent reflections and of the rest of\n" " their axial rows - a real screw leaves the first far below the second. The columns are\n" - " the centering evidence; they say little about screws (the merged sigma floor pins I/sigma).\n"; + " the centering evidence; they say little about screws.\n"; if (result.best_space_group.has_value()) { os << "Best space group: " << result.best_space_group->short_name(); diff --git a/image_analysis/scale_merge/SearchSpaceGroup.h b/image_analysis/scale_merge/SearchSpaceGroup.h index 9326ddfb..765dcd73 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.h +++ b/image_analysis/scale_merge/SearchSpaceGroup.h @@ -47,8 +47,8 @@ struct SpaceGroupCandidateScore { double present_mean_i_over_sigma = 0.0; // Screw evidence, as median E^2 = I/(shell): the reflections a screw axis extinguishes against // the rest of that same axial row. This is the comparison the screw test makes, and the only one - // that means anything for it - the I/sigma means above are dominated by the merged-sigma floor - // (sigma >= b|I| pins I/sigma at ISa), so they read the same for absent and present alike. + // that means anything for it - the I/sigma means above are dominated by the error model's + // intensity-proportional term, so they read much the same for absent and present alike. double screw_absent_median_e_squared = 0.0; double screw_row_median_e_squared = 0.0; bool consistent = false; // absent class confirmed weak (few violations) diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index e7a39f82..8564a81c 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1382,7 +1382,7 @@ static int RunRugnux(int argc, char **argv) { merge_engine.FilterByImageCC(experiment.GetScalingSettings().GetMinCCForImage() > 0.0); // Fit the (a, b) error model from symmetry-mate scatter before merging, exactly as the full // pipeline does (Rugnux.cpp). Without this the offline --scale merge would use the identity - // model and produce much worse stills intensities (no (b*I)^2 systematic term, no sigma floor). + // model and produce much worse stills intensities (no (b*I)^2 systematic term at all). merge_engine.RefineErrorModel(reflections); if (merge_engine.ErrorModelActive()) logger.Info("Error model: a={:.3f} b={:.3f} ISa={:.1f} chi2={:.2f}", merge_engine.ErrorModelA(), diff --git a/tests/SearchSpaceGroupTwinTest.cpp b/tests/SearchSpaceGroupTwinTest.cpp index e9a46e88..46bc8814 100644 --- a/tests/SearchSpaceGroupTwinTest.cpp +++ b/tests/SearchSpaceGroupTwinTest.cpp @@ -76,8 +76,8 @@ namespace { // model's b matching the systematic scatter that is actually there. {"sigmas 1.7x too small", 1.7, 0.05, std::nullopt}, // The other way a fitted error model misses: the statistical sigmas come out somewhat too - // LARGE while b - the asymptotic I/sigma ceiling, ISa = 1/b - is fitted 3x too optimistic, so - // the systematic scatter present is 3x what the merged sigmas' floor admits. + // LARGE while b - the asymptotic per-observation I/sigma, ISa = 1/b - is fitted 3x too + // optimistic, so the systematic scatter present is 3x what the merged sigmas admit. {"ISa 3x too optimistic", 0.6, 0.02, 0.06}, }; @@ -259,8 +259,8 @@ TEST_CASE("SearchSpaceGroup on a perfect merohedral twin returns one of the two } } -// THE property this harness exists for. Multiplicity changes only the sigmas - the random part of a -// merged sigma averages down as 1/sqrt(n) while the systematic floor b*|I| does not - so it changes +// THE property this harness exists for. Multiplicity changes only the sigmas - a merged sigma averages +// down as 1/sqrt(n) while the systematic error the crystal carries does not - so it changes // how well the SAME crystal is measured, never what its symmetry is. A symmetry decision that moves // when the same crystal is merged 2x instead of 18x is a defect of the criterion, not a property of // the data. diff --git a/tests/SyntheticMergedReflections.h b/tests/SyntheticMergedReflections.h index cc34e58d..395c487a 100644 --- a/tests/SyntheticMergedReflections.h +++ b/tests/SyntheticMergedReflections.h @@ -32,12 +32,11 @@ // alpha = 0.5 the two are identical and the twin is indistinguishable from real symmetry. // Setting the true group to the SUPERgroup instead gives the untwinned high-symmetry control. // -// * the MERGE MULTIPLICITY, modelled the way the real merge behaves (Merge.h, -// SigmaWithSystematicFloor): the random part of the merged sigma averages down as -// 1/sqrt(multiplicity) while the systematic part (b*I - absorption, partiality, beam flicker; -// correlated across a reflection's repeats) does not, so the merged sigma is -// max(sigma_statistical, b*|I|). Multiplicity therefore changes the sigmas but NOT the physics, -// and no symmetry decision may depend on it. +// * the MERGE MULTIPLICITY, modelled the way the real merge behaves: the error model gives one +// observation sigma^2 = sigma_counting^2 + (b*I)^2 and the inverse-variance merge of n of them +// divides that by n, so the merged sigma is sqrt(sigma_counting^2 + (b*I)^2)/sqrt(multiplicity). +// Multiplicity therefore changes the sigmas but NOT the physics, and no symmetry decision may +// depend on it. // // * an ERROR-MODEL MISCALIBRATION - real merged sigmas come out under-estimated (~1.7x), which is // what pushes the merge's reduced chi^2 to ~3 and switches SearchSpaceGroup between its @@ -201,7 +200,7 @@ namespace jfjoch_test { // Statistical error of one observation, and of the merge of n of them. const double sigma_one = std::sqrt(i_obs + p.background_variance); // Systematic error: a property of the reflection, identical in every observation - // of it, so it survives the merge - this is what the b*|I| sigma floor models. + // of it, so it survives the merge (and the merged sigma does not know about it). const double systematic = p.true_systematic_b.value_or(p.error_model_b) * i_obs * gauss(rng); @@ -218,18 +217,18 @@ namespace jfjoch_test { const double i_half = i_obs + systematic + sigma_stat_half * gauss(rng); r.I_half[half] = static_cast(i_half); r.sigma_half[half] = static_cast( - std::max(sigma_stat_half, p.error_model_b * std::abs(i_half)) / - p.sigma_miscalibration); + std::hypot(sigma_one, p.error_model_b * i_half) / + std::sqrt(static_cast(n_half[half])) / p.sigma_miscalibration); sum_n_i += n_half[half] * i_half; } const double i_merged = sum_n_i / n_obs; - const double sigma_stat = sigma_one / std::sqrt(static_cast(n_obs)); r.I = static_cast(i_merged); - // Merge.h SigmaWithSystematicFloor, then thrown off by the error-model - // miscalibration. - r.sigma = static_cast( - std::max(sigma_stat, p.error_model_b * std::abs(i_merged)) / p.sigma_miscalibration); + // The error model on one observation, averaged down by the merge, then thrown off + // by the error-model miscalibration. + r.sigma = static_cast(std::hypot(sigma_one, p.error_model_b * i_merged) / + std::sqrt(static_cast(n_obs)) / + p.sigma_miscalibration); merged.push_back(r); } From 97dbbc50b4cf844cb77fa19ca8cb3ba710d761b2 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 9 Aug 2026 20:58:35 +0200 Subject: [PATCH 2/2] Merging: fit the error model on the reflections the cutoff keeps The (a, b) fit ran over the whole merged range and the automatic resolution cutoff was applied afterwards, so the sigma correction applied to the reflections that survive was calibrated largely on reflections that do not. Measured on one dataset: a = 0.286 fitted over 843k reflections, 22k written. A manual --scaling-high-resolution already restricts the population at ingest, so only the automatic path was affected. Fit over the full range, merge, read the cutoff from that merge, refit (a, b) on the samples the cutoff keeps, merge again. The circularity resolves by direction: the cutoff comes from CC1/2, a correlation of the two half-set means, which the sigma scale barely moves, so the cutoff can be read first and the sigmas calibrated on the population it chose. One refinement, not an iteration; one extra merge pass. Note this is invisible to the rotation battery, which passes an explicit high resolution limit matched to XDS and so never exercises the automatic cutoff. With a manual limit the fitted (a, b) are byte-identical to before. The equivalent defect in the stills / offline --scale path is untouched; it is a different engine and needs its own validation. Co-Authored-By: Claude Opus 5 (1M context) --- .../scale_merge/RotationScaleMerge.cpp | 350 ++++++++++-------- 1 file changed, 194 insertions(+), 156 deletions(-) diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 9e5d055d..f77c3a7a 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -1675,67 +1675,16 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool std::vector reject_median(n_groups, NAN); double error_model_a = 1.0, error_model_b = 0.0, error_model_chi2 = 0.0; bool error_model_active = false; - { - struct Sample { double s2, I2, dev2; }; - std::vector samples; - std::vector cnt(n_groups, 0); // per-group usable count (both paths; feeds reject-median) - bool did_gpu = false; -#ifdef JFJOCH_USE_CUDA - if (use_gpu_merge) { - const int nf = static_cast(fulls.size()); - std::vector gs2(nf), gI2(nf), gdev2(nf); - std::vector gvalid(nf); - gpu_->MergeEmSamples(for_search, min_partiality, em_mean.data(), cnt.data(), - gs2.data(), gI2.data(), gdev2.data(), gvalid.data()); - samples.reserve(nf); - for (int i = 0; i < nf; ++i) - if (gvalid[i]) samples.push_back({gs2[i], gI2[i], gdev2[i]}); - did_gpu = true; - } -#endif - if (!did_gpu) { - // Per-group inverse-variance mean over usable fulls (>=2 obs), and the leverage-corrected samples. - std::vector sw(n_groups, 0.0), swI(n_groups, 0.0); - for (const auto &o : fulls) { - if (!usable_merge(o)) continue; - const double sigma_corr = static_cast(o.sigma) * o.corr; - const double w = 1.0 / (sigma_corr * sigma_corr); - sw[o.group] += w; swI[o.group] += w * (static_cast(o.I) * o.corr); cnt[o.group]++; - } - for (int g = 0; g < n_groups; ++g) - if (cnt[g] >= 2 && sw[g] > 0.0) em_mean[g] = swI[g] / sw[g]; - - samples.reserve(fulls.size()); - for (const auto &o : fulls) { - if (!usable_merge(o) || cnt[o.group] < 2) continue; - const double mean = em_mean[o.group]; - if (!std::isfinite(mean)) continue; - const double sigma_corr = static_cast(o.sigma) * o.corr; - const double s2 = sigma_corr * sigma_corr; - const double w = 1.0 / s2; - const double factor = 1.0 - w / sw[o.group]; - if (factor < 0.05) continue; - const double resid = static_cast(o.I) * o.corr - mean; - samples.push_back({s2, mean * mean, resid * resid / factor}); - } - } - - // Per-group outlier-rejection median of I*corr (host both paths - a per-group median is awkward on - // the GPU; cheap here, cnt >= 2 filter from the em pass). Fed to the merge accumulate. - if (reject_outliers) { - std::vector> iv(n_groups); - for (const auto &o : fulls) - if (usable_merge(o) && cnt[o.group] >= 2) - iv[o.group].push_back(o.I * o.corr); - for (int g = 0; g < n_groups; ++g) - if (!iv[g].empty()) { - std::nth_element(iv[g].begin(), iv[g].begin() + iv[g].size() / 2, iv[g].end()); - reject_median[g] = iv[g][iv[g].size() / 2]; - } - } - constexpr int n_bins = 16; - // Fit (a, b) from the intensity-binned median deviations. Factored into a lambda so it can be - // re-run on a misfit-free pool below; takes the samples by value (it sorts them in place). + // One leverage-corrected sample per usable full: its raw variance, its group's mean intensity, its + // squared deviation from that mean - and the resolution it sits at, because the fit is re-run below + // over the samples that survive the automatic resolution cutoff. + struct Sample { double s2, I2, dev2; float d; }; + std::vector samples; + constexpr int n_bins = 16; + // Fit (a, b) from the intensity-binned median deviations of a pool of samples, then refit on that + // pool's misfit-free subset. A lambda because the pool changes once the cutoff below is known. + auto fit_error_model = [&](const std::vector &pool) { + // Takes the samples by value (it sorts them in place). auto fit_ab = [&](std::vector smp) { if (smp.size() < static_cast(8 * n_bins)) return; @@ -1778,22 +1727,80 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool error_model_chi2 = chi2.empty() ? 0.0 : median_of(chi2) / CHI2_1_MEDIAN; } }; - fit_ab(samples); + fit_ab(pool); // Refit on a misfit-free pool: the merge drops symmetry outliers (|I - median| > reject_nsigma * // sigma) from the merged intensity, so drop the equivalent samples (dev2 > reject_nsigma^2 * model // variance) from the error-model fit too, keeping the fitted sigmas consistent with the reflections - // that actually survive. Operates on the shared `samples`, so CPU and GPU stay bit-identical. + // that actually survive. Operates on the shared samples, so CPU and GPU stay bit-identical. if (reject_outliers && error_model_active) { const double ns2 = reject_nsigma * reject_nsigma, b2 = error_model_b * error_model_b; std::vector kept; - kept.reserve(samples.size()); - for (const auto &s : samples) { + kept.reserve(pool.size()); + for (const auto &s : pool) { const double v = error_model_a * s.s2 + b2 * s.I2; if (v > 0.0 && s.dev2 <= ns2 * v) kept.push_back(s); } - if (kept.size() >= static_cast(8 * n_bins) && kept.size() < samples.size()) + if (kept.size() >= static_cast(8 * n_bins) && kept.size() < pool.size()) fit_ab(std::move(kept)); } + }; + { + std::vector cnt(n_groups, 0); // per-group usable count (both paths; feeds reject-median) + bool did_gpu = false; +#ifdef JFJOCH_USE_CUDA + if (use_gpu_merge) { + const int nf = static_cast(fulls.size()); + std::vector gs2(nf), gI2(nf), gdev2(nf); + std::vector gvalid(nf); + gpu_->MergeEmSamples(for_search, min_partiality, em_mean.data(), cnt.data(), + gs2.data(), gI2.data(), gdev2.data(), gvalid.data()); + samples.reserve(nf); + for (int i = 0; i < nf; ++i) + if (gvalid[i]) samples.push_back({gs2[i], gI2[i], gdev2[i], fulls[i].d}); + did_gpu = true; + } +#endif + if (!did_gpu) { + // Per-group inverse-variance mean over usable fulls (>=2 obs), and the leverage-corrected samples. + std::vector sw(n_groups, 0.0), swI(n_groups, 0.0); + for (const auto &o : fulls) { + if (!usable_merge(o)) continue; + const double sigma_corr = static_cast(o.sigma) * o.corr; + const double w = 1.0 / (sigma_corr * sigma_corr); + sw[o.group] += w; swI[o.group] += w * (static_cast(o.I) * o.corr); cnt[o.group]++; + } + for (int g = 0; g < n_groups; ++g) + if (cnt[g] >= 2 && sw[g] > 0.0) em_mean[g] = swI[g] / sw[g]; + + samples.reserve(fulls.size()); + for (const auto &o : fulls) { + if (!usable_merge(o) || cnt[o.group] < 2) continue; + const double mean = em_mean[o.group]; + if (!std::isfinite(mean)) continue; + const double sigma_corr = static_cast(o.sigma) * o.corr; + const double s2 = sigma_corr * sigma_corr; + const double w = 1.0 / s2; + const double factor = 1.0 - w / sw[o.group]; + if (factor < 0.05) continue; + const double resid = static_cast(o.I) * o.corr - mean; + samples.push_back({s2, mean * mean, resid * resid / factor, o.d}); + } + } + + // Per-group outlier-rejection median of I*corr (host both paths - a per-group median is awkward on + // the GPU; cheap here, cnt >= 2 filter from the em pass). Fed to the merge accumulate. + if (reject_outliers) { + std::vector> iv(n_groups); + for (const auto &o : fulls) + if (usable_merge(o) && cnt[o.group] >= 2) + iv[o.group].push_back(o.I * o.corr); + for (int g = 0; g < n_groups; ++g) + if (!iv[g].empty()) { + std::nth_element(iv[g].begin(), iv[g].begin() + iv[g].size() / 2, iv[g].end()); + reject_median[g] = iv[g][iv[g].size() / 2]; + } + } + fit_error_model(samples); } // Asymptotic I/sigma. ISa is by definition the I -> infinity limit of the signal-to-noise, i.e. the // reproducibility of the strongest reflections (Diederichs, Acta Cryst. D66 (2010), 733-740). The @@ -1803,8 +1810,10 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // fractional scatter of its symmetry mates estimates the systematic term, and the robust median over // 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) { + double error_model_b_asymptotic = 0.0; + auto estimate_asymptote = [&]() { + error_model_b_asymptotic = error_model_b; + if (!error_model_active) return; struct GroupScatter { double sum = 0, sum_sq = 0, sum_var = 0; int n = 0; }; std::vector gs(n_groups); for (const auto &o : fulls) { @@ -1863,17 +1872,8 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // 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. - 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; - if (error_model_active) - logger.Info("Error model: a={:.3f} b={:.3f} ISa={:.1f} chi2={:.2f}", error_model_a, error_model_b, - isa_reported, error_model_chi2); + }; + estimate_asymptote(); auto corrected_sigma = [&](float I_corr, float sigma_corr, int g) -> float { if (!error_model_active) return sigma_corr; @@ -1882,94 +1882,132 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool + (error_model_b * I_for_b) * (error_model_b * I_for_b); return v > 0.0 ? static_cast(std::sqrt(v)) : sigma_corr; }; - // ---- Merge: per-group inverse-variance sums with corrected sigma + deterministic half sets. ---- + // ---- Merge: per-group inverse-variance sums with corrected sigma + deterministic half sets, then + // export. In a lambda because the automatic resolution cutoff below re-runs it once the error + // model has been refitted on the reflections that survive the cut. ---- struct Accum { double swI = 0, sw = 0, swIh[2] = {0, 0}, swh[2] = {0, 0}; size_t nh[2] = {0, 0}; float d = NAN; }; - std::vector acc(n_groups); - size_t reject_count = 0; - std::vector rejected_obs(fulls.size(), 0); // per-full outlier-rejected flag (both paths) - bool did_gpu_acc = false; -#ifdef JFJOCH_USE_CUDA - if (use_gpu_merge) { - std::vector aswI(n_groups), asw(n_groups), aswIh0(n_groups), aswIh1(n_groups), - aswh0(n_groups), aswh1(n_groups), ad(n_groups); - std::vector anh0(n_groups), anh1(n_groups), arej(n_groups); - gpu_->MergeAccum(error_model_a, error_model_b, error_model_active, - reject_outliers, reject_nsigma, reject_median.data(), - aswI.data(), asw.data(), aswIh0.data(), aswIh1.data(), - aswh0.data(), aswh1.data(), anh0.data(), anh1.data(), ad.data(), arej.data(), - rejected_obs.data()); - for (int g = 0; g < n_groups; ++g) { - Accum &a = acc[g]; - a.swI = aswI[g]; a.sw = asw[g]; a.swIh[0] = aswIh0[g]; a.swIh[1] = aswIh1[g]; - a.swh[0] = aswh0[g]; a.swh[1] = aswh1[g]; - a.nh[0] = static_cast(anh0[g]); a.nh[1] = static_cast(anh1[g]); - a.d = static_cast(ad[g]); - reject_count += static_cast(arej[g]); - } - did_gpu_acc = true; - } -#endif - if (!did_gpu_acc) - for (const auto &o : fulls) { - if (!usable_merge(o)) continue; - const float I_corr = o.I * o.corr; - float sigma_corr = o.sigma * o.corr; - sigma_corr = corrected_sigma(I_corr, sigma_corr, o.group); - if (reject_outliers && error_model_active && std::isfinite(reject_median[o.group]) - && std::fabs(I_corr - reject_median[o.group]) > reject_nsigma * sigma_corr) { - ++reject_count; - rejected_obs[&o - fulls.data()] = 1; - continue; - } - const double w = 1.0 / (static_cast(sigma_corr) * sigma_corr); - const double wI = w * I_corr; - const int half = HalfForImage(o.frame); - auto &a = acc[o.group]; - a.swI += wI; a.sw += w; - a.swIh[half] += wI; a.swh[half] += w; a.nh[half]++; - if (!std::isfinite(a.d) && std::isfinite(o.d) && o.d > 0.0f) a.d = o.d; - } - // ---- Export merged reflections (+ resolution-shell R-free flags). ---- + std::vector acc; Result result; - result.isa = isa_reported; - std::vector merged_I(n_groups, NAN); - float d_min = std::numeric_limits::max(), d_max = 0.0f; - for (int g = 0; g < n_groups; ++g) { - const auto &a = acc[g]; - if (a.sw <= 0.0) continue; - MergedReflection mr{}; - mr.h = group_h[g]; mr.k = group_k[g]; mr.l = group_l[g]; - mr.I = static_cast(a.swI / a.sw); - // 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(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(a.swIh[i] / a.swh[i]); - mr.sigma_half[i] = static_cast(1.0 / std::sqrt(a.swh[i])); + std::vector merged_I; + size_t reject_count = 0; + std::vector rejected_obs; // per-full outlier-rejected flag (both paths) + auto run_merge = [&]() { + acc.assign(n_groups, Accum{}); + result.merged.clear(); + merged_I.assign(n_groups, NAN); + reject_count = 0; + rejected_obs.assign(fulls.size(), 0); + bool did_gpu_acc = false; +#ifdef JFJOCH_USE_CUDA + if (use_gpu_merge) { + std::vector aswI(n_groups), asw(n_groups), aswIh0(n_groups), aswIh1(n_groups), + aswh0(n_groups), aswh1(n_groups), ad(n_groups); + std::vector anh0(n_groups), anh1(n_groups), arej(n_groups); + gpu_->MergeAccum(error_model_a, error_model_b, error_model_active, + reject_outliers, reject_nsigma, reject_median.data(), + aswI.data(), asw.data(), aswIh0.data(), aswIh1.data(), + aswh0.data(), aswh1.data(), anh0.data(), anh1.data(), ad.data(), arej.data(), + rejected_obs.data()); + for (int g = 0; g < n_groups; ++g) { + Accum &a = acc[g]; + a.swI = aswI[g]; a.sw = asw[g]; a.swIh[0] = aswIh0[g]; a.swIh[1] = aswIh1[g]; + a.swh[0] = aswh0[g]; a.swh[1] = aswh1[g]; + a.nh[0] = static_cast(anh0[g]); a.nh[1] = static_cast(anh1[g]); + a.d = static_cast(ad[g]); + reject_count += static_cast(arej[g]); } + did_gpu_acc = true; } - if (!std::isfinite(a.d) || a.d <= 0.0f) continue; - d_min = std::min(d_min, a.d); - d_max = std::max(d_max, a.d); - merged_I[g] = mr.I; - result.merged.push_back(mr); - } +#endif + if (!did_gpu_acc) + for (const auto &o : fulls) { + if (!usable_merge(o)) continue; + const float I_corr = o.I * o.corr; + float sigma_corr = o.sigma * o.corr; + sigma_corr = corrected_sigma(I_corr, sigma_corr, o.group); + if (reject_outliers && error_model_active && std::isfinite(reject_median[o.group]) + && std::fabs(I_corr - reject_median[o.group]) > reject_nsigma * sigma_corr) { + ++reject_count; + rejected_obs[&o - fulls.data()] = 1; + continue; + } + const double w = 1.0 / (static_cast(sigma_corr) * sigma_corr); + const double wI = w * I_corr; + const int half = HalfForImage(o.frame); + auto &a = acc[o.group]; + a.swI += wI; a.sw += w; + a.swIh[half] += wI; a.swh[half] += w; a.nh[half]++; + if (!std::isfinite(a.d) && std::isfinite(o.d) && o.d > 0.0f) a.d = o.d; + } + // ---- Export merged reflections. ---- + for (int g = 0; g < n_groups; ++g) { + const auto &a = acc[g]; + if (a.sw <= 0.0) continue; + MergedReflection mr{}; + mr.h = group_h[g]; mr.k = group_k[g]; mr.l = group_l[g]; + mr.I = static_cast(a.swI / a.sw); + // 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(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(a.swIh[i] / a.swh[i]); + mr.sigma_half[i] = static_cast(1.0 / std::sqrt(a.swh[i])); + } + } + if (!std::isfinite(a.d) || a.d <= 0.0f) continue; + merged_I[g] = mr.I; + result.merged.push_back(mr); + } + }; + run_merge(); // Automatic high-resolution cutoff (post-merge): trim the written reflections + reported shells to - // the CC1/2 fall-off. The scaling, combine and error model above already ran over the full range, - // and the per-image _process.h5 is written elsewhere from the partials, so no data is lost. A manual + // the CC1/2 fall-off. The scaling and combine above ran over the full range, and the per-image + // _process.h5 is written elsewhere from the partials, so no data is lost. A manual // --scaling-high-resolution (d_min_limit) wins; the P1 search merge (for_search) is never cut, so // the space-group search still sees the full range. const std::optional effective_d_min = ApplyResolutionCutoff( result.merged, d_min_limit, resolution_cutoff_method, resolution_cc_target, for_search, logger); + // The error model has to be calibrated on the reflections that are kept, not on the ones that are + // thrown away: on a default run the cut can remove the majority of the measured range, and a fit + // spanning it is dominated by reflections that are not written (measured on one run: a = 0.28, + // b = 0.159 over the full range against a = 0.42, b = 0.119 over the kept one, on 40% of the + // observations the whole-range fit saw). A manual limit needs nothing here - it already + // restricted the observations at ingest - so this is the automatic cut catching up with it. + // The circularity is resolved by direction: the cutoff is read from the provisional merge, and + // CC1/2 is a correlation of the two half-set MEANS, which the sigma scale barely moves; the sigmas + // are then calibrated on the population the cutoff chose. One refinement, not an iteration. + if (effective_d_min && effective_d_min != d_min_limit) { + std::vector in_range; + in_range.reserve(samples.size()); + for (const auto &s : samples) + if (s.d >= *effective_d_min) in_range.push_back(s); + fit_error_model(in_range); + estimate_asymptote(); + run_merge(); + std::erase_if(result.merged, [&](const MergedReflection &m) { + return std::isfinite(m.d) && m.d < *effective_d_min; + }); + } + + // 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. + constexpr double MIN_ASYMPTOTIC_B = 0.01; // ISa cap 100 + result.isa = error_model_b_asymptotic >= MIN_ASYMPTOTIC_B ? 1.0 / error_model_b_asymptotic : 0.0; + if (error_model_active) + logger.Info("Error model: a={:.3f} b={:.3f} ISa={:.1f} chi2={:.2f}", error_model_a, error_model_b, + result.isa, error_model_chi2); + AssignRfreeFlags(result.merged, x.GetSpaceGroupNumber().value_or(1), rfree_fraction); // French-Wilson (F, and F(+)/F(-) from the anomalous split) is deferred until after the anomalous // accumulator below has attached I(+)/I(-), so the two hands get their amplitudes in one pass.