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<int32_t>(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 <I/sigma> (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) <noreply@anthropic.com>
This commit is contained in:
2026-08-06 19:06:00 +02:00
co-authored by Claude Opus 5
parent 17eb6ef091
commit 06b8c8ed66
5 changed files with 26 additions and 7 deletions
+4 -1
View File
@@ -23,7 +23,10 @@ ResolutionShells::ResolutionShells(float d_min, float d_max, int32_t nshells)
}
std::optional<int32_t> 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<int32_t>(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)