From 7bbf072ad2e8366dbc1643615c23aa16e5e6877e Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 11 Aug 2026 06:14:46 +0200 Subject: [PATCH] Scaling: do not report an ISa that was never measured The error model's systematic term b is identified only by the spread of I^2/sigma^2 across the intensity bins the fit uses, and those bins hold equal COUNTS. So when fewer reflections are strong than one bin holds - a sixteenth of the pool - the top bin's median sits at an intensity where b cannot be measured at all, and the fit hands it the bins' own noise-selection slope instead: sorting noise by its group mean squared makes dev2 rise with I2 even when the true b is zero, and with no strong bin to out-vote it that slope becomes b. The result is not a small error. On the battery's weakest crystal, 2.2% of whose fulls reach I/sigma 2, the fit returns b = 5.6 - sigma -> 2*I at the strong end - and since corrected_sigma applies b at the GROUP MEAN, sigma^2 = a*sigma^2 + (b*mean)^2 is a per-group constant that caps merged |I/sigma| at sqrt(n)/b. The cap lands at 2.3, so 98.8% of merged reflections come out below 3 and the reported ISa is 0.50, on data whose CC1/2 is 99.3% at multiplicity 18.7. XDS fits 6.13 from the same images. Feeding XDS's own scaled observations through this estimator returns 0.84, so it is the estimator and not the data; synthetic data built with b = 0 and 1.8% strong reproduces a = 0.51 and ISa 0.50 to two digits, and recovers the truth as soon as the strong fraction passes one bin. So refuse to report what was not measured: when the strongest bin's own (I/sigma)^2 is below 4, fit a alone, hold b at zero and warn that ISa is unmeasured. The threshold is not delicate - the two crystals it fires on sit at 0.22 and 0.84 while the next crystal in the battery is at 31.7 and a healthy one at 342, so anything from 4 to 25 selects the same two. Full 38-crystal rotation battery: it fires on those two crystals and no others, and space groups are unchanged at 35/38. Dropping the spurious term also fixes the merge weights it had been distorting - on the worse of the two, R_meas 19.2 -> 13.6%, low-resolution R_meas 13.8 -> 6.9% against XDS's 14.1%, CC1/2 98.7 -> 100.0%, with chi2 1.11 on the one-parameter model. Two further crystals move slightly; the guard never fires on either, and they are marginal crystals of the kind whose two-pass branch any recompilation can shift. This reports the parameter as unmeasured rather than clamping it to something plausible, because the honest statement is that the data do not reach far enough for a systematic error to be seen - not that there is none. Co-Authored-By: Claude Opus 5 (1M context) --- docs/CHANGELOG.md | 1 + .../scale_merge/RotationScaleMerge.cpp | 49 +++++++++++++++++-- 2 files changed, 46 insertions(+), 4 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 04796300..b4637fad 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -5,6 +5,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as * rugnux: The error-model **a** and **b** are reported in XDS's convention, and `_reflns.jfjoch_diffrn_ISa` now carries the whole-range `1/sqrt(a*b)` that XDS's ISa denotes; the strong-reflection asymptote moves to `_reflns.jfjoch_diffrn_ISa_asymptotic`. **A file written by an earlier version carries the asymptote under the old name.** * Bragg integration: the background ring's outer radius default changes from 10 px to **13 px**, which roughly doubles the pixels behind each background estimate; the signal disk is unchanged. +* Scaling: when too few reflections are strong enough to constrain the error model's systematic term **b**, it is now held at zero and **ISa is reported as unmeasured** rather than as the spurious value the fit would otherwise extrapolate. * Scaling: the partials of one rocking event now share a single exact-Bragg angle, instead of each solving its own from its frame's refined orientation, which had been feeding per-frame refinement noise into partiality in proportion to 1/zeta. * Bragg integration: signal pixels shared with a neighbouring reflection are now dropped from the profile fit (`--overlap off|reject|exclude`, **default `exclude`**), so a crowded reflection no longer reads its neighbour's flux as its own; `reject` instead discards a reflection whose cleanly-observed profile fraction is below `--overlap-minpk` (default 0.75). * rugnux: The background ring can be **elongated radially per reflection** for broadband data (`--integration-stencil `, default 0 = the fixed circular ring), by `k` times the beam's radial streak; the `r1` signal box stays circular. diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index e494c8af..f9e43a73 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -1795,6 +1795,7 @@ 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; + bool error_model_b_unmeasured = false; // b had no leverage; ISa is not reported // 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. @@ -1822,6 +1823,22 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool bI2.push_back(median_of(vI2)); bd2.push_back(median_of(vd2) / CHI2_1_MEDIAN); } + // `b` is identified ONLY by the spread of I^2/sigma^2 across the bins, and the bins hold equal + // COUNTS - so when fewer reflections are strong than one bin holds (1/16 of the pool), the top + // bin's median sits at an intensity where b cannot be measured, and the fit assigns it the + // bins' own noise-selection slope instead. Measured: a crystal with 2.2% of its fulls above + // I/sigma 2 reports b = 5.6, i.e. sigma -> 2*I at the strong end and ISa 0.50, while XDS gets + // 6.13 from the same images; synthetic data with b_true = 0 and 1.8% strong reproduces that to + // two digits, and recovers the truth as soon as the strong fraction passes 1/16. sigma^2 is + // then a*sigma^2 + (b*mean)^2 with a per-group constant, which caps merged |I/sigma| at + // sqrt(n)/b - reporting an impossible systematic error on data whose CC1/2 is 99%. + // + // So refuse to report what was not measured: if the strongest bin's own (I/sigma)^2 is small, + // fit `a` alone and leave b at zero. The threshold is not delicate - over the rotation battery + // the two crystals this fires on sit at 0.22 and 0.84 while the next is 31.7 and a healthy one + // is 342, so anything from 4 to 25 selects the same two and nothing else changes. + constexpr double B_LEVER_MIN = 4.0; + const bool b_measurable = bI2.back() > B_LEVER_MIN * bs2.back(); std::vector bd2_sorted = bd2; const double dev2_floor = std::max(1e-30, 1e-3 * median_of(bd2_sorted)); double Ass = 0, AsI = 0, AII = 0, Bs = 0, BI = 0; @@ -1833,10 +1850,26 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool Bs += wgt * s2 * d2; BI += wgt * I2 * d2; } const double det = Ass * AII - AsI * AsI; - if (det > 1e-10 * Ass * AII) { + if (!b_measurable) { + // One parameter, same weights: dev2 = a*sigma^2 alone. + if (Ass > 0.0) { + error_model_a = std::clamp(Bs / Ass, 0.25, 100.0); + error_model_b = 0.0; + error_model_b_unmeasured = true; + error_model_active = true; + std::vector chi2; + chi2.reserve(smp.size()); + for (const auto &s : smp) { + const double v = error_model_a * s.s2; + if (v > 0.0) chi2.push_back(s.dev2 / v); + } + error_model_chi2 = chi2.empty() ? 0.0 : median_of(chi2) / CHI2_1_MEDIAN; + } + } else if (det > 1e-10 * Ass * AII) { error_model_a = std::clamp((Bs * AII - BI * AsI) / det, 0.25, 100.0); const double b2 = std::max((Ass * BI - AsI * Bs) / det, 0.0); error_model_b = std::sqrt(b2); + error_model_b_unmeasured = false; error_model_active = true; std::vector chi2; chi2.reserve(smp.size()); @@ -2136,9 +2169,17 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // is what XDS's ISa means, and the strong-reflection asymptote above, which XDS has no // equivalent of and which can only ever be the more optimistic of the two. const auto em = ToXdsErrorModel(error_model_a, error_model_b); - logger.Info("Error model (XDS convention): a={:.3f} b={:.3e} ISa={:.1f} " - "(strong-reflection asymptote {:.1f}) chi2={:.2f}", - em.a, em.b, em.isa, result.isa_asymptotic, error_model_chi2); + if (error_model_b_unmeasured) + logger.Warning("Error model (XDS convention): a={:.3f}, b NOT MEASURABLE - fewer than one " + "intensity bin's worth of reflections are strong enough to constrain it, so " + "b is held at 0 and ISa is not reported. chi2={:.2f}. This says the data do " + "not reach far enough for a systematic error to be seen, not that there is " + "none; a resolution range matched to the signal would measure it", + em.a, error_model_chi2); + else + logger.Info("Error model (XDS convention): a={:.3f} b={:.3e} ISa={:.1f} " + "(strong-reflection asymptote {:.1f}) chi2={:.2f}", + em.a, em.b, em.isa, result.isa_asymptotic, error_model_chi2); } AssignRfreeFlags(result.merged, x.GetSpaceGroupNumber().value_or(1), rfree_fraction);