From 06b8c8ed663c4e7e87c3e3b194be95933f65d193 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 6 Aug 2026 19:06:00 +0200 Subject: [PATCH] Merge statistics: count the observations the merge kept, not the ones it walked Whenever the merge-time ice-ring mask dropped a band, the per-shell observation count and hence the reported multiplicity were wrong. On one crystal the lowest resolution shell read 40780 observations over 1932 unique reflections - 21.1x - where the truth is 27007 and 13.98x, and the overall redundancy read 12.52 against 12.29. Only counts were affected: intensities, sigmas, R_meas, CC1/2, completeness and ISa were right throughout, because a masked group carries merged_I = NaN and never enters those sums. It looked like double counting and was not - it is a MOVE. Two independent faults, both in three lines: total_obs rides on the R_meas re-walk, whose filter deliberately ignores the ring mask (and, on a search pass, the ice flag) so that R_meas is computed on the same reflections either way. RmeasUsable therefore differs from MergeUsable by exactly those two tests, and the observations they admit were being counted against a `unique` that excludes them. On the GPU path that count is binned by the GROUP's resolution, and a group every one of whose observations is masked never has one written - acc[g].d stays NaN. ResolutionShells::GetShell(NaN) then returned shell 0 rather than nothing: NaN fails both bound comparisons, falls through to the arithmetic, and static_cast(NaN) is INT_MIN, which the clamp maps to 0. So the masked ring's observations were re-labelled into the lowest-resolution shell, four shells from the ring they came from. The two paths disagreeing on the same run is what settled it: with the mask on, the GPU statistics gave shell 0 = 752 and the CPU statistics 423, while the merged intensities were identical. Count the merged population instead - acc[g].nh, which the merge already accumulates per group - and guard the CPU increment with usable_merge. The rnusable skip stays: any group present in the merged output has at least one observation passing MergeUsable, and MergeUsable is a subset of RmeasUsable, so it cannot drop a group that contributes to `unique`. With the mask off and for_search false the two predicates are identical, so this is provably inert on every shipped configuration - demonstrated on four configurations, including one where ice handling is active but the mask does not fire: the statistics blocks are unchanged. (The reflection lists differ in the last ulp on 3-12% of lines, but so do two runs of the same binary; that is the known rotation nondeterminism, and the statistics block is what is stable.) The NaN guard also removes a silent contamination nobody was looking for. Four call sites validate a resolution with `d <= 0`, which NaN passes: the Wilson-B fit and per-shell (CalcISigma), the per-image resolution plot (SpotUtils) and the shell Wilson prior (FrenchWilson) were all binning non-finite d into their lowest-resolution shell. French-Wilson now falls back to the global mean rather than to that shell's, which is the worst prior available. Co-Authored-By: Claude Opus 5 (1M context) --- common/ResolutionShells.cpp | 5 ++++- image_analysis/scale_merge/Merge.cpp | 6 +++++- image_analysis/scale_merge/RotationScaleMerge.cpp | 13 +++++++++++-- image_analysis/scale_merge/RotationScaleMergeGPU.cu | 6 ++++-- image_analysis/scale_merge/RotationScaleMergeGPU.h | 3 ++- 5 files changed, 26 insertions(+), 7 deletions(-) diff --git a/common/ResolutionShells.cpp b/common/ResolutionShells.cpp index 828a64c7..cac11f04 100644 --- a/common/ResolutionShells.cpp +++ b/common/ResolutionShells.cpp @@ -23,7 +23,10 @@ ResolutionShells::ResolutionShells(float d_min, float d_max, int32_t nshells) } std::optional ResolutionShells::GetShell(float d) const { - if (d <= d_min || d > d_max) + // NaN fails every comparison, so without the explicit test it would fall through to the + // arithmetic below, where static_cast(NaN) is INT_MIN and the clamp turns it into + // shell 0 - silently binning "no resolution" as the lowest-resolution shell. + if (!std::isfinite(d) || d <= d_min || d > d_max) return {}; if (d == d_max) diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 4a5c834c..dd7c2147 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -633,7 +633,11 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m continue; const int s = *shell; if (s >= 0 && s < n_shells) { - acc[s].total_obs++; + // Multiplicity counts what the merge kept. This walk deliberately applies a wider + // filter than the merge so R_meas is computed on the same reflections either way, + // but a masked-ring reflection is not in `unique` and must not be counted against it. + if (!IsMaskedRing(r)) + acc[s].total_obs++; const auto key = generator(r).pack(); const auto mit = merged_I.find(key); if (mit != merged_I.end()) { diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index f3dd6840..34df9a66 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -2017,7 +2017,13 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool if (rnusable[g] == 0) continue; const auto shell = shells.GetShell(acc[g].d); if (!shell || *shell < 0 || *shell >= n_shells) continue; - sa[*shell].total_obs += rnusable[g]; + // Count the MERGED population, not the R_meas one. The R_meas re-walk deliberately + // ignores the ring mask (and, on a search pass, the ice flag), so its count includes + // observations that never entered `unique` - which inflates the reported multiplicity + // of whatever shell they land in. acc[g].nh is what actually went into this group's + // mean, and it is zero for a masked group. (acc[g].d is NaN for such a group, so + // GetShell above already declines it; this is the same statement made where it counts.) + sa[*shell].total_obs += static_cast(acc[g].nh[0] + acc[g].nh[1]); if (std::isfinite(merged_I[g]) && rn[g] > 0) { auto &r = rmeas[g]; r.sum_abs_dev = rabsdev[g]; r.sum_I = rsumI[g]; r.n = rn[g]; r.shell = *shell; @@ -2037,7 +2043,10 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool if (!std::isfinite(I_corr) || !std::isfinite(sigma_corr) || sigma_corr <= 0.0f) continue; const auto shell = shells.GetShell(o.d); if (!shell || *shell < 0 || *shell >= n_shells) continue; - sa[*shell].total_obs++; + // Only what the merge kept counts towards multiplicity - see the GPU branch above. The + // R_meas accumulation below keeps its own, wider filter. + if (usable_merge(o)) + sa[*shell].total_obs++; if (std::isfinite(merged_I[o.group])) { auto &r = rmeas[o.group]; r.sum_abs_dev += std::fabs(static_cast(I_corr) - merged_I[o.group]); diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.cu b/image_analysis/scale_merge/RotationScaleMergeGPU.cu index 676f3243..f02400fb 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.cu +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.cu @@ -558,8 +558,10 @@ namespace { } } - // One thread per group: R_meas accumulators (sum|I_corr - merged_I|, sum_I, n) + usable count for the - // per-shell total_observations. Mirrors MergeAndStats' R_meas re-walk (looser, cell-only filter). + // One thread per group: R_meas accumulators (sum|I_corr - merged_I|, sum_I, n) + the count of + // observations this looser walk accepted. Mirrors MergeAndStats' R_meas re-walk (cell-only filter). + // That count is NOT the per-shell total_observations - it is wider than the merge, so it would + // over-report multiplicity on a masked ring; the host uses it only to skip empty groups. __global__ void MergeRmeasKernel(MergeParams p) { for (int g = blockIdx.x * blockDim.x + threadIdx.x; g < p.n_groups; g += gridDim.x * blockDim.x) { const int lo = p.gstart[g], hi = lo + p.gcount[g]; diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.h b/image_analysis/scale_merge/RotationScaleMergeGPU.h index b5c3d37a..714c44df 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.h +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.h @@ -82,7 +82,8 @@ public: double *swh0, double *swh1, int32_t *nh0, int32_t *nh1, double *d_out, int32_t *rejected, uint8_t *rejected_obs); - // Per-group R_meas accumulators (sum|I_corr-merged_I|, sum_I, n, and the usable count for per-shell + // Per-group R_meas accumulators (sum|I_corr-merged_I|, sum_I, n, and the count this looser walk + // accepted - which the host uses only to skip empty groups, NOT as the per-shell // total_observations); merged_I is uploaded. All arrays length n_groups. void MergeRmeas(const double *merged_I, double *absdev, double *sumI, int32_t *n, int32_t *nusable);