diff --git a/common/ScalingSettings.cpp b/common/ScalingSettings.cpp index a91aa1ff..866d48a4 100644 --- a/common/ScalingSettings.cpp +++ b/common/ScalingSettings.cpp @@ -115,6 +115,15 @@ ScalingSettings &ScalingSettings::MinCCForImage(double input) { return *this; } +double ScalingSettings::GetOutlierRejectNsigma() const { + return outlier_reject_nsigma; +} + +ScalingSettings &ScalingSettings::OutlierRejectNsigma(double input) { + outlier_reject_nsigma = input; // <= 0 disables; no upper bound (large = effectively off) + return *this; +} + double ScalingSettings::GetMinPartiality() const { return min_partiality; } diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h index f7c724f8..6d548cfc 100644 --- a/common/ScalingSettings.h +++ b/common/ScalingSettings.h @@ -23,6 +23,7 @@ class ScalingSettings { std::optional wedge_for_scaling; double min_partiality = 0.02; double min_cc_for_image = 0.0; + double outlier_reject_nsigma = 6.0; // per-observation merge outlier rejection (XDS/DIALS-style); <= 0 disables double rfree_fraction = 0.05; IntensityFormat intensity_format = IntensityFormat::MTZ; @@ -37,6 +38,7 @@ public: ScalingSettings& HighResolutionLimit_A(double limit); ScalingSettings& MinPartiality(double min_partiality); ScalingSettings& MinCCForImage(double min_cc_for_image); + ScalingSettings& OutlierRejectNsigma(double input); ScalingSettings& RfreeFraction(double input); ScalingSettings& FileFormat(IntensityFormat input); @@ -63,6 +65,7 @@ public: [[nodiscard]] double GetMinPartiality() const; [[nodiscard]] double GetMinCCForImage() const; + [[nodiscard]] double GetOutlierRejectNsigma() const; [[nodiscard]] double GetRfreeFraction() const; [[nodiscard]] IntensityFormat GetFileFormat() const; diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 08f875ba..d604e092 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -31,7 +31,9 @@ MergeOnTheFly::MergeOnTheFly(const DiffractionExperiment &x) ? std::optional(scaling_settings.GetMinCCForImage()) : std::nullopt), min_partiality(scaling_settings.GetMinPartiality()), - generator(scaling_settings.GetMergeFriedel(), space_group_number) { + generator(scaling_settings.GetMergeFriedel(), space_group_number), + reject_outliers(scaling_settings.GetOutlierRejectNsigma() > 0.0), + reject_nsigma(scaling_settings.GetOutlierRejectNsigma()) { } MergeOnTheFly &MergeOnTheFly::ReferenceCell(const std::optional &cell) { @@ -172,8 +174,8 @@ void MergeOnTheFly::RefineErrorModel(const std::vector &outc const double mean = sum_wI / sum_w; error_model_mean_I[key] = static_cast(mean); // Robust centre for outlier rejection: the median intensity (resists the very - // outliers the inverse-variance mean is being protected from). - { + // outliers the inverse-variance mean is being protected from). Only when active. + if (reject_outliers) { std::vector iv; iv.reserve(obs.size()); for (const auto &o: obs) diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index a1ea2776..ab9877f4 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -109,8 +109,7 @@ public: [[nodiscard]] double ErrorModelA() const { return error_model_a; } [[nodiscard]] double ErrorModelB() const { return error_model_b; } - // Enable per-observation outlier rejection (call after RefineErrorModel, before AddImage). - void SetRejectOutliers(double nsigma) { reject_outliers = true; reject_nsigma = nsigma; } + // Outlier rejection (driven by ScalingSettings::GetOutlierRejectNsigma) reports its count. [[nodiscard]] size_t RejectedCount() const { return reject_count; } void AddImage(const IntegrationOutcome& outcome, bool cc_mask = false); diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index ef847d86..93ade850 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -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 High resolution limit for spot finding (default: no limit)" << std::endl; std::cout << " --min-partiality Minimum partiality to accept reflection (default: 0.02)" << std::endl; + std::cout << " --reject-outliers Per-observation merge outlier rejection, N sigma from the per-reflection median (default: 6; 0 disables)" << std::endl; std::cout << " --min-image-cc Per-image CC limit in percent (default: no limit)" << std::endl; std::cout << " --scaling-iterations Number of scaling iterations with no reference data (default: 3)" << std::endl; std::cout << " --scaling-output 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 d_min_scale_merge; std::optional integration_radius_arg; // "r1" or "r1,r2,r3" + std::optional 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::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);