diff --git a/docs/CPU_DATA_ANALYSIS_DECISIONS.md b/docs/CPU_DATA_ANALYSIS_DECISIONS.md index dc3f1c73b..ed5c56e79 100644 --- a/docs/CPU_DATA_ANALYSIS_DECISIONS.md +++ b/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 diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 2ac99041c..79a02bd2c 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/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). diff --git a/image_analysis/scale_merge/WilsonOutliers.cpp b/image_analysis/scale_merge/WilsonOutliers.cpp index ef2c9e83c..909d1c6c8 100644 --- a/image_analysis/scale_merge/WilsonOutliers.cpp +++ b/image_analysis/scale_merge/WilsonOutliers.cpp @@ -189,7 +189,7 @@ WilsonOutlierResult WilsonOutliers(const std::vector &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(tested.size()); @@ -210,16 +210,18 @@ WilsonOutlierResult WilsonOutliers(const std::vector &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(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]; diff --git a/image_analysis/scale_merge/WilsonOutliers.h b/image_analysis/scale_merge/WilsonOutliers.h index 03e33872f..317f4b1e5 100644 --- a/image_analysis/scale_merge/WilsonOutliers.h +++ b/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(h) * S + k) * S + l; +} + struct WilsonOutlierResult { std::vector rejected; // per observation std::vector e2; // I / (epsilon _shell); NaN outside the tested range diff --git a/image_analysis/scale_merge/WilsonOutliersGPU.cu b/image_analysis/scale_merge/WilsonOutliersGPU.cu index f3fd947e7..789ff2a46 100644 --- a/image_analysis/scale_merge/WilsonOutliersGPU.cu +++ b/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]; diff --git a/tests/WilsonOutliersTest.cpp b/tests/WilsonOutliersTest.cpp index b13052ccc..38f9d963b 100644 --- a/tests/WilsonOutliersTest.cpp +++ b/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 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(I), static_cast(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;