From 4c5da20eda87662deebd5ed2bd587d5427204948 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 14 Jul 2026 19:52:46 +0200 Subject: [PATCH] rugnux: add --partiality-uncertainty stills merge term (partiality-model error) A stills reflection recorded at partiality p carries a systematic intensity error ~(dp/p) that is proportional to and grows as p falls; plain counting sigma misses it, so strong low-p partials are over-trusted in the merge. Add sigma^2 += (c**(1-p))^2 in MergeOnTheFly::CorrectedSigma - the stills-partiality analog of the existing rotation --capture-uncertainty term ((1-captured_fraction)*I in RotationScaleMerge). Inert when partiality==1 (no --still-partiality), and gated on a real systematic (error_model_b > 1, i.e. ISa < 1) so it fires on strong/medium stills but auto-skips weak counting-limited data where it would only over-concentrate the merge and hurt. Opt-in via --partiality-uncertainty (default 0; ~2.5 recommended with --still-partiality). Prototyped (Python replica of the merge): lyso8 CC1/2 +1.7, CCref +4.8, R-free -0.023; LOV CCref +2; harmful on weak OCP (hence the b>1 gate). Full in-binary validation across all serial targets pending. Co-Authored-By: Claude Opus 4.8 (1M context) --- common/ScalingSettings.cpp | 9 +++++++++ common/ScalingSettings.h | 3 +++ image_analysis/scale_merge/Merge.cpp | 21 +++++++++++++++++---- image_analysis/scale_merge/Merge.h | 2 +- rugnux/rugnux_cli.cpp | 8 ++++++++ 5 files changed, 38 insertions(+), 5 deletions(-) diff --git a/common/ScalingSettings.cpp b/common/ScalingSettings.cpp index f3341e1d..c57ff626 100644 --- a/common/ScalingSettings.cpp +++ b/common/ScalingSettings.cpp @@ -191,6 +191,15 @@ double ScalingSettings::GetCaptureUncertaintyCoeff() const { return capture_uncertainty_coeff; } +ScalingSettings &ScalingSettings::PartialityUncertaintyCoeff(double input) { + partiality_uncertainty_coeff = input; + return *this; +} + +double ScalingSettings::GetPartialityUncertaintyCoeff() const { + return partiality_uncertainty_coeff; +} + ScalingSettings &ScalingSettings::MinCapturedFraction(double input) { if (input < 0.0) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h index 1056bb1b..9995ecf2 100644 --- a/common/ScalingSettings.h +++ b/common/ScalingSettings.h @@ -29,6 +29,7 @@ class ScalingSettings { // a fraction f<1 of its rocking curve is extrapolated, and the unobserved (1-f) carries a // systematic error ~coeff*(1-f)*I that plain counting sigma misses. 0 = off (baseline). double capture_uncertainty_coeff = 0.0; + double partiality_uncertainty_coeff = 0.0; // Full-level captured-fraction floor for the rot3d combine: drop a reconstructed full whose rocking // curve was only fractionally captured (sum of its partials' partiality < this). Distinct from // min_partiality, which gates individual partials; this gates the assembled full. 0 = off (baseline). @@ -75,6 +76,7 @@ public: ScalingSettings& MinPartiality(double min_partiality); ScalingSettings& ForcedMosaicity(std::optional input); ScalingSettings& CaptureUncertaintyCoeff(double input); + ScalingSettings& PartialityUncertaintyCoeff(double input); ScalingSettings& MinCapturedFraction(double input); ScalingSettings& MinCCForImage(double min_cc_for_image); ScalingSettings& OutlierRejectNsigma(double input); @@ -112,6 +114,7 @@ public: [[nodiscard]] double GetMinPartiality() const; [[nodiscard]] std::optional GetForcedMosaicity() const; [[nodiscard]] double GetCaptureUncertaintyCoeff() const; + [[nodiscard]] double GetPartialityUncertaintyCoeff() const; [[nodiscard]] double GetMinCapturedFraction() const; [[nodiscard]] double GetMinCCForImage() const; [[nodiscard]] double GetOutlierRejectNsigma() const; diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 405a9b26..916e72e9 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); + sigma_corr = CorrectedSigma(I_corr, sigma_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 @@ -131,7 +131,7 @@ void MergeOnTheFly::AddImage(const IntegrationOutcome &outcome, int64_t image_id } } -float MergeOnTheFly::CorrectedSigma(float I_corr, float sigma_corr, uint64_t hkl_key) const { +float MergeOnTheFly::CorrectedSigma(float I_corr, float sigma_corr, uint64_t hkl_key, float partiality) const { if (!error_model_active) return sigma_corr; @@ -140,8 +140,21 @@ 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; - const double v = error_model_a * static_cast(sigma_corr) * sigma_corr - + (error_model_b * I_for_b) * (error_model_b * I_for_b); + double v = error_model_a * static_cast(sigma_corr) * sigma_corr + + (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 + // error ~ (dp/p) that is proportional to and grows as p falls - plain counting sigma misses it, + // so strong low-p partials would otherwise be over-trusted. This is the stills-partiality analog of + // the rotation --capture-uncertainty term ((1-captured_fraction)*I in RotationScaleMerge). Inert when + // partiality == 1 (no stills partiality model). Gated on a real systematic (error_model_b > 1, i.e. + // ISa < 1): on weak counting-limited data (small b) it would only over-concentrate the merge and hurt. + const double c = scaling_settings.GetPartialityUncertaintyCoeff(); + if (c > 0.0 && error_model_b > 1.0) { + const double one_minus_p = std::max(0.0, std::min(1.0, 1.0 - static_cast(partiality))); + const double t = c * I_for_b * one_minus_p; + v += t * t; + } return (v > 0.0) ? static_cast(std::sqrt(v)) : sigma_corr; } diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 9c0ee572..32f5c58a 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -109,7 +109,7 @@ 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) const; + [[nodiscard]] float CorrectedSigma(float I_corr, float sigma_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/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index a31bdaa8..2ce2dab7 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -94,6 +94,7 @@ void print_usage() { std::cout << " --resolution-shells Number of resolution shells in the reported statistics table (default: 10)" << std::endl; std::cout << " --min-partiality Minimum partiality to accept reflection (default: 0.02)" << std::endl; std::cout << " --capture-uncertainty rot3d: systematic sigma ~num*(1-captured_fraction)*I on under-captured fulls (default: 1.0 for rot3d, 0 otherwise)" << std::endl; + std::cout << " --partiality-uncertainty stills: extra merge sigma ~num*(1-partiality)* on partials (use with --still-partiality; auto-gated to error-model b>1 / ISa<1; default 0, ~2.5 recommended)" << std::endl; std::cout << " --min-captured-fraction rot3d: drop a combined full whose rocking curve was captured below this fraction (edge-of-sweep truncated fulls) (default: 0.7 for rotation, 0 otherwise; 0 = off)" << std::endl; std::cout << " --mosaicity Diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed" << std::endl; std::cout << " --reject-outliers Per-observation merge outlier rejection, N sigma from the per-reflection median (default: 6 for rot3d, XDS/DIALS-style; 0 = off)" << std::endl; @@ -159,6 +160,7 @@ enum { OPT_STILL_PARTIALITY, OPT_SCALE_FULLS, OPT_CAPTURE_UNCERTAINTY, + OPT_PARTIALITY_UNCERTAINTY, OPT_MIN_CAPTURED_FRACTION, OPT_MOSAICITY, OPT_SMOOTH_G, @@ -239,6 +241,7 @@ static option long_options[] = { {"max-spots", required_argument, nullptr, OPT_MAX_SPOTS}, {"min-partiality", required_argument, nullptr, OPT_MIN_PARTIALITY}, {"capture-uncertainty", required_argument, nullptr, OPT_CAPTURE_UNCERTAINTY}, + {"partiality-uncertainty", required_argument, nullptr, OPT_PARTIALITY_UNCERTAINTY}, {"min-captured-fraction", required_argument, nullptr, OPT_MIN_CAPTURED_FRACTION}, {"mosaicity", required_argument, nullptr, OPT_MOSAICITY}, {"min-image-cc", required_argument, nullptr, OPT_MIN_IMAGE_CC}, @@ -475,6 +478,7 @@ int main(int argc, char **argv) { double min_partiality = 0.02; std::optional min_captured_fraction_arg; // explicit --min-captured-fraction; default depends on rotation std::optional capture_uncertainty_arg; // explicit --capture-uncertainty; default depends on rot3d + std::optional partiality_uncertainty_arg; // --partiality-uncertainty (stills partiality merge sigma) std::optional forced_mosaicity_arg; // diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed double min_image_cc = 0.0; int64_t scaling_iter = 3; @@ -730,6 +734,9 @@ int main(int argc, char **argv) { case OPT_CAPTURE_UNCERTAINTY: capture_uncertainty_arg = parse_double_arg(optarg, "--capture-uncertainty", logger); break; + case OPT_PARTIALITY_UNCERTAINTY: + partiality_uncertainty_arg = parse_double_arg(optarg, "--partiality-uncertainty", logger); + break; case OPT_MIN_CAPTURED_FRACTION: min_captured_fraction_arg = parse_double_arg(optarg, "--min-captured-fraction", logger); break; @@ -1306,6 +1313,7 @@ int main(int argc, char **argv) { // over-extrapolated under-captured fulls and, with the mosaicity fix, lifts rotation ISa/anomalous // substantially. Off for non-rot3d (no combine). An explicit --capture-uncertainty always wins. scaling_settings.CaptureUncertaintyCoeff(capture_uncertainty_arg.value_or(rotation_indexing ? 1.0 : 0.0)); + scaling_settings.PartialityUncertaintyCoeff(partiality_uncertainty_arg.value_or(0.0)); scaling_settings.ForcedMosaicity(forced_mosaicity_arg); scaling_settings.MinCCForImage(min_image_cc / 100.0); // --min-image-cc is in percent; the setting is a fraction scaling_settings.OutlierRejectNsigma(