From cbc6a851577a63bc80b2ee96a86de6a18f59cf71 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 22 Jul 2026 22:50:25 +0200 Subject: [PATCH] Debias stills merge with expected-variance weighting 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 . 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) --- common/ScalingSettings.cpp | 9 +++++++++ common/ScalingSettings.h | 11 +++++++++++ image_analysis/scale_merge/Merge.cpp | 23 ++++++++++++++++++++--- image_analysis/scale_merge/Merge.h | 3 ++- rugnux/RugnuxCommandLine.cpp | 2 ++ rugnux/rugnux_cli.cpp | 9 +++++++++ 6 files changed, 53 insertions(+), 4 deletions(-) diff --git a/common/ScalingSettings.cpp b/common/ScalingSettings.cpp index 01477cb9..288a447e 100644 --- a/common/ScalingSettings.cpp +++ b/common/ScalingSettings.cpp @@ -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"); diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h index 74105de7..c541cc61 100644 --- a/common/ScalingSettings.h +++ b/common/ScalingSettings.h @@ -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 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; diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index b67484ef..0b7ec2cb 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -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 &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(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 ). Rebuild the + // signal part at the reflection's EXPECTED intensity instead - decompose out the background part + // and re-add corr* - 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(sigma_corr) * sigma_corr; + if (scaling_settings.GetExpectedVarianceMerge()) { + const double bkg_var = std::max(0.0, a_var - static_cast(image_scale_corr) * I_corr); + const double base = bkg_var + static_cast(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 diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 441252c9..fdcd42c4 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -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 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 diff --git a/rugnux/RugnuxCommandLine.cpp b/rugnux/RugnuxCommandLine.cpp index 5e748279..3a0071a5 100644 --- a/rugnux/RugnuxCommandLine.cpp +++ b/rugnux/RugnuxCommandLine.cpp @@ -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.) diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 42711113..e9fcfbc4 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -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 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 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 space_group_number; std::optional 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);