Wilson outlier test: a Friedel mate does not confirm an improbable observation

The Wilson test keeps an improbable observation whose reflection has other
large observations: they confirm each other, which is how a genuinely strong
reflection survives. A spot of another lattice (satellite crystal, salt or ice
crystallite) that lands on a predicted reflection defeats this: that lattice
diffracts centrosymmetrically too, and half a turn later the rotation brings
-h and the spot's own Friedel mate to the same mirrored detector position, so
the artefact always arrives as a Friedel pair that vouches for itself.

Seen on a P212121 rotation set: one 1.04 A reflection measured three times -
h on one frame and -h half a turn later, both ~5 px off the prediction on the
same detector row, compact single-frame spots of 3000 and 10000 counts (500x
the shell mean), and a symmetry mate elsewhere that reads ~0. The pair
confirmed each other, the median of the three was then an artefact and the
median test deleted the true measurement, and the merged intensity alone took
the last-shell CC1/2 from ~0.18 to 0.80 (XDS: 0.08, both flagged as misfits).

An improbable observation of the exactly opposite indices is now neither a
confirmation nor company; a probable one still testifies against it. The
observation carries a key of its own unreduced indices with key(-h) = -key(h).
Same change on the device path (still bit-identical to the host).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01UqXYdL4qBGYRiSXSGXfCLu

docs/CPU_DATA_ANALYSIS_DECISIONS.md

@@ -76,7 +76,7 @@ Beside them, on rotation data with `--finalist-ledger`, the run reports **twin-i
 
 Merging applies a per-observation median-based $N\sigma$ cut (`--reject-outliers`, default 6σ on rotation, 10σ on stills). On the rotation path the median is weighted by each full's counting variance taken at the reflection's **expected** intensity, as the merge weights are (§10.4): an observation's own Poisson variance falls with its own intensity, so own-variance weights hand the median to whichever equivalents came out low — on a strongly absorbing crystal, where some frames read equivalents ten to a hundred times down, the median sat on the absorbed ones and the test removed correct measurements far above it. The band about the median keeps the counting part linear, $\pm N\sigma_\text{counting}$, but takes the systematic part as a **factor**, $e^{\pm N s}$ about the expected intensity: the error model's systematic term is multiplicative — absorption, illuminated volume, a scale that is off — and as likely to be a factor $1/f$ as $f$, while a band symmetric in $I$ reaches far below the median and only a little above it once that term is large. $s$ is the spread of $\ln(I/\text{median})$ measured on the merge's own fulls (their weighted median of $|\ln(I/\text{median})|$, weighted by $(\langle I\rangle/\sigma_\text{counting})^2$ so the fulls whose ratio counting noise does not blur carry it), with the error model's $b$ as the fallback where it cannot be measured; to first order in $Ns$ the band is the linear one. The same $N\sigma$ cut is fed back into the error model: after an initial $a,b$ fit the parameters are re-fit once on the reflections that survive rejection (dropping any whose squared deviation exceeds $N^2\,[a\,\sigma^2 + (b\,\langle I\rangle)^2]$), so the calibrated errors describe the reflections that actually enter the merge rather than the pre-rejection pool. On **stills** the band is the plain linear one, $\pm N\sqrt{a\,\sigma^2 + (b\,\langle I\rangle)^2}$ about the median: a still's observations of one reflection are partials and per-crystal scales, many of them near zero, so the spread of $\ln(I/\text{median})$ is wide and a band drawn from it let through a powder grain or a neighbour's spot in the box at hundreds of $\sigma$ — on a high-multiplicity run a single such observation was a weak reflection's whole merged intensity, and they took the low-resolution shells' $\mathrm{CC}_{1/2}$ down to 0.6. The error model's $b$ is that spread as it is. $N$ is wider than on rotation because a still's partials scatter with heavier tails than the band's normal core; what has to go sits tens to hundreds of $\sigma$ out.
 
