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

rugnux: add --partiality-uncertainty stills merge term (partiality-model error)

Browse Source 
A stills reflection recorded at partiality p carries a systematic intensity error ~(dp/p) that is
proportional to <I> 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*<I>*(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 <c> (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) <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-07-14 19:52:46 +02:00
co-authored by Claude Opus 4.8
parent 2f5ed411e2
commit 4c5da20eda
5 changed files with 38 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files
common/ScalingSettings.cpp
+9
View File
@@ -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,
common/ScalingSettings.h
+3
View File
@@ -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<double> 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<double> GetForcedMosaicity() const;
[[nodiscard]] double GetCaptureUncertaintyCoeff() const;
[[nodiscard]] double GetPartialityUncertaintyCoeff() const;
[[nodiscard]] double GetMinCapturedFraction() const;
[[nodiscard]] double GetMinCCForImage() const;
[[nodiscard]] double GetOutlierRejectNsigma() const;
image_analysis/scale_merge/Merge.cpp
+17 -4
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);
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<double>(sigma_corr) * sigma_corr
+ (error_model_b * I_for_b) * (error_model_b * I_for_b);
double v = error_model_a * static_cast<double>(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 <I> 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<double>(partiality)));
const double t = c * I_for_b * one_minus_p;
v += t * t;
}
return (v > 0.0) ? static_cast<float>(std::sqrt(v)) : sigma_corr;
}
image_analysis/scale_merge/Merge.h
+1 -1
View File
@@ -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<uint64_t, float> 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
rugnux/rugnux_cli.cpp
+8
View File
@@ -94,6 +94,7 @@ void print_usage() {
std::cout << " --resolution-shells <num> Number of resolution shells in the reported statistics table (default: 10)" << std::endl;
std::cout << " --min-partiality <num> Minimum partiality to accept reflection (default: 0.02)" << std::endl;
std::cout << " --capture-uncertainty <num> 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 <num> stills: extra merge sigma ~num*(1-partiality)*<I> 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 <num> 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 <num> Diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed" << std::endl;
std::cout << " --reject-outliers <num> 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<double> min_captured_fraction_arg; // explicit --min-captured-fraction; default depends on rotation
std::optional<double> capture_uncertainty_arg; // explicit --capture-uncertainty; default depends on rot3d
std::optional<double> partiality_uncertainty_arg; // --partiality-uncertainty (stills partiality merge sigma)
std::optional<double> 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(