mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions 4 Packages Projects Releases 67 Wiki Activity

Debias stills merge with expected-variance weighting

Browse Source 
The serial-stills merge (MergeOnTheFly::CorrectedSigma) weighted each
observation by 1/sigma^2 using the observation's OWN sigma. Below ~1
photon the Poisson signal part of that sigma correlates with the
observation's up/down fluctuation, so the inverse-variance mean is
biased low: an up-fluctuated observation acquires a larger sigma and is
over-downweighted. The rotation combine (RotationScaleMerge::
process_rawrun) already avoids this by rebuilding the signal variance at
the pooled estimate; the stills path did not.

Decompose each observation's variance into a background/read part (kept
per-observation) and a Poisson signal part, and rebuild the signal part
at the reflection's expected <I>. Bit-identical when an observation sits
at its reflection mean; only weak-shell weights move. Now default on, so
the stills path matches the rotation path;
--no-expected-variance-merge restores the old observed-sigma weighting.

Validated by paired refinement (phenix, 5 free-set seeds, byte-identical
free flags across arms): R-free-neutral on strong lysozyme and lower
R-free on weak serial-stills data checked against an independent
deposited model (6/6 seeds). The CC1/2 dip on strong data reflects
precision, not accuracy. Applies to both offline rugnux and the online
broker stills merge.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-07-22 22:50:25 +02:00
co-authored by Claude Opus 4.8
parent c52886c8ff
commit cbc6a85157
6 changed files with 53 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files
common/ScalingSettings.cpp
+9
View File
@@ -160,6 +160,15 @@ bool ScalingSettings::GetStillsModulation() const {
return stills_modulation;
}
ScalingSettings &ScalingSettings::ExpectedVarianceMerge(bool input) {
expected_variance_merge = input;
return *this;
}
bool ScalingSettings::GetExpectedVarianceMerge() const {
return expected_variance_merge;
}
ScalingSettings &ScalingSettings::SmoothGDegrees(double input) {
if (input < 0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Smooth-G range must be non-negative");
common/ScalingSettings.h
+11
View File
@@ -51,6 +51,15 @@ class ScalingSettings {
// rotation path fits its own modulation via RotationScaleMerge). Enabled by rugnux --stills-modulation.
bool stills_modulation = false;
// Expected-variance merge weighting for the STILLS merge (MergeOnTheFly). When combining a reflection's
// redundant observations by inverse variance, rebuild the Poisson signal part of each observation's
// variance at the reflection's EXPECTED <I> instead of the observation's own intensity. Weighting by an
// observation's own sigma^2 biases the inverse-variance mean low at <1 photon (an up-fluctuated
// observation gets a larger sigma and is over-downweighted). Default on - it mirrors the rotation combine
// (RotationScaleMerge::process_rawrun), which already does this, and is R-free-neutral on strong data and
// better on weak. --no-expected-variance-merge restores the old observed-sigma weighting.
bool expected_variance_merge = true;
// Smooth the per-frame scale G across frames (centered moving average of log G) before the rot3d
// combine, so a rocking event's partials share a consistent scale. Given as a ROTATION RANGE in
// degrees (like XDS DELPHI), converted to an odd frame window from the oscillation step; this keeps
@@ -92,6 +101,7 @@ public:
ScalingSettings& AbsorptionIter(int input);
ScalingSettings& CorrectionSurfaces(bool input);
ScalingSettings& StillsModulation(bool input);
ScalingSettings& ExpectedVarianceMerge(bool input);
ScalingSettings& SmoothGDegrees(double input);
ScalingSettings& RelativeBDegrees(double input);
@@ -131,6 +141,7 @@ public:
[[nodiscard]] int GetAbsorptionIter() const;
[[nodiscard]] bool GetCorrectionSurfaces() const;
[[nodiscard]] bool GetStillsModulation() const;
[[nodiscard]] bool GetExpectedVarianceMerge() const;
[[nodiscard]] double GetSmoothGDegrees() const;
[[nodiscard]] double GetRelativeBDegrees() const;
image_analysis/scale_merge/Merge.cpp
+20 -3
View File
@@ -93,7 +93,7 @@ void MergeOnTheFly::AddImage(const IntegrationOutcome &outcome, int64_t image_id
continue;
auto hkl = generator(r);
auto hkl_key = hkl.pack();
sigma_corr = CorrectedSigma(I_corr, sigma_corr, hkl_key, r.partiality);
sigma_corr = CorrectedSigma(I_corr, sigma_corr, r.image_scale_corr, hkl_key, r.partiality);
// Robust outlier rejection: drop this observation if it sits more than
// reject_nsigma error-model sigmas from the reflection's median. Needs the active
@@ -262,7 +262,8 @@ double MergeOnTheFly::RefineModulation(std::vector<IntegrationOutcome> &outcomes
return gain / std::max(base, 1e-30); // held-out gain fraction (caller logs)
}
float MergeOnTheFly::CorrectedSigma(float I_corr, float sigma_corr, uint64_t hkl_key, float partiality) const {
float MergeOnTheFly::CorrectedSigma(float I_corr, float sigma_corr, float image_scale_corr,
uint64_t hkl_key, float partiality) const {
if (!error_model_active)
return sigma_corr;
@@ -271,7 +272,23 @@ float MergeOnTheFly::CorrectedSigma(float I_corr, float sigma_corr, uint64_t hkl
const auto it = error_model_mean_I.find(hkl_key);
const double I_for_b = (it != error_model_mean_I.end()) ? it->second : I_corr;
double v = error_model_a * static_cast<double>(sigma_corr) * sigma_corr
// Base variance for the a*sigma^2 term. A weak observation's sigma^2 = a background/read part plus a
// Poisson signal part proportional to its OWN intensity; weighting the merge by 1/sigma^2 with that
// per-observation sigma biases the inverse-variance mean low at <1 photon (an up-fluctuated observation
// gets a larger sigma and is over-downweighted, so the weighted mean drifts below <I>). Rebuild the
// signal part at the reflection's EXPECTED intensity <I> instead - decompose out the background part
// and re-add corr*<I> - so the weight no longer correlates with the observation's own fluctuation. This
// mirrors the rotation combine in RotationScaleMerge::process_rawrun and is bit-identical when the
// observation sits at its reflection mean.
double a_var = static_cast<double>(sigma_corr) * sigma_corr;
if (scaling_settings.GetExpectedVarianceMerge()) {
const double bkg_var = std::max(0.0, a_var - static_cast<double>(image_scale_corr) * I_corr);
const double base = bkg_var + static_cast<double>(image_scale_corr) * std::max(0.0, I_for_b);
if (base > 0.0)
a_var = base;
}
double v = error_model_a * a_var
+ (error_model_b * I_for_b) * (error_model_b * I_for_b);
// Partiality-model uncertainty: a reflection recorded at fraction p carries a systematic intensity
image_analysis/scale_merge/Merge.h
+2 -1
View File
@@ -129,7 +129,8 @@ class MergeOnTheFly {
// observations), so it inflates sigma without biasing the inverse-variance weights -
// using the per-observation I_i instead would over-weight down-fluctuated points.
std::unordered_map<uint64_t, float> error_model_mean_I;
[[nodiscard]] float CorrectedSigma(float I_corr, float sigma_corr, uint64_t hkl_key, float partiality) const;
[[nodiscard]] float CorrectedSigma(float I_corr, float sigma_corr, float image_scale_corr,
uint64_t hkl_key, float partiality) const;
// Optional per-observation outlier rejection: drop observations whose corrected
// intensity lies more than reject_nsigma error-model sigmas from the reflection's
rugnux/RugnuxCommandLine.cpp
+2
View File
@@ -153,6 +153,8 @@ std::string RugnuxCommandLine(const ProcessConfig &config,
add("--partiality-uncertainty", num(sc.GetPartialityUncertaintyCoeff()));
if (sc.GetStillsModulation())
args.emplace_back("--stills-modulation");
if (!sc.GetExpectedVarianceMerge())
args.emplace_back("--no-expected-variance-merge");
// When merging, the CLI skips the large _process.h5 unless asked; emit the flag when it is
// wanted so a copied command matches the GUI's "Save _process.h5" choice. (write_merged has
// no CLI equivalent - the CLI always writes the .mtz/.cif when merging.)
rugnux/rugnux_cli.cpp
+9
View File
@@ -94,6 +94,7 @@ void print_usage() {
std::cout << " --relative-b[=deg] rot3d: fit a per-batch relative-B (beyond the single decay slope) over deg-degree batches; cross-validated (default: 10 deg when bare; off otherwise)" << std::endl;
std::cout << " --no-scaling-corrections rot3d: disable the (default-on) decay + absorption correction surfaces fitted on the fulls after scale-fulls" << std::endl;
std::cout << " --stills-modulation stills: fit a detector-plane modulation (flat-field) surface over the merged reflections (cross-validated; experimental, default off)" << std::endl;
std::cout << " --no-expected-variance-merge stills: disable the default expected-variance merge weighting (which rebuilds each weak observation's signal variance at the reflection mean to de-bias the inverse-variance merge); restores observed-sigma weighting" << std::endl;
std::cout << " -A, --anomalous Anomalous mode (don't merge Friedel pairs)" << std::endl;
std::cout << " -B, --refine-bfactor Refine per image B-factor (stills only)" << std::endl;
std::cout << " --scaling-high-resolution <num> High resolution limit for scaling/merging (manual override; default: no limit)" << std::endl;
@@ -176,6 +177,7 @@ enum {
OPT_RELATIVE_B,
OPT_NO_SCALING_CORRECTIONS,
OPT_STILLS_MODULATION,
OPT_NO_EXPECTED_VARIANCE_MERGE,
OPT_DETECT_ICE_RINGS,
OPT_NO_SCALE_FULLS,
OPT_WRITE_PROCESS_H5,
@@ -223,6 +225,7 @@ static option long_options[] = {
{"relative-b", optional_argument, nullptr, OPT_RELATIVE_B},
{"no-scaling-corrections", no_argument, nullptr, OPT_NO_SCALING_CORRECTIONS},
{"stills-modulation", no_argument, nullptr, OPT_STILLS_MODULATION},
{"no-expected-variance-merge", no_argument, nullptr, OPT_NO_EXPECTED_VARIANCE_MERGE},
{"refine", required_argument, nullptr, 'r'},
{"two-pass-rotation", optional_argument, nullptr, 'R'},
@@ -506,6 +509,7 @@ int main(int argc, char **argv) {
std::optional<double> relative_b_deg_arg; // --relative-b[=deg]; per-batch relative-B width, 0 (off) unless given
bool no_scaling_corrections = false; // --no-scaling-corrections: disable rot3d decay+absorption surfaces
bool stills_modulation_flag = false; // --stills-modulation: detector-plane flat-field surface for stills
bool no_expected_variance_merge = false; // --no-expected-variance-merge: restore observed-sigma stills merge weighting
bool anomalous_mode = false;
std::optional<int64_t> space_group_number;
std::optional<UnitCell> fixed_reference_unit_cell;
@@ -792,6 +796,9 @@ int main(int argc, char **argv) {
case OPT_STILLS_MODULATION:
stills_modulation_flag = true;
break;
case OPT_NO_EXPECTED_VARIANCE_MERGE:
no_expected_variance_merge = true;
break;
case OPT_NO_SCALING_CORRECTIONS:
no_scaling_corrections = true;
break;
@@ -1074,6 +1081,7 @@ int main(int argc, char **argv) {
(experiment.GetGoniometer().has_value() && !force_still) ? 0.7 : 0.0));
scaling_settings.MinCCForImage(min_image_cc / 100.0); // --min-image-cc is percent; the setting is a fraction
scaling_settings.StillsModulation(stills_modulation_flag);
scaling_settings.ExpectedVarianceMerge(!no_expected_variance_merge);
scaling_settings.OutlierRejectNsigma(
outlier_reject_nsigma.value_or(
(experiment.GetGoniometer().has_value() && !force_still) ? REJECT_OUTLIERS_DEFAULT_NSIGMA : 0.0));
@@ -1386,6 +1394,7 @@ int main(int argc, char **argv) {
if (no_scaling_corrections)
scaling_settings.CorrectionSurfaces(false);
scaling_settings.StillsModulation(stills_modulation_flag);
scaling_settings.ExpectedVarianceMerge(!no_expected_variance_merge);
if (d_min_scale_merge)
scaling_settings.HighResolutionLimit_A(d_min_scale_merge.value());
if (resolution_cutoff_method) scaling_settings.ResolutionCutoff(*resolution_cutoff_method);