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

Merge: make outlier rejection a ScalingSettings parameter + CLI flag, default 6

Browse Source 
Promote the per-observation merge outlier rejection from the temporary PR_REJECT
env knob to a real setting: ScalingSettings::OutlierRejectNsigma (default 6, <=0
disables), driven into MergeOnTheFly via the constructor, with a --reject-outliers
CLI flag in jfjoch_process. Default-on at 6 sigma matches XDS (MISFITS/REJECT) and
DIALS (normalised-deviation test); it validated on the full-jet R-free (0.2625 ->
0.2585). Applies to both the PixelRefine and classical merge paths.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-06-15 12:42:25 +02:00
co-authored by Claude Opus 4.8
parent a47b376dc3
commit 6f733d74c2
5 changed files with 30 additions and 12 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tools/jfjoch_process.cpp
+12 -7
View File
@@ -72,6 +72,7 @@ void print_usage() {
std::cout << " -w, --wedge[=num] Refine image wedge during scaling with starting wedge value" << std::endl;
std::cout << " --scaling-high-resolution <num> High resolution limit for spot finding (default: no limit)" << std::endl;
std::cout << " --min-partiality <num> Minimum partiality to accept reflection (default: 0.02)" << std::endl;
std::cout << " --reject-outliers <num> Per-observation merge outlier rejection, N sigma from the per-reflection median (default: 6; 0 disables)" << std::endl;
std::cout << " --min-image-cc <num> Per-image CC limit in percent (default: no limit)" << std::endl;
std::cout << " --scaling-iterations <num> Number of scaling iterations with no reference data (default: 3)" << std::endl;
std::cout << " --scaling-output <txt> Output format for scaling results mtz|cif|txt (default: mtz)" << std::endl;
@@ -97,7 +98,8 @@ enum {
OPT_REDO_ROTATION_SPOTS,
OPT_FORCE_ROTATION_LATTICE,
OPT_BANDWIDTH,
OPT_INTEGRATION_RADIUS
OPT_INTEGRATION_RADIUS,
OPT_REJECT_OUTLIERS
};
static option long_options[] = {
@@ -135,6 +137,7 @@ static option long_options[] = {
{"scaling-output", required_argument, nullptr, OPT_SCALING_OUTPUT},
{"bandwidth", required_argument, nullptr, OPT_BANDWIDTH},
{"integration-radius", required_argument, nullptr, OPT_INTEGRATION_RADIUS},
{"reject-outliers", required_argument, nullptr, OPT_REJECT_OUTLIERS},
{nullptr, 0, nullptr, 0}
};
@@ -316,6 +319,7 @@ int main(int argc, char **argv) {
float d_min_spot_finding = 1.5;
std::optional<float> d_min_scale_merge;
std::optional<std::string> integration_radius_arg; // "r1" or "r1,r2,r3"
std::optional<double> outlier_reject_nsigma; // merge per-observation outlier rejection
if (argc == 1) {
print_usage();
@@ -508,6 +512,9 @@ int main(int argc, char **argv) {
case OPT_INTEGRATION_RADIUS:
integration_radius_arg = optarg;
break;
case OPT_REJECT_OUTLIERS:
outlier_reject_nsigma = std::stod(optarg);
break;
case OPT_MIN_IMAGE_CC:
min_image_cc = std::stod(optarg);
break;
@@ -674,6 +681,8 @@ int main(int argc, char **argv) {
scaling_settings.RotationWedgeForScaling(wedge_for_scaling);
scaling_settings.MinPartiality(min_partiality);
scaling_settings.MinCCForImage(min_image_cc);
if (outlier_reject_nsigma)
scaling_settings.OutlierRejectNsigma(*outlier_reject_nsigma);
scaling_settings.FileFormat(intensity_format);
experiment.ImportScalingSettings(scaling_settings);
@@ -1034,15 +1043,11 @@ int main(int argc, char **argv) {
const double isa = (b > 0.0) ? 1.0 / b : std::numeric_limits<double>::infinity();
logger.Info("Error model: sigma'^2 = {:.3f} sigma^2 + ({:.4f} I)^2 ISa = {:.1f}", a, b, isa);
}
if (const char *rj = std::getenv("PR_REJECT")) { // TEMP A/B: per-observation outlier rejection
const double nsig = std::atof(rj);
merge_engine.SetRejectOutliers(nsig);
logger.Info("Outlier rejection enabled: dropping observations > {:.1f} sigma from the per-reflection median", nsig);
}
for (auto &i : indexer.GetIntegrationOutcome())
merge_engine.AddImage(i);
if (merge_engine.RejectedCount() > 0)
logger.Info("Outlier rejection removed {} observations", merge_engine.RejectedCount());
logger.Info("Outlier rejection (>{:.1f} sigma from the per-reflection median) removed {} observations",
experiment.GetScalingSettings().GetOutlierRejectNsigma(), merge_engine.RejectedCount());
auto merged_reflections = merge_engine.ExportReflections();
auto merged_statistics = merge_engine.MergeStats(merged_reflections, indexer.GetIntegrationOutcome(), reference_data);