Serial stills: an observation with a negative background variance is no measurement
A frame whose background region reads below zero - a detector module whose pedestal is off on that
frame - integrates reflections with a negative background variance, a sigma floored at one photon and a
net intensity of thousands (I/sigma ~8000). The unscaled merge the per-image scales are fitted against
weights such an observation ~100x above any real one, so it takes its reflection over; every frame that
sees that reflection and finds nothing there is then pulled towards a scale of 0 (the proportional model
error caps such an observation's normalised residual at ~3.3 sigma, so the Cauchy fit never rejects it
and its pull grows as G falls), lands under the 1/50 floor and is dropped, and the partiality refit
against the survivors drops more.
On a serial run whose detector shows this on one module every 324th frame, always on the same residue
of the frame number modulo 4, two same-arm quarters of the run (frames mod 4 = 0 and = 2) merged to
CC1/2 0.78 and 0.28. The collapsed frames were good frames (CC with the deposited model's Fcalc 0.29
against 0.34 for the kept ones, the same hand) whose CC with the scaling reference was 0.03: the
reference, not the frames, was wrong. Measured (frames mod 4 = 0 / = 2, run's geometry):
before: 446 / 1355 frames dropped, R_split 32.6 / 64.0%, <I^2>/<I>^2 2.08 / 2.74, <|L|> 0.475 / 0.788,
CC of the merge with the model Fcalc (d >= 2.5 A) 0.842 / 0.381
after: 446 / 412 dropped, R_split 32.6 / 35.9%, 2.08 / 2.17, 0.475 / 0.478, 0.842 / 0.816
The whole run: 9081 -> 3081 frames dropped, CC1/2 0.59 -> 0.955, R_split 18.6 -> 14.9%, Fcalc CC
0.879 -> 0.884 (d >= 2.5 A); two runs without the defect are unchanged.
AcceptReflection now refuses such an observation (a NaN variance, from older files, is kept), so it
leaves every stills scale fit, merge and ambiguity correlation, and the stills merge logs how many
observations on how many frames were left out.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK
(cherry picked from commit bd5c59915fd7698afcc0ef81722e562c78c5ee3a)
This commit is contained in:
1 parent
97395c6d1a
commit
ced02f1f78
4 files changed
+75
No files matched your search
@@ -343,6 +343,8 @@ The partiality applied is fixed by the data type and scaling stage, not chosen f
|
||||
|
||||
Reflections below a minimum partiality are left out of the merge (`--min-partiality`, default 0.02 on rotation and 0.5 on stills): the empirical rocking curve of a still has a floor far out on its tail, from mis-oriented crystals, and dividing such an observation by a small model $p$ makes it biased, not merely noisy. On stills this holds for the merges the scale and the partiality are refined against. The final merge, which has the fitted error model, takes every observation: the model's partiality term $c\,((1-p)/p)\,\langle I\rangle$ (§10.5) weights a low-partiality observation by what it is worth. A reflection - each Bijvoet hand under `-A` - is reported only if at least one of its observations reaches $p \ge 0.2$; one seen only far out on its rocking curve rests on the partiality model alone. The error model is fitted on the observations at $p \ge 0.2$.
|
||||
|
||||
An observation whose background variance is below zero is left out of the stills scale fits and merges, and the log counts them (`Scaling: N observations on M frames are left out`). Pixels counting photons cannot give one: it is what a background region reading below zero leaves, as when a detector module's pedestal is off on a frame, and integration then floors the reflection's $\sigma$ at one photon while its net intensity, measured against the negative background, reaches thousands. Such an observation weighs a hundred times more than any real one, takes its reflection over in an unscaled merge, and every frame that sees that reflection and finds nothing there is fitted towards a scale of zero, below the floor of the per-frame scales (§10.3), and dropped. On a serial run where one module failed this way every 324th frame, always in the same residue of the frame number modulo 4, the frames of that residue lost about half of their number to the floor and those of the other residues next to none, so two halves of one run merged to CC1/2 0.28 and 0.78.
|
||||
|
||||
### 10.3 Smoothing of per-frame scales
|
||||
|
||||
The per-frame scales $G_i$ are fit by robust (Cauchy) inverse-variance-weighted ratios; there is no explicit $G\approx1$ prior. For rotation datasets, optional smoothing enforces the expectation that scale and mosaicity vary slowly across a sweep: **after** the per-frame fit, $\log G_i$ (and the mosaicity) are replaced by a centred **moving average** over a window spanning a configurable rotation range (XDS DELPHI-like; `--smooth-g`, default 5° for rot3d, off otherwise). It is a post-fit smoothing pass, not a curvature penalty inside the least-squares objective. (The per-frame scale refitted on the combined fulls, §10.6, is smoothed differently: by a penalised smoother whose smoothness is chosen by cross-validation.)
|
||||
|
||||
@@ -96,6 +96,13 @@ bool AcceptReflection(const Reflection &r, std::optional<double> d_min_limit, st
|
||||
return false;
|
||||
if (!std::isfinite(r.sigma) || r.sigma <= 0.0)
|
||||
return false;
|
||||
// A background variance below zero is no measurement either. It is what a background region reading
|
||||
// below zero leaves - pixels that are not counting photons, such as a module whose pedestal is off on
|
||||
// that frame - and the sigma of such a reflection falls to its floor of a photon while its net intensity
|
||||
// reaches thousands: an observation a hundred times more trusted than any real one, which takes its
|
||||
// reflection over in every merge and drags every frame that sees that reflection towards a scale of 0.
|
||||
if (r.var_bkg < 0.0f)
|
||||
return false;
|
||||
return true;
|
||||
}
|
||||
|
||||
@@ -113,6 +120,8 @@ bool AcceptReflection(const Reflection &r, double d_min_limit, double d_max_limi
|
||||
return false;
|
||||
if (!std::isfinite(r.sigma) || r.sigma <= 0.0)
|
||||
return false;
|
||||
if (r.var_bkg < 0.0f) // see above
|
||||
return false;
|
||||
return true;
|
||||
}
|
||||
|
||||
@@ -708,6 +708,22 @@ OtherCellForm LeaveOutMinorCellForm(std::vector<IntegrationOutcome> &outcomes, c
|
||||
return {forms.minor_frames, cells.size(), forms.minor_cell, forms.minor_frames};
|
||||
}
|
||||
|
||||
// The observations AcceptReflection refuses for a background variance below zero (see there), counted so the
|
||||
// log says how many there were and on how many frames: a module that stopped counting on a frame gives dozens.
|
||||
static void LogNegativeVariance(const std::vector<IntegrationOutcome> &outcomes, Logger &logger) {
|
||||
size_t observations = 0, frames = 0;
|
||||
for (const auto &o : outcomes) {
|
||||
size_t n = 0;
|
||||
for (const auto &r : o.reflections)
|
||||
n += r.var_bkg < 0.0f;
|
||||
observations += n;
|
||||
frames += n > 0;
|
||||
}
|
||||
if (observations > 0)
|
||||
logger.Info("Scaling: {} observations on {} frames are left out - their background variance is below zero, "
|
||||
"which pixels counting photons cannot give", observations, frames);
|
||||
}
|
||||
|
||||
StillsMergeResult StillsScaleMerge(const DiffractionExperiment &experiment, std::vector<IntegrationOutcome> &outcomes,
|
||||
const std::optional<UnitCell> &reference_cell,
|
||||
const std::optional<StillsRunWidth> &run_width,
|
||||
@@ -733,6 +749,8 @@ StillsMergeResult StillsScaleMerge(const DiffractionExperiment &experiment, std:
|
||||
for (const auto &r : o.reflections)
|
||||
f.corr.push_back(r.image_scale_corr);
|
||||
}
|
||||
if (!for_search)
|
||||
LogNegativeVariance(outcomes, logger);
|
||||
{
|
||||
phase("Scaling images");
|
||||
TimingMark scale_mark(logger, "stills scaling images");
|
||||
|
||||
@@ -127,6 +127,14 @@ TEST_CASE("AcceptReflection_ResolutionLimits") {
|
||||
r.overloaded = true;
|
||||
CHECK_FALSE(AcceptReflection(r, std::nullopt, std::nullopt));
|
||||
CHECK_FALSE(AcceptReflection(r, 0.0, 0.0));
|
||||
|
||||
// Nor is one whose background variance is below zero; an unknown one (NaN, older files) is kept.
|
||||
r.overloaded = false;
|
||||
r.var_bkg = -1.0f;
|
||||
CHECK_FALSE(AcceptReflection(r, std::nullopt, std::nullopt));
|
||||
CHECK_FALSE(AcceptReflection(r, 0.0, 0.0));
|
||||
r.var_bkg = NAN;
|
||||
CHECK(AcceptReflection(r, std::nullopt, std::nullopt));
|
||||
}
|
||||
|
||||
// --- Completeness denominator --------------------------------------------------------------------
|
||||
@@ -1559,6 +1567,44 @@ TEST_CASE("ScaleOnTheFly_TooFewReflectionsLeavesTheImageOut") {
|
||||
CHECK(std::isnan(r.image_scale_corr));
|
||||
}
|
||||
|
||||
// A frame whose background region reads below zero (a module not counting photons) integrates reflections with
|
||||
// a negative background variance, a sigma at its floor of one photon and an intensity of thousands. Merged,
|
||||
// such an observation takes its reflection over, and every frame that sees the reflection - and finds
|
||||
// nothing there - is then fitted towards a scale of 0. Left out, the frames keep their own scales.
|
||||
TEST_CASE("ScaleOnTheFly_NegativeBackgroundVarianceDoesNotCollapseTheScales") {
|
||||
DiffractionExperiment x;
|
||||
x.SetSpaceGroup(gemmi::get_spacegroup_by_name("P 1"));
|
||||
constexpr int frames = 20, n = 60;
|
||||
std::vector<IntegrationOutcome> outcomes(frames);
|
||||
for (int i = 0; i < frames; ++i) {
|
||||
const float G = 1.0f + 0.1f * static_cast<float>(i % 5);
|
||||
for (int h = 1; h <= n; ++h) {
|
||||
Reflection r{};
|
||||
r.h = h; r.k = 1; r.l = 1;
|
||||
r.d = 3.0f;
|
||||
// the first ten are beyond the crystal's diffraction: nothing there
|
||||
r.I = h <= 10 ? 0.0f : G * 100.0f * static_cast<float>(1 + h % 7);
|
||||
r.var_bkg = 50.0f;
|
||||
r.sigma = std::sqrt(r.I + r.var_bkg);
|
||||
r.partiality = 1.0f;
|
||||
r.prescaling_corr = 1.0f;
|
||||
r.image_scale_corr = 1.0f;
|
||||
if (i == 0 && h <= 10) { // the broken module of frame 0
|
||||
r.I = 1e5f;
|
||||
r.var_bkg = -1e4f;
|
||||
r.sigma = 1.0f;
|
||||
}
|
||||
outcomes[i].reflections.push_back(r);
|
||||
}
|
||||
}
|
||||
ScaleOnTheFly(x, MergeAll(x, outcomes, 1, std::nullopt)).Scale(outcomes, 1);
|
||||
for (int i = 1; i < frames; ++i) {
|
||||
REQUIRE(outcomes[i].image_scale_g.has_value());
|
||||
CHECK(*outcomes[i].image_scale_g / *outcomes[1].image_scale_g
|
||||
== Catch::Approx((1.0 + 0.1 * (i % 5)) / 1.1).epsilon(0.01));
|
||||
}
|
||||
}
|
||||
|
||||
// Stills with the Friedel hands kept apart report CC_anom (between frame-whole half-sets) and SigAno: a
|
||||
// consistent anomalous difference far above the noise correlates close to 1.
|
||||
TEST_CASE("MergeOnTheFly_StillsAnomalousStatistics") {
|
||||
|
||||
Reference in new issue
Block a user