diff --git a/common/ScalingSettings.cpp b/common/ScalingSettings.cpp index c3e75ee5..4cdae567 100644 --- a/common/ScalingSettings.cpp +++ b/common/ScalingSettings.cpp @@ -209,4 +209,37 @@ ScalingSettings &ScalingSettings::ScalingRegularize(bool input) { bool ScalingSettings::GetScalingRegularize() const { return scaling_regularize; -} \ No newline at end of file +} + +ScalingSettings &ScalingSettings::ResolutionCutoff(ResolutionCutoffMethod input) { + resolution_cutoff = input; + return *this; +} + +ResolutionCutoffMethod ScalingSettings::GetResolutionCutoff() const { + return resolution_cutoff; +} + +ScalingSettings &ScalingSettings::ResolutionCCTarget(double input) { + if (input <= 0.0 || input >= 1.0) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Resolution CC target must be between 0 and 1"); + resolution_cc_target = input; + return *this; +} + +double ScalingSettings::GetResolutionCCTarget() const { + return resolution_cc_target; +} + +ScalingSettings &ScalingSettings::ReportShellCount(int input) { + if (input < 1) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Number of report shells must be at least 1"); + report_shell_count = input; + return *this; +} + +int ScalingSettings::GetReportShellCount() const { + return report_shell_count; +} diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h index 3a4256d3..8545d886 100644 --- a/common/ScalingSettings.h +++ b/common/ScalingSettings.h @@ -8,6 +8,11 @@ enum class IntensityFormat { Text, mmCIF, MTZ}; +// How the high-resolution cutoff for the written reflections and the reported shell table is chosen +// when no explicit --scaling-high-resolution is given. Off = keep the full (detector-edge) range; +// CCHalfLogistic = fit the CC1/2 fall-off and cut one shell past cc_target (DIALS-style, generous). +enum class ResolutionCutoffMethod { Off, CCHalfLogistic }; + class ScalingSettings { bool refine_b = false; double max_b = 200.0; @@ -45,6 +50,13 @@ class ScalingSettings { double rfree_fraction = 0.05; IntensityFormat intensity_format = IntensityFormat::mmCIF; + // Automatic high-resolution cutoff for the written reflections + reported shells (not the scaling + // or the error model, and not the per-image _process.h5). Applied only when no explicit + // high_resolution_limit_A is set - that manual limit always wins and disables the auto-cut. + ResolutionCutoffMethod resolution_cutoff = ResolutionCutoffMethod::CCHalfLogistic; + double resolution_cc_target = 0.30; // CC1/2 value defining the fall-off limit before the +1 shell + int report_shell_count = 10; // number of resolution shells in the reported statistics table + bool scaling_regularize = false; public: ScalingSettings& RefineB(bool input); @@ -66,6 +78,10 @@ public: ScalingSettings& FileFormat(IntensityFormat input); ScalingSettings& ScalingRegularize(bool input); + ScalingSettings& ResolutionCutoff(ResolutionCutoffMethod input); + ScalingSettings& ResolutionCCTarget(double input); + ScalingSettings& ReportShellCount(int input); + [[nodiscard]] bool GetRefineB() const; [[nodiscard]] bool GetRefineWedge() const; @@ -96,4 +112,8 @@ public: [[nodiscard]] double GetRfreeFraction() const; [[nodiscard]] IntensityFormat GetFileFormat() const; [[nodiscard]] bool GetScalingRegularize() const; + + [[nodiscard]] ResolutionCutoffMethod GetResolutionCutoff() const; + [[nodiscard]] double GetResolutionCCTarget() const; + [[nodiscard]] int GetReportShellCount() const; }; diff --git a/image_analysis/scale_merge/CMakeLists.txt b/image_analysis/scale_merge/CMakeLists.txt index fe81bb17..dcac4989 100644 --- a/image_analysis/scale_merge/CMakeLists.txt +++ b/image_analysis/scale_merge/CMakeLists.txt @@ -9,6 +9,8 @@ ADD_LIBRARY(JFJochScaleMerge ScaleOnTheFly.h RotationScaleMerge.cpp RotationScaleMerge.h + ResolutionCutoff.cpp + ResolutionCutoff.h HKLKey.cpp HKLKey.h ScalingResult.h diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index c924b6f1..9228fb67 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -572,11 +572,13 @@ void CalcPossibleReflections(int space_group_number , MergeStatistics MergeOnTheFly::MergeStats(const std::vector &merged, const std::vector &integration_outcome, - const std::vector &reference) { + const std::vector &reference, + std::optional d_min_override) { - constexpr int n_shells = 10; + const int n_shells = scaling_settings.GetReportShellCount(); - auto d_min_limit_A = scaling_settings.GetHighResolutionLimit_A(); + auto d_min_limit_A = d_min_override.has_value() + ? d_min_override : scaling_settings.GetHighResolutionLimit_A(); std::unordered_map reference_intensities; if (!reference.empty()) { diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 0ffe8fc7..9c0ee572 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -153,9 +153,13 @@ public: [[nodiscard]] std::vector DeltaCChalfReject(const std::vector &outcomes, double nsigma) const; + // d_min_override, when set, is the effective high-resolution limit for the shell table (used for + // the automatic resolution cutoff computed by the caller); otherwise the manual + // ScalingSettings high-resolution limit stands. The number of shells is ScalingSettings::ReportShellCount. MergeStatistics MergeStats(const std::vector &merged, const std::vector &reflections, - const std::vector &reference = {}); + const std::vector &reference = {}, + std::optional d_min_override = std::nullopt); std::vector ExportReflections(); }; diff --git a/image_analysis/scale_merge/ResolutionCutoff.cpp b/image_analysis/scale_merge/ResolutionCutoff.cpp new file mode 100644 index 00000000..65df0f44 --- /dev/null +++ b/image_analysis/scale_merge/ResolutionCutoff.cpp @@ -0,0 +1,174 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "ResolutionCutoff.h" + +#include +#include +#include + +#include "../../common/CorrelationCoefficient.h" + +namespace { + // Fine CC1/2 bins for the fit (finer than the 10 reported shells, per the design). A bin needs + // this many merged reflections for its CC1/2 to be trusted. + constexpr int N_FIT_BINS = 25; + constexpr int MIN_BIN_COUNT = 10; + constexpr int MIN_FIT_BINS = 5; // need at least this many usable bins for a fit + constexpr int EXTEND_BINS_PAST_FALLOFF = 2; // bins kept beyond the first sub-target bin + constexpr double SHELLS_FOR_EXTENSION = 10.0; // "+1 shell" = one 10-shell width in s (report-independent) + + double Logistic(double s, double k, double s0) { + return 1.0 / (1.0 + std::exp(k * (s - s0))); + } + + // Weighted (equal-weight) sum of squared residuals of the logistic against the binned CC1/2. + double FitSSE(const std::vector &s, const std::vector &cc, double k, double s0) { + double sse = 0.0; + for (size_t i = 0; i < s.size(); ++i) { + const double r = cc[i] - Logistic(s[i], k, s0); + sse += r * r; + } + return sse; + } +} + +ResolutionCutoffResult ComputeCCHalfLogisticCutoff(const std::vector &merged, + double cc_target, Logger &logger) { + ResolutionCutoffResult result; + if (!(cc_target > 0.0 && cc_target < 1.0)) { + result.note = "invalid CC target"; + return result; + } + + // s = 1/d^2 range over the merged reflections that carry a half-set pair. + double s_lo = std::numeric_limits::max(), s_hi = 0.0; + double d_data_min = std::numeric_limits::max(); + for (const auto &m : merged) { + if (!(m.d > 0.0f) || !std::isfinite(m.I_half[0]) || !std::isfinite(m.I_half[1])) + continue; + const double s = 1.0 / (static_cast(m.d) * m.d); + s_lo = std::min(s_lo, s); + s_hi = std::max(s_hi, s); + d_data_min = std::min(d_data_min, static_cast(m.d)); + } + if (!(s_lo < s_hi)) { + result.note = "no half-set data for CC1/2 fit"; + return result; + } + + // Bin CC1/2 against s (equal width in s, matching the reporting shells which are equal in 1/d^2). + const double bin_w = (s_hi - s_lo) / N_FIT_BINS; + std::vector bin_cc(N_FIT_BINS); + std::vector bin_n(N_FIT_BINS, 0); + for (const auto &m : merged) { + if (!(m.d > 0.0f) || !std::isfinite(m.I_half[0]) || !std::isfinite(m.I_half[1])) + continue; + const double s = 1.0 / (static_cast(m.d) * m.d); + int b = static_cast((s - s_lo) / bin_w); + b = std::clamp(b, 0, N_FIT_BINS - 1); + bin_cc[b].Add(m.I_half[0], m.I_half[1]); + ++bin_n[b]; + } + + // Usable bins (enough counts), in ascending-s order, with their bin-centre s. + std::vector s_bin, cc_bin; + for (int b = 0; b < N_FIT_BINS; ++b) { + if (bin_n[b] < MIN_BIN_COUNT) continue; + const double cc = bin_cc[b].GetCC(); + if (!std::isfinite(cc)) continue; + s_bin.push_back(s_lo + (b + 0.5) * bin_w); + cc_bin.push_back(cc); + } + if (static_cast(s_bin.size()) < MIN_FIT_BINS) { + result.note = "too few usable CC1/2 bins"; + return result; + } + + // Restrict to the contiguous fall-off from low res: keep bins up to a couple past the first one + // that drops below cc_target, so a high-res noise blip cannot pull the fit back up. If the lowest + // bin is already below cc_target there is no low-res plateau to anchor on - bail out. + if (cc_bin.front() < cc_target) { + result.note = "no low-resolution CC1/2 plateau"; + return result; + } + size_t keep = cc_bin.size(); + for (size_t i = 0; i < cc_bin.size(); ++i) { + if (cc_bin[i] < cc_target) { + keep = std::min(cc_bin.size(), i + 1 + EXTEND_BINS_PAST_FALLOFF); + break; + } + } + s_bin.resize(keep); + cc_bin.resize(keep); + if (static_cast(s_bin.size()) < MIN_FIT_BINS) { + result.note = "too few CC1/2 bins in the fall-off region"; + return result; + } + + // Fit the logistic by a grid search over (k>0, s0) then a local coordinate-descent refine + // (dependency-free; the fall-off is smooth and the grid lands close). s0 spans the s range; k + // spans transitions from very gradual to very sharp relative to that range. + const double s_range = s_hi - s_lo; + double best_k = 0.0, best_s0 = 0.0, best_sse = std::numeric_limits::max(); + constexpr int N_S0 = 60, N_K = 40; + const double k_min = 2.0 / s_range, k_max = 200.0 / s_range; + for (int ik = 0; ik < N_K; ++ik) { + const double k = k_min * std::pow(k_max / k_min, static_cast(ik) / (N_K - 1)); + for (int is = 0; is < N_S0; ++is) { + const double s0 = s_lo + s_range * static_cast(is) / (N_S0 - 1); + const double sse = FitSSE(s_bin, cc_bin, k, s0); + if (sse < best_sse) { best_sse = sse; best_k = k; best_s0 = s0; } + } + } + + double k = best_k, s0 = best_s0; + double step_s0 = s_range / N_S0, step_k = best_k * 0.5; + for (int iter = 0; iter < 200; ++iter) { + bool improved = false; + for (const double ds : {step_s0, -step_s0}) { + const double sse = FitSSE(s_bin, cc_bin, k, s0 + ds); + if (sse < best_sse) { best_sse = sse; s0 += ds; improved = true; } + } + for (const double dk : {step_k, -step_k}) { + const double kt = k + dk; + if (kt <= 0.0) continue; + const double sse = FitSSE(s_bin, cc_bin, kt, s0); + if (sse < best_sse) { best_sse = sse; k = kt; improved = true; } + } + if (!improved) { step_s0 *= 0.5; step_k *= 0.5; } + if (step_s0 < 1e-6 * s_range && step_k < 1e-6 * best_k) break; + } + + // s where the fitted CC1/2 crosses cc_target, then "one shell too far". + const double s_cross = s0 + std::log(1.0 / cc_target - 1.0) / k; + const double delta_s = s_range / SHELLS_FOR_EXTENSION; + const double s_final = s_cross + delta_s; + + // No cut if the fall-off is beyond the measured edge (CC1/2 still healthy at the highest s). + if (s_final >= s_hi) { + result.note = "CC1/2 does not fall off within the measured range"; + return result; + } + // Low-resolution floor: never cut into good low-res data. A fit that puts the cutoff within two + // shells of the lowest-res data is not a real fall-off - keep the full range and warn. + if (s_final <= s_lo + 2.0 * delta_s) { + logger.Warning("Resolution cutoff fit landed at low resolution (degenerate CC1/2 fall-off); " + "keeping the full resolution range"); + result.note = "degenerate low-resolution fit"; + return result; + } + + double d_cut = 1.0 / std::sqrt(s_final); + d_cut = std::max(d_cut, d_data_min); // cannot cut beyond the highest-resolution reflection + + // A cut that is not meaningfully coarser than the data edge is a no-op. + if (d_cut <= d_data_min * 1.001) { + result.note = "CC1/2 healthy to the detector edge"; + return result; + } + + result.d_cut = d_cut; + result.note = "CC1/2 logistic fall-off, +1 shell"; + return result; +} diff --git a/image_analysis/scale_merge/ResolutionCutoff.h b/image_analysis/scale_merge/ResolutionCutoff.h new file mode 100644 index 00000000..ec5b9a16 --- /dev/null +++ b/image_analysis/scale_merge/ResolutionCutoff.h @@ -0,0 +1,29 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include + +#include "../../common/Logger.h" +#include "../../common/Reflection.h" // MergedReflection + +// Automatic high-resolution cutoff from the CC1/2 fall-off of the merged half-sets (DIALS-style). +// The merge itself, the error model and the per-image _process.h5 are left untouched - only the +// written reflections and the reported shell table should be trimmed to the returned d_cut. +// +// Method (see docs/rugnux_resolution_cutoff_design.md): bin CC1/2 against s = 1/d^2 in fine bins, +// fit a logistic CC1/2(s) = 1/(1+exp(k*(s-s0))) to the contiguous-from-low-res fall-off, take the s +// where the fit crosses cc_target, then extend by one mean (10-shell) shell width in s ("one shell +// too far", generous). d_cut is nullopt when the fit is degenerate (too few bins, flat/non-monotone) +// or CC1/2 never falls below cc_target inside the measured range - the caller then keeps the full +// range. cc_target in (0,1); merged must carry finite I_half[0]/I_half[1] to contribute. +struct ResolutionCutoffResult { + std::optional d_cut; // high-resolution limit (A); nullopt => keep the full range + std::string note; // human-readable description of the decision (for logging) +}; + +ResolutionCutoffResult ComputeCCHalfLogisticCutoff(const std::vector &merged, + double cc_target, Logger &logger); diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index f38b90e4..13f07362 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -16,6 +16,7 @@ #include #include "HKLKey.h" +#include "ResolutionCutoff.h" #include "../../common/CorrelationCoefficient.h" #include "../../common/CrystalLattice.h" #include "../../common/Definitions.h" @@ -159,6 +160,9 @@ RotationScaleMerge::RotationScaleMerge(const DiffractionExperiment &experiment, rfree_fraction = s.GetRfreeFraction(); scale_fulls = s.GetScaleFulls(); scaling_iter = std::max(1, scaling_iterations); + resolution_cutoff_method = s.GetResolutionCutoff(); + resolution_cc_target = s.GetResolutionCCTarget(); + report_shell_count = s.GetReportShellCount(); if (const auto forced = s.GetForcedMosaicity(); forced.has_value()) mosaicity_deg = *forced; else @@ -1033,6 +1037,34 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool result.merged.push_back(mr); } + // Automatic high-resolution cutoff (post-merge): trim the written reflections + reported shells to + // the CC1/2 fall-off. The scaling, combine and error model above already ran over the full range, + // and the per-image _process.h5 is written elsewhere from the partials, so no data is lost. A manual + // --scaling-high-resolution (d_min_limit) wins; the P1 search merge (for_search) is never cut, so + // the space-group search still sees the full range. + std::optional effective_d_min = d_min_limit; + if (!effective_d_min && !for_search && resolution_cutoff_method == ResolutionCutoffMethod::CCHalfLogistic) { + const auto rc = ComputeCCHalfLogisticCutoff(result.merged, resolution_cc_target, logger); + if (rc.d_cut) { + effective_d_min = rc.d_cut; + logger.Info("Auto resolution cutoff: {:.2f} A ({}; override with --scaling-high-resolution)", + *rc.d_cut, rc.note); + } else { + logger.Info("Auto resolution cutoff: none ({}); keeping the full resolution range", rc.note); + } + } + if (effective_d_min) { + result.merged.erase(std::remove_if(result.merged.begin(), result.merged.end(), + [&](const MergedReflection &m) { return std::isfinite(m.d) && m.d < *effective_d_min; }), + result.merged.end()); + // Recompute the merged resolution span for the R-free binning below over the trimmed set. + d_min = std::numeric_limits::max(); d_max = 0.0f; + for (const auto &m : result.merged) { + if (!std::isfinite(m.d) || m.d <= 0.0f) continue; + d_min = std::min(d_min, m.d); d_max = std::max(d_max, m.d); + } + } + if (rfree_fraction > 0.0 && !result.merged.empty() && d_min < d_max && d_min > 0.0f) { constexpr int n_shells = 20; ResolutionShells shells(d_min * 0.999f, d_max * 1.001f, n_shells); @@ -1049,12 +1081,12 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool if (reject_count > 0) logger.Info("Merge outlier rejection: dropped {} observations", reject_count); - // ---- Statistics (10 shells): completeness, multiplicity, , R_meas, CC1/2. ---- - constexpr int n_shells = 10; + // ---- Statistics (report_shell_count shells): completeness, multiplicity, , R_meas, CC1/2. ---- + const int n_shells = report_shell_count; float sd_min = std::numeric_limits::max(), sd_max = 0.0f; for (const auto &m : result.merged) { if (!std::isfinite(m.d) || m.d <= 0.0f) continue; - if (d_min_limit && m.d < *d_min_limit) continue; + if (effective_d_min && m.d < *effective_d_min) continue; sd_min = std::min(sd_min, m.d); sd_max = std::max(sd_max, m.d); } if (!(sd_min < sd_max && sd_min > 0.0f)) diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index f5ab0ae7..f40b046c 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -98,6 +98,11 @@ private: bool scale_fulls = true; double mosaicity_deg = 0.1; float ice_half_width_q = 0.0f; + // Automatic high-resolution cutoff for the written reflections + reported shells (post-merge; the + // scaling, combine and error model always run over the full range). Manual d_min_limit wins. + ResolutionCutoffMethod resolution_cutoff_method = ResolutionCutoffMethod::Off; + double resolution_cc_target = 0.30; + int report_shell_count = 10; // Flat buffers, allocated once by Ingest() and reused across Run() calls. std::vector partials; // all per-frame partials, grouped by frame diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 716c7711..80b5dd9c 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -28,6 +28,7 @@ #include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h" #include "../image_analysis/scale_merge/Merge.h" #include "../image_analysis/scale_merge/RotationScaleMerge.h" +#include "../image_analysis/scale_merge/ResolutionCutoff.h" #include "../image_analysis/scale_merge/ScalingResult.h" #include "../image_analysis/scale_merge/SearchSpaceGroup.h" #include "../image_analysis/scale_merge/TwinningAnalysis.h" @@ -499,8 +500,32 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { ScaleMergeResult out; out.merged = merge_engine.ExportReflections(); + + // Automatic high-resolution cutoff (post-merge): trim the written reflections + reported + // shells to the CC1/2 fall-off. The merge, scaling and error model above ran over the full + // range, and the _process.h5 is written from the per-image outcomes, so no data is lost. A + // manual --scaling-high-resolution wins; the P1 search merge (for_search) is never cut so the + // space-group search still sees the full range. (Rotation is cut inside RotationScaleMerge.) + const auto &cut_ss = experiment_.GetScalingSettings(); + std::optional effective_d_min = cut_ss.GetHighResolutionLimit_A(); + if (!effective_d_min && !for_search + && cut_ss.GetResolutionCutoff() == ResolutionCutoffMethod::CCHalfLogistic) { + const auto rc = ComputeCCHalfLogisticCutoff(out.merged, cut_ss.GetResolutionCCTarget(), logger); + if (rc.d_cut) { + effective_d_min = rc.d_cut; + logger.Info("Auto resolution cutoff: {:.2f} A ({}; override with --scaling-high-resolution)", + *rc.d_cut, rc.note); + } else { + logger.Info("Auto resolution cutoff: none ({}); keeping the full resolution range", rc.note); + } + } + if (effective_d_min) + out.merged.erase(std::remove_if(out.merged.begin(), out.merged.end(), + [&](const MergedReflection &m) { return std::isfinite(m.d) && m.d < *effective_d_min; }), + out.merged.end()); + phase("Computing statistics"); - out.statistics = merge_engine.MergeStats(out.merged, merge_input, config_.reference_data); + out.statistics = merge_engine.MergeStats(out.merged, merge_input, config_.reference_data, effective_d_min); result.error_model_isa = merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0; logger.Info("Merge complete ({} unique reflections, {})", out.merged.size(), label); return out; diff --git a/tools/rugnux_cli.cpp b/tools/rugnux_cli.cpp index 02ac0c6d..bb3563ad 100644 --- a/tools/rugnux_cli.cpp +++ b/tools/rugnux_cli.cpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include @@ -24,6 +25,7 @@ #include "../image_analysis/scale_merge/Merge.h" #include "../image_analysis/scale_merge/ScaleOnTheFly.h" #include "../image_analysis/scale_merge/RotationScaleMerge.h" +#include "../image_analysis/scale_merge/ResolutionCutoff.h" #include "../image_analysis/scale_merge/TwinningAnalysis.h" #include "../rugnux/Rugnux.h" @@ -79,7 +81,10 @@ void print_usage() { std::cout << " --smooth-g[=deg] rot3d: smooth per-frame scale G over a deg-degree rotation range (XDS DELPHI-like) before the combine (default: 5 for rot3d; 0 = off)" << 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" << std::endl; - std::cout << " --scaling-high-resolution High resolution limit for scaling/merging (default: no limit)" << std::endl; + std::cout << " --scaling-high-resolution High resolution limit for scaling/merging (manual override; default: no limit)" << std::endl; + std::cout << " --resolution-cutoff Automatic high-resolution cutoff for the written reflections + reported shells: cc-logistic|off (default: cc-logistic; ignored when --scaling-high-resolution is set)" << std::endl; + std::cout << " --resolution-cc-target CC1/2 target defining the cc-logistic fall-off (default: 0.30)" << std::endl; + 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 << " --mosaicity Diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed" << std::endl; @@ -123,6 +128,9 @@ enum { OPT_MIN_IMAGE_CC, OPT_SCALING_ITERATIONS, OPT_SCALING_HIGH_RESOLUTION, + OPT_RESOLUTION_CUTOFF, + OPT_RESOLUTION_CC_TARGET, + OPT_RESOLUTION_SHELLS, OPT_SCALING_OUTPUT, OPT_SINGLE_PASS_ROTATION, OPT_REDO_ROTATION_SPOTS, @@ -213,6 +221,9 @@ static option long_options[] = { {"min-image-cc", required_argument, nullptr, OPT_MIN_IMAGE_CC}, {"scaling-iterations", required_argument, nullptr, OPT_SCALING_ITERATIONS}, {"scaling-high-resolution", required_argument, nullptr, OPT_SCALING_HIGH_RESOLUTION}, + {"resolution-cutoff", required_argument, nullptr, OPT_RESOLUTION_CUTOFF}, + {"resolution-cc-target", required_argument, nullptr, OPT_RESOLUTION_CC_TARGET}, + {"resolution-shells", required_argument, nullptr, OPT_RESOLUTION_SHELLS}, {"scaling-output", required_argument, nullptr, OPT_SCALING_OUTPUT}, {"bandwidth", required_argument, nullptr, OPT_BANDWIDTH}, {"integration-radius", required_argument, nullptr, OPT_INTEGRATION_RADIUS}, @@ -448,6 +459,9 @@ int main(int argc, char **argv) { float d_min_spot_finding = 1.5; std::optional d_min_scale_merge; + std::optional resolution_cutoff_method; // --resolution-cutoff cc-logistic|off + std::optional resolution_cc_target; // --resolution-cc-target + std::optional report_shell_count; // --resolution-shells std::optional integration_radius_arg; // "r1" or "r1,r2,r3" std::optional integrator_mode; // --integrator boxsum|gaussian|empirical std::optional outlier_reject_nsigma; // merge per-observation outlier rejection @@ -683,6 +697,27 @@ int main(int argc, char **argv) { case OPT_SCALING_HIGH_RESOLUTION: d_min_scale_merge = atof(optarg); break; + case OPT_RESOLUTION_CUTOFF: + if (strcmp(optarg, "cc-logistic") == 0) + resolution_cutoff_method = ResolutionCutoffMethod::CCHalfLogistic; + else if (strcmp(optarg, "off") == 0) + resolution_cutoff_method = ResolutionCutoffMethod::Off; + else { + logger.Error("Invalid --resolution-cutoff value: {} (expected cc-logistic|off)", optarg); + print_usage(); + exit(EXIT_FAILURE); + } + break; + case OPT_RESOLUTION_CC_TARGET: + resolution_cc_target = parse_double_arg(optarg, "--resolution-cc-target", logger); + break; + case OPT_RESOLUTION_SHELLS: + report_shell_count = atoi(optarg); + if (report_shell_count.value() < 1) { + logger.Error("Invalid --resolution-shells value: {} (must be >= 1)", report_shell_count.value()); + exit(EXIT_FAILURE); + } + break; case OPT_SCALING_OUTPUT: if (strcmp(optarg, "mtz") == 0) { intensity_format = IntensityFormat::MTZ; @@ -882,6 +917,9 @@ int main(int argc, char **argv) { ScalingSettings scaling_settings; if (d_min_scale_merge) scaling_settings.HighResolutionLimit_A(d_min_scale_merge.value()); + if (resolution_cutoff_method) scaling_settings.ResolutionCutoff(*resolution_cutoff_method); + if (resolution_cc_target) scaling_settings.ResolutionCCTarget(*resolution_cc_target); + if (report_shell_count) scaling_settings.ReportShellCount(*report_shell_count); scaling_settings.MergeFriedel(!anomalous_mode); scaling_settings.RefineB(refine_bfactor); scaling_settings.MinPartiality(min_partiality); @@ -933,7 +971,29 @@ int main(int argc, char **argv) { for (size_t i = 0; i < reflections.size(); ++i) merge_engine.AddImage(reflections[i], static_cast(i)); merged_reflections = merge_engine.ExportReflections(); - merged_statistics = merge_engine.MergeStats(merged_reflections, reflections, reference_data); + + // Automatic high-resolution cutoff (post-merge), matching the full-analysis path: a manual + // --scaling-high-resolution wins, otherwise trim the written reflections + reported shells + // to the CC1/2 fall-off. (Rotation is cut inside RotationScaleMerge above.) + const auto &cut_ss = experiment.GetScalingSettings(); + std::optional effective_d_min = cut_ss.GetHighResolutionLimit_A(); + if (!effective_d_min + && cut_ss.GetResolutionCutoff() == ResolutionCutoffMethod::CCHalfLogistic) { + const auto rc = ComputeCCHalfLogisticCutoff(merged_reflections, cut_ss.GetResolutionCCTarget(), + logger); + if (rc.d_cut) { + effective_d_min = rc.d_cut; + logger.Info("Auto resolution cutoff: {:.2f} A ({}; override with --scaling-high-resolution)", + *rc.d_cut, rc.note); + } + } + if (effective_d_min) + merged_reflections.erase(std::remove_if(merged_reflections.begin(), merged_reflections.end(), + [&](const MergedReflection &m) { return std::isfinite(m.d) && m.d < *effective_d_min; }), + merged_reflections.end()); + + merged_statistics = merge_engine.MergeStats(merged_reflections, reflections, reference_data, + effective_d_min); error_model_isa = merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0; } @@ -1123,6 +1183,9 @@ int main(int argc, char **argv) { scaling_settings.SmoothGDegrees(smooth_g_deg_arg.value_or(rotation_indexing ? SMOOTH_G_DEFAULT_DEG : 0.0)); if (d_min_scale_merge) scaling_settings.HighResolutionLimit_A(d_min_scale_merge.value()); + if (resolution_cutoff_method) scaling_settings.ResolutionCutoff(*resolution_cutoff_method); + if (resolution_cc_target) scaling_settings.ResolutionCCTarget(*resolution_cc_target); + if (report_shell_count) scaling_settings.ReportShellCount(*report_shell_count); scaling_settings.MergeFriedel(!anomalous_mode); scaling_settings.RefineB(refine_bfactor); scaling_settings.MinPartiality(min_partiality);