-The median needs three observations, so a reflection measured once or twice — typically one good observation and one artefact (a hot pixel, a zinger) after Friedel merging — is never tested by it, and an artefact that looks precise (thousands of counts, a small relative $\sigma$) can outweigh mates from weak frames and become the median itself. Every observation of the written merge is therefore also judged by **Wilson statistics**. Beyond 4 Å its $E^2 = I/(\varepsilon\,\langle I/\varepsilon\rangle_\mathrm{shell})$ is compared with a bound set by a budget of $\alpha = 0.01$ expected false rejections per dataset: with $N$ observations tested, $\ln(2N/\alpha)$ for acentrics and twice that for centrics, and the observation's lower confidence limit $I - z\sigma$ (with its own $\sigma$, and $z$ from the same budget) must exceed it, so noise in weak shells does not trigger it; a shell whose $\langle I\rangle$ is not itself established at that significance is not judged. The bound is widened by the tail scale of the dataset's own intensity distribution, measured on its well-measured observations (1 for a Wilson crystal; larger under a pseudo-translation or anisotropy). An improbable singleton is dropped; an improbable observation with company only when most of the reflection's other observations are probable and it disagrees with their mean beyond the errors — several large observations confirm each other. An observation whose signal disk lost pixels to the mask or to saturation is never counted as such a witness, since its profile estimate may be low: a strong reflection is not replaced by its clipped mate. The count is `OBSERVATIONS_REJECTED_WILSON=` in the report (included in `OBSERVATIONS_REJECTED=`), and the developer report lists each observation with its image and detector position.
+The median needs three observations, so a reflection measured once or twice — typically one good observation and one artefact (a hot pixel, a zinger) after Friedel merging — is never tested by it, and an artefact that looks precise (thousands of counts, a small relative $\sigma$) can outweigh mates from weak frames and become the median itself. Every observation of the written merge is therefore also judged by **Wilson statistics**. Beyond 4 Å its $E^2 = I/(\varepsilon\,\langle I/\varepsilon\rangle_\mathrm{shell})$ is compared with a bound set by a budget of $\alpha = 0.01$ expected false rejections per dataset: with $N$ observations tested, $\ln(2N/\alpha)$ for acentrics and twice that for centrics, and the observation's lower confidence limit $I - z\sigma$ (with its own $\sigma$, and $z$ from the same budget) must exceed it, so noise in weak shells does not trigger it; a shell whose $\langle I\rangle$ is not itself established at that significance is not judged. The bound is widened by the tail scale of the dataset's own intensity distribution, measured on its well-measured observations (1 for a Wilson crystal; larger under a pseudo-translation or anisotropy). An improbable singleton is dropped; an improbable observation with company only when most of the reflection's other observations are probable and it disagrees with their mean beyond the errors — several large observations confirm each other. A Friedel mate does not confirm: an improbable observation of the exactly opposite indices is not counted as company at all. A spot of another lattice (a satellite crystal, a crystallite of salt or ice) that lands on a predicted reflection is the commonest large observation that is not the crystal's, and that lattice diffracts centrosymmetrically too: half a turn later the rotation brings $-h$ and the spot's own Friedel mate to the same mirrored detector position, so the two always arrive together and would otherwise vouch for each other. A probable Friedel mate still testifies against the observation. An observation whose signal disk lost pixels to the mask or to saturation is never counted as such a witness, since its profile estimate may be low: a strong reflection is not replaced by its clipped mate. The count is `OBSERVATIONS_REJECTED_WILSON=` in the report (included in `OBSERVATIONS_REJECTED=`), and the developer report lists each observation with its image and detector position.
 
 ### 13.4 Automatic resolution cutoff
 

image_analysis/scale_merge/RotationScaleMerge.cpp

@@ -6102,7 +6102,8 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
                     const Obs &o = fulls[wilson_full[j]];
                     const float I_corr = o.I * o.corr;
                     wobs[j] = {I_corr, model_sigma(o, I_corr, o.sigma * o.corr), o.d, group_epsilon[g],
-                               group_centric[g] != 0, o.clipped != 0, merge_friedel ? g : pair_of_group[g]};
+                               group_centric[g] != 0, o.clipped != 0, merge_friedel ? g : pair_of_group[g],
+                               HKLFriedelKey(o.h, o.k, o.l)};
                 }
             });
             // On the device where the merge is: the same test to the bit (WilsonOutliersGPU.h).

image_analysis/scale_merge/WilsonOutliers.cpp

@@ -189,7 +189,7 @@ WilsonOutlierResult WilsonOutliers(const std::vector<WilsonObservation> &obs, do
     // An improbable observation measured once is rejected. One with company is rejected when most of
     // its reflection's other observations are probable, unclipped witnesses and it disagrees with their
     // mean beyond the errors - for a pair, the other member. Where most are large, the reflection is
-    // (Aimless's "keep if most observations are large").
+    // (Aimless's "keep if most observations are large"). An improbable Friedel mate is not company.
     // Each reflection is decided on its own observations alone, so the reflections go side by side:
     // where each one's stretch of `tested` starts, then the stretches in blocks.
     const int n_tested = static_cast<int>(tested.size());
@@ -210,16 +210,18 @@ WilsonOutlierResult WilsonOutliers(const std::vector<WilsonObservation> &obs, do
             for (size_t q = lo; q < hi; ++q) {
                 const size_t i = tested[q];
                 if (!improbable[i]) continue;
-                size_t n_probable = 0;
+                const int64_t mate = -obs[i].hkl;   // its Friedel mate (see WilsonOutliers.h)
+                size_t n_probable = 0, n_others = 0;
                 double sw = 0.0, swI = 0.0;
                 for (size_t r = lo; r < hi; ++r) {
                     const auto &o = obs[tested[r]];
-                    if (r == q || improbable[tested[r]] || o.clipped) continue;
+                    if (r == q || (improbable[tested[r]] && mate != 0 && o.hkl == mate)) continue;
+                    ++n_others;
+                    if (improbable[tested[r]] || o.clipped) continue;
                     const double w = 1.0 / (static_cast<double>(o.sigma) * o.sigma);
                     sw += w; swI += w * o.I;
                     ++n_probable;
                 }
-                const size_t n_others = hi - lo - 1;
                 bool reject = n_others == 0;
                 if (2 * n_probable > n_others) {
                     const auto &o = obs[i];

image_analysis/scale_merge/WilsonOutliers.h

@@ -36,6 +36,14 @@
 // kept, which is how a genuinely strong reflection survives. An observation whose signal disk lost
 // pixels to the mask or to saturation is not counted as a probable witness: its profile estimate of
 // what remained may be low, and a strong reflection must not be replaced by its clipped mate.
+//
+// Confirmation needs a second, independent measurement, and a Friedel mate is not one against the
+// commonest cause of a large observation that is not the crystal's: a spot of another lattice (a
+// satellite crystal, a crystallite of salt or ice) that lands on a predicted reflection. That lattice
+// is centrosymmetric in diffraction too, and half a turn later the rotation brings both -h and the
+// spot's own Friedel mate to the same mirrored place on the detector, so the two always come together.
+// An improbable observation of the exact opposite indices therefore neither confirms the observation
+// nor counts among its company; a probable one still testifies against it.
 struct WilsonObservation {
     float I;          // scaled intensity
     float sigma;      // its own standard uncertainty, at its own intensity
@@ -44,8 +52,15 @@ struct WilsonObservation {
     bool centric;
     bool clipped;     // lost pixels to the mask or to saturation: may be low, so never a witness
     int32_t unit;     // the reflection the observation measures; < 0 = not part of the test
+    int64_t hkl = 0;  // HKLFriedelKey of its own (unreduced) indices; 0 = not known
 };
 
+// A key of the indices with key(-h) = -key(h), so a Friedel mate is the observation of the negated key.
+inline int64_t HKLFriedelKey(int32_t h, int32_t k, int32_t l) {
+    constexpr int64_t S = int64_t{1} << 20;
+    return (static_cast<int64_t>(h) * S + k) * S + l;
+}
+
 struct WilsonOutlierResult {
     std::vector<uint8_t> rejected;   // per observation
     std::vector<float> e2;           // I / (epsilon <I/epsilon>_shell); NaN outside the tested range

image_analysis/scale_merge/WilsonOutliersGPU.cu

@@ -177,17 +177,19 @@ namespace {
         for (int q = lo; q < hi; ++q) {
             const int i = tested[q];
             if (!improbable[i]) continue;
-            size_t n_probable = 0;
+            const int64_t mate = -obs[i].hkl;
+            size_t n_probable = 0, n_others = 0;
             double sw = 0.0, swI = 0.0;
             for (int r = lo; r < hi; ++r) {
                 const WilsonObservation &o = obs[tested[r]];
-                if (r == q || improbable[tested[r]] || o.clipped) continue;
+                if (r == q || (improbable[tested[r]] && mate != 0 && o.hkl == mate)) continue;
+                ++n_others;
+                if (improbable[tested[r]] || o.clipped) continue;
                 const double w = 1.0 / __dmul_rn(double(o.sigma), double(o.sigma));
                 sw = __dadd_rn(sw, w);
                 swI = __dadd_rn(swI, __dmul_rn(w, double(o.I)));
                 ++n_probable;
             }
-            const size_t n_others = hi - lo - 1;
             bool reject = n_others == 0;
             if (2 * n_probable > n_others) {
                 const WilsonObservation &o = obs[i];

tests/WilsonOutliersTest.cpp

@@ -66,6 +66,27 @@ TEST_CASE("WilsonOutliers: two large observations of one reflection confirm each
     CHECK(r.n_rejected == 0);
 }
 
+TEST_CASE("WilsonOutliers: a Friedel pair of large observations does not confirm itself", "[wilson_outliers]") {
+    auto obs = WilsonPopulation(10000, 2, 1000.0, 100.0, 2);
+    // h and -h, both high - another lattice's spot and its own Friedel mate - and a symmetry mate of h
+    // (other indices, measured elsewhere) that reads the crystal's ordinary intensity.
+    const int32_t unit = 10000;
+    obs.push_back({60000.0f, 300.0f, 2.5f, 1.0f, false, false, unit, HKLFriedelKey(5, -3, 7)});
+    obs.push_back({61000.0f, 300.0f, 2.5f, 1.0f, false, false, unit, HKLFriedelKey(-5, 3, -7)});
+    obs.push_back({900.0f, 40.0f, 2.5f, 1.0f, false, false, unit, HKLFriedelKey(-5, 3, 7)});
+    auto r = WilsonOutliers(obs, 0.01);
+    CHECK(r.rejected[obs.size() - 3] == 1);
+    CHECK(r.rejected[obs.size() - 2] == 1);
+    CHECK(r.rejected[obs.size() - 1] == 0);
+    CHECK(r.n_rejected == 2);
+
+    // The same two large observations under indices that are not each other's negation (h and a
+    // symmetry mate) confirm each other, as before.
+    obs[obs.size() - 2].hkl = HKLFriedelKey(5, -3, -7);
+    r = WilsonOutliers(obs, 0.01);
+    CHECK(r.n_rejected == 0);
+}
+
 TEST_CASE("WilsonOutliers: the symmetry enhancement factor keeps an axial reflection", "[wilson_outliers]") {
     auto obs = WilsonPopulation(10000, 2, 1000.0, 100.0, 3);
     // Measured once, 12x <I/epsilon> at epsilon 4: 48x the shell mean, but only E^2 = 12.
@@ -166,7 +187,7 @@ TEST_CASE("WilsonOutliers: the device test is the host test", "[wilson_outliers]
                 double I = I_true + sigma * gauss(rng);
                 if (uni(rng) < 0.002) I *= 200.0;
                 obs.push_back({static_cast<float>(I), static_cast<float>(sigma), d, epsilon, centric,
-                               uni(rng) < 0.05, u});
+                               uni(rng) < 0.05, u, HKLFriedelKey(m % 2 ? -1 : 1, m % 2 ? -2 : 2, m % 3 ? 3 : -3)});
             }
         }
         obs[3].unit = -1;