diff --git a/common/ScalingSettings.cpp b/common/ScalingSettings.cpp index f0cd7ac5..1044cdd8 100644 --- a/common/ScalingSettings.cpp +++ b/common/ScalingSettings.cpp @@ -133,6 +133,15 @@ bool ScalingSettings::GetCombine3D() const { return combine_3d; } +ScalingSettings &ScalingSettings::ScaleFulls(bool input) { + scale_fulls = input; + return *this; +} + +bool ScalingSettings::GetScaleFulls() const { + return scale_fulls; +} + double ScalingSettings::GetMinPartiality() const { return min_partiality; } diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h index cef40bbb..39e6cb0d 100644 --- a/common/ScalingSettings.h +++ b/common/ScalingSettings.h @@ -30,6 +30,10 @@ class ScalingSettings { // model, which stays Rotation - see Combine3D). bool combine_3d = false; + // Scale fulls: after the rot3d combine, refit a per-frame scale on the combined fulls (XDS + // order). A no-op without rot3d. + bool scale_fulls = false; + double rfree_fraction = 0.05; IntensityFormat intensity_format = IntensityFormat::MTZ; @@ -45,6 +49,7 @@ public: ScalingSettings& MinCCForImage(double min_cc_for_image); ScalingSettings& OutlierRejectNsigma(double input); ScalingSettings& Combine3D(bool input); + ScalingSettings& ScaleFulls(bool input); ScalingSettings& RfreeFraction(double input); ScalingSettings& FileFormat(IntensityFormat input); @@ -73,6 +78,7 @@ public: [[nodiscard]] double GetMinCCForImage() const; [[nodiscard]] double GetOutlierRejectNsigma() const; [[nodiscard]] bool GetCombine3D() const; + [[nodiscard]] bool GetScaleFulls() const; [[nodiscard]] double GetRfreeFraction() const; [[nodiscard]] IntensityFormat GetFileFormat() const; diff --git a/image_analysis/scale_merge/CMakeLists.txt b/image_analysis/scale_merge/CMakeLists.txt index 6681329d..190040c7 100644 --- a/image_analysis/scale_merge/CMakeLists.txt +++ b/image_analysis/scale_merge/CMakeLists.txt @@ -5,8 +5,6 @@ ADD_LIBRARY(JFJochScaleMerge Merge.h ScaleOnTheFly.cpp ScaleOnTheFly.h - GlobalScale.cpp - GlobalScale.h Combine3D.cpp Combine3D.h HKLKey.cpp diff --git a/image_analysis/scale_merge/GlobalScale.cpp b/image_analysis/scale_merge/GlobalScale.cpp deleted file mode 100644 index c01c322b..00000000 --- a/image_analysis/scale_merge/GlobalScale.cpp +++ /dev/null @@ -1,489 +0,0 @@ -// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#include "GlobalScale.h" - -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -namespace { - struct HKLKey { - int64_t h = 0; - int64_t k = 0; - int64_t l = 0; - bool is_positive = true; // only relevant if opt.merge_friedel == false - - bool operator==(const HKLKey &o) const noexcept { - return h == o.h && k == o.k && l == o.l && is_positive == o.is_positive; - } - }; - - struct HKLKeyHash { - size_t operator()(const HKLKey &key) const noexcept { - auto mix = [](uint64_t x) { - x ^= x >> 33; - x *= 0xff51afd7ed558ccdULL; - x ^= x >> 33; - x *= 0xc4ceb9fe1a85ec53ULL; - x ^= x >> 33; - return x; - }; - const uint64_t a = static_cast(key.h); - const uint64_t b = static_cast(key.k); - const uint64_t c = static_cast(key.l); - const uint64_t d = static_cast(key.is_positive ? 1 : 0); - return static_cast(mix(a) ^ (mix(b) << 1) ^ (mix(c) << 2) ^ (mix(d) << 3)); - } - }; - - inline double SafeSigma(double s, double min_sigma) { - if (!std::isfinite(s) || s <= 0.0) - return min_sigma; - return std::max(s, min_sigma); - } - - inline double SafeD(double d) { - if (!std::isfinite(d) || d <= 0.0) - return std::numeric_limits::quiet_NaN(); - return d; - } - - inline int SafeToInt(int64_t x) { - if (x < std::numeric_limits::min() || x > std::numeric_limits::max()) - throw std::out_of_range("HKL index out of int range for Gemmi"); - return static_cast(x); - } - - inline double SafeInv(double x, double fallback) { - if (!std::isfinite(x) || x == 0.0) - return fallback; - return 1.0 / x; - } - - inline HKLKey CanonicalizeHKLKey(const Reflection &r, const GlobalScaleOptions &opt) { - HKLKey key{}; - key.h = r.h; - key.k = r.k; - key.l = r.l; - key.is_positive = true; - - if (!opt.space_group.has_value()) { - if (!opt.merge_friedel) { - const HKLKey neg{-r.h, -r.k, -r.l, true}; - const bool pos = std::tie(key.h, key.k, key.l) >= std::tie(neg.h, neg.k, neg.l); - if (!pos) { - key.h = -key.h; - key.k = -key.k; - key.l = -key.l; - key.is_positive = false; - } - } - return key; - } - - const gemmi::SpaceGroup &sg = *opt.space_group; - const gemmi::GroupOps gops = sg.operations(); - const gemmi::ReciprocalAsu rasu(&sg); - - const gemmi::Op::Miller in{{SafeToInt(r.h), SafeToInt(r.k), SafeToInt(r.l)}}; - const auto [asu_hkl, sign_plus] = rasu.to_asu_sign(in, gops); - - key.h = asu_hkl[0]; - key.k = asu_hkl[1]; - key.l = asu_hkl[2]; - key.is_positive = opt.merge_friedel ? true : sign_plus; - return key; - } - - struct IntensityRotResidual { - IntensityRotResidual(const Reflection &r, double sigma_obs, double wedge_deg) - : Iobs_(static_cast(r.I)), - weight_(SafeInv(sigma_obs, 1.0)), - delta_phi_(r.delta_phi_deg), - lp_(SafeInv(r.rlp, 1.0)), - c1_(r.zeta / std::sqrt(2.0)) { - } - - template - bool operator()(const T *const G, - const T *const mosaicity, - const T *const Itrue, - const T *const wedge, - T *residual) const { - T partiality; - if (mosaicity[0] >= 0.0) { - const T half_wedge = wedge[0] / T(2.0); - const T arg_plus = T(delta_phi_ + half_wedge) * T(c1_) / mosaicity[0]; - const T arg_minus = T(delta_phi_ - half_wedge) * T(c1_) / mosaicity[0]; - partiality = (ceres::erf(arg_plus) - ceres::erf(arg_minus)) / T(2.0); - } else - partiality = T(1.0); - - const T Ipred = G[0] * partiality * T(lp_) * Itrue[0]; - residual[0] = (Ipred - T(Iobs_)) * T(weight_); - return true; - } - - double Iobs_; - double weight_; - double delta_phi_; - double lp_; - double c1_; - }; - - struct IntensityStillResidual { - IntensityStillResidual(const Reflection &r, double sigma_obs) - : Iobs_(static_cast(r.I)), - weight_(SafeInv(sigma_obs, 1.0)), - lp_(SafeInv(r.rlp, 1.0)), - dist_ewald_sq_(r.dist_ewald * r.dist_ewald) { - } - - template - bool operator()(const T *const G, - const T *const R, - const T *const Itrue, - T *residual) const { - const T partiality = ceres::exp(-T(dist_ewald_sq_)/R[0]); - const T Ipred = G[0] * partiality * T(lp_) * Itrue[0]; - residual[0] = (Ipred - T(Iobs_)) * T(weight_); - return true; - } - - double Iobs_; - double weight_; - double lp_; - double dist_ewald_sq_; - }; - - struct IntensityFixedResidual { - IntensityFixedResidual(const Reflection &r, double sigma_obs, double partiality) - : Iobs_(static_cast(r.I)), - weight_(SafeInv(sigma_obs, 1.0)), - corr_(partiality * SafeInv(r.rlp, 1.0)) - {} - - template - bool operator()(const T *const G, - const T *const Itrue, - T *residual) const { - const T Ipred = T(corr_) * G[0] * Itrue[0]; - residual[0] = (Ipred - T(Iobs_)) * T(weight_); - return true; - } - - double Iobs_; - double weight_; - double corr_; - }; - - struct ScaleRegularizationResidual { - explicit ScaleRegularizationResidual(double sigma_k) - : inv_sigma_(SafeInv(sigma_k, 1.0)) { - } - - template - bool operator()(const T *const k, T *residual) const { - residual[0] = (k[0] - T(1.0)) * T(inv_sigma_); - return true; - } - - double inv_sigma_; - }; - - struct SmoothnessRegularizationResidual { - explicit SmoothnessRegularizationResidual(double sigma) - : inv_sigma_(SafeInv(sigma, 1.0)) { - } - - template - bool operator()(const T *const k0, - const T *const k1, - const T *const k2, - T *residual) const { - residual[0] = (ceres::log(k0[0]) + ceres::log(k2[0]) - T(2.0) * ceres::log(k1[0])) * T(inv_sigma_); - return true; - } - - double inv_sigma_; - }; - - struct ObsRef { - const Reflection *r = nullptr; - int img_id = 0; - int hkl_slot = -1; - double sigma = 0.0; - }; - - void scale(const GlobalScaleOptions &opt, - std::vector &g, - std::vector &mosaicity, - std::vector &R_sq, - const std::vector &image_slot_used, - bool rotation_crystallography, - size_t nhkl, - const std::vector &obs) { - ceres::Problem problem; - - std::vector Itrue(nhkl, 0.0); - - // Initialize Itrue from per-HKL median of observed intensities - { - std::vector > per_hkl_I(nhkl); - for (const auto &o: obs) { - per_hkl_I[o.hkl_slot].push_back(static_cast(o.r->I)); - } - for (int h = 0; h < nhkl; ++h) { - auto &v = per_hkl_I[h]; - if (v.empty()) { - Itrue[h] = std::max(opt.min_sigma, 1e-6); - continue; - } - std::nth_element(v.begin(), v.begin() + static_cast(v.size() / 2), v.end()); - double med = v[v.size() / 2]; - if (!std::isfinite(med) || med <= opt.min_sigma) - med = opt.min_sigma; - Itrue[h] = med; - } - } - - double wedge = opt.wedge_deg.value_or(0.0); - - for (const auto &o: obs) { - switch (opt.partiality_model) { - case GlobalScaleOptions::PartialityModel::Rotation: { - auto *cost = new ceres::AutoDiffCostFunction( - new IntensityRotResidual(*o.r, o.sigma, opt.wedge_deg.value_or(0.0))); - problem.AddResidualBlock(cost, - nullptr, - &g[o.img_id], - &mosaicity[o.img_id], - &Itrue[o.hkl_slot], - &wedge); - } - break; - case GlobalScaleOptions::PartialityModel::Still: { - auto *cost = new ceres::AutoDiffCostFunction( - new IntensityStillResidual(*o.r, o.sigma)); - problem.AddResidualBlock(cost, - nullptr, - &g[o.img_id], - &R_sq[o.img_id], - &Itrue[o.hkl_slot]); - } - break; - case GlobalScaleOptions::PartialityModel::Unity: { - auto *cost = new ceres::AutoDiffCostFunction( - new IntensityFixedResidual(*o.r, o.sigma, 1.0)); - problem.AddResidualBlock(cost, - nullptr, - &g[o.img_id], - &Itrue[o.hkl_slot]); - } - break; - case GlobalScaleOptions::PartialityModel::Fixed: { - auto *cost = new ceres::AutoDiffCostFunction( - new IntensityFixedResidual(*o.r, o.sigma, o.r->partiality)); - problem.AddResidualBlock(cost, - nullptr, - &g[o.img_id], - &Itrue[o.hkl_slot]); - } - break; - } - } - - for (int i = 0; i < g.size(); ++i) { - if (image_slot_used[i]) { - auto *cost = new ceres::AutoDiffCostFunction( - new ScaleRegularizationResidual(0.05)); - problem.AddResidualBlock(cost, nullptr, &g[i]); - } - } - - if (rotation_crystallography) { - if (opt.smoothen_g) { - for (int i = 0; i < g.size() - 2; ++i) { - if (image_slot_used[i] && image_slot_used[i + 1] && image_slot_used[i + 2]) { - auto *cost = new ceres::AutoDiffCostFunction( - new SmoothnessRegularizationResidual(0.05)); - - problem.AddResidualBlock(cost, nullptr, &g[i], &g[i + 1], &g[i + 2]); - } - } - } - - if (opt.smoothen_mos && opt.partiality_model == GlobalScaleOptions::PartialityModel::Rotation) { - for (int i = 0; i < mosaicity.size() - 2; ++i) { - if (image_slot_used[i] && image_slot_used[i + 1] && image_slot_used[i + 2]) { - auto *cost = new ceres::AutoDiffCostFunction( - new SmoothnessRegularizationResidual(0.05)); - problem.AddResidualBlock(cost, nullptr, &mosaicity[i], &mosaicity[i + 1], &mosaicity[i + 2]); - } - } - } - } - - if (opt.partiality_model == GlobalScaleOptions::PartialityModel::Still) { - for (int i = 0; i < R_sq.size(); ++i) { - if (image_slot_used[i]) { - problem.SetParameterLowerBound(&R_sq[i], 0, 1e-9); - problem.SetParameterUpperBound(&R_sq[i], 0, 1.0); - } - } - } - - // Scaling factors must be always positive - for (int i = 0; i < g.size(); i++) { - if (image_slot_used[i]) - problem.SetParameterLowerBound(&g[i], 0, 1e-12); - } - - // Mosaicity refinement + bounds - if (opt.partiality_model == GlobalScaleOptions::PartialityModel::Rotation) { - for (int i = 0; i < mosaicity.size(); ++i) { - if (image_slot_used[i]) { - problem.SetParameterLowerBound(&mosaicity[i], 0, opt.mosaicity_min_deg); - problem.SetParameterUpperBound(&mosaicity[i], 0, opt.mosaicity_max_deg); - } - } - if (!opt.refine_wedge) - problem.SetParameterBlockConstant(&wedge); - else - problem.SetParameterLowerBound(&wedge, 0, 0.0); - } - - // use all available threads - unsigned int hw = std::thread::hardware_concurrency(); - if (hw == 0) - hw = 1; // fallback - - ceres::Solver::Options options; - - options.linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY; - options.minimizer_progress_to_stdout = false; - options.max_num_iterations = opt.max_num_iterations; - options.max_solver_time_in_seconds = opt.max_solver_time_s; - options.num_threads = static_cast(hw); - options.function_tolerance = 1e-4; - - ceres::Solver::Summary summary; - ceres::Solve(options, &problem, &summary); - std::cout << summary.BriefReport() << std::endl; - } - - void proc_obs(const std::vector > &observations, - const GlobalScaleOptions &opt, - std::vector &image_slot_used, - std::vector &obs, - std::unordered_map &hklToSlot - ) { - for (int i = 0; i < observations.size(); i++) { - for (const auto &r: observations[i]) { - const double d = SafeD(r.d); - if (!std::isfinite(d)) - continue; - if (!std::isfinite(r.I)) - continue; - - if (opt.d_min_limit_A > 0.0 && d < opt.d_min_limit_A) - continue; - - if (!std::isfinite(r.zeta) || r.zeta <= 0.0f) - continue; - if (!std::isfinite(r.rlp) || r.rlp == 0.0f) - continue; - - const double sigma = SafeSigma(r.sigma, opt.min_sigma); - - const int img_id = i / opt.image_cluster; - image_slot_used[img_id] = 1; - - int hkl_slot; - try { - const HKLKey key = CanonicalizeHKLKey(r, opt); - auto it = hklToSlot.find(key); - if (it == hklToSlot.end()) { - hkl_slot = static_cast(hklToSlot.size()); - hklToSlot.emplace(key, hkl_slot); - } else { - hkl_slot = it->second; - } - } catch (...) { - continue; - } - - ObsRef o; - o.r = &r; - o.img_id = img_id; - o.hkl_slot = hkl_slot; - o.sigma = sigma; - obs.push_back(o); - } - } - } -} // namespace - -GlobalScaleResult GlobalScaleCeres(const std::vector > &observations, - const GlobalScaleOptions &opt) { - if (opt.image_cluster <= 0) - throw std::invalid_argument("image_cluster must be positive"); - - const bool rotation_crystallography = opt.wedge_deg.has_value(); - - size_t nrefl = 0; - for (const auto &i: observations) - nrefl += i.size(); - - std::vector obs; - obs.reserve(nrefl); - - std::unordered_map hklToSlot; - hklToSlot.reserve(nrefl); - - size_t n_image_slots = observations.size() / opt.image_cluster + - (observations.size() % opt.image_cluster > 0 ? 1 : 0); - - std::vector image_slot_used(n_image_slots, 0); - - proc_obs(observations, opt, image_slot_used, obs, hklToSlot); - - const int nhkl = static_cast(hklToSlot.size()); - - std::vector g(n_image_slots, 1.0); - std::vector mosaicity(n_image_slots, opt.mosaicity_init_deg); - std::vector R_sq(n_image_slots, 0.001 * 0.001); - - for (int i = 0; i < n_image_slots; i++) { - if (!image_slot_used[i]) { - mosaicity[i] = NAN; - g[i] = NAN; - R_sq[i] = NAN; - } else if (opt.mosaicity_init_deg_vec.size() > i && std::isfinite(opt.mosaicity_init_deg_vec[i])) { - mosaicity[i] = opt.mosaicity_init_deg_vec[i]; - } - } - - scale(opt, g, mosaicity, R_sq, image_slot_used, rotation_crystallography, nhkl, obs); - - GlobalScaleResult out; - out.image_scale_g.resize(observations.size(), NAN); - out.mosaicity_deg.resize(observations.size(), NAN); - for (int i = 0; i < observations.size(); i++) { - size_t img_slot = i / opt.image_cluster; - if (image_slot_used[img_slot]) { - out.image_scale_g[i] = static_cast(g[img_slot]); - out.mosaicity_deg[i] = static_cast(mosaicity[img_slot]); - } - } - return out; -} diff --git a/image_analysis/scale_merge/GlobalScale.h b/image_analysis/scale_merge/GlobalScale.h deleted file mode 100644 index 85618efe..00000000 --- a/image_analysis/scale_merge/GlobalScale.h +++ /dev/null @@ -1,54 +0,0 @@ -// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#pragma once - -#include -#include -#include - -#include "../../common/Reflection.h" -#include "gemmi/symmetry.hpp" - -// Global (joint) scaling: instead of the alternating merge-then-scale-each-image loop, refine all -// per-image scales G_i, per-image mosaicities, the shared per-HKL true intensities Itrue and the -// global wedge together in a single Ceres problem. Returns only the per-image scale factors; the -// merge afterwards is done by MergeOnTheFly (Merge.h), the same as the alternating path. -struct GlobalScaleOptions { - int max_num_iterations = 100; - double max_solver_time_s = 60.0; - double min_sigma = 1e-3; - - // Symmetry canonicalization of HKL prior to scaling. If not set, raw HKL are used as-is. - std::optional space_group; - - // If true, treat Friedel mates as equivalent. If false, keep them separate (anomalous). - bool merge_friedel = true; - - // Rotation (wedge) partiality model. Set the wedge to enable it; leave unset for stills. - std::optional wedge_deg; - - double mosaicity_init_deg = 0.17; - double mosaicity_min_deg = 1e-3; - double mosaicity_max_deg = 2.0; - std::vector mosaicity_init_deg_vec; - - bool smoothen_g = true; - bool smoothen_mos = true; - - double d_min_limit_A = 0.0; - - int64_t image_cluster = 1; - - bool refine_wedge = false; - - enum class PartialityModel {Fixed, Rotation, Unity, Still} partiality_model = PartialityModel::Fixed; -}; - -struct GlobalScaleResult { - std::vector image_scale_g; ///< per input image (NaN if the image had no usable reflections) - std::vector mosaicity_deg; ///< per input image (NaN unless the rotation model is used) -}; - -GlobalScaleResult GlobalScaleCeres(const std::vector>& observations, - const GlobalScaleOptions& opt = {}); diff --git a/process/JFJochProcess.cpp b/process/JFJochProcess.cpp index 27cc22b1..59070af3 100644 --- a/process/JFJochProcess.cpp +++ b/process/JFJochProcess.cpp @@ -27,7 +27,6 @@ #include "../image_analysis/image_preprocessing/ImagePreprocessorCPU.h" #include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h" #include "../image_analysis/scale_merge/Merge.h" -#include "../image_analysis/scale_merge/GlobalScale.h" #include "../image_analysis/scale_merge/ScaleOnTheFly.h" #include "../image_analysis/scale_merge/SearchSpaceGroup.h" #include "../image_analysis/scale_merge/Combine3D.h" @@ -57,69 +56,6 @@ namespace { return ret; } - // Global (joint) scaling: refine all per-image scales/mosaicities and the shared per-HKL - // intensities in one Ceres problem (GlobalScaleCeres), then write the result back into each - // reflection's image_scale_corr exactly as ScaleOnTheFly does - recompute the rotation - // partiality from the refined per-image mosaicity, then image_scale_corr = rlp/(partiality*G). - // The downstream rot3d combine and MergeOnTheFly are unchanged, so every metric stays - // comparable with the alternating path; only how G and mosaicity are found differs. - void RunGlobalScaling(const DiffractionExperiment &experiment, - std::vector &outcomes, Logger &logger) { - const auto &s = experiment.GetScalingSettings(); - const bool rotation = experiment.GetPartialityModel() == PartialityModel::Rotation; - - std::vector> observations; - observations.reserve(outcomes.size()); - for (const auto &o: outcomes) - observations.push_back(o.reflections); - - GlobalScaleOptions opt; - opt.merge_friedel = s.GetMergeFriedel(); - if (const auto sg = experiment.GetGemmiSpaceGroup()) - opt.space_group = *sg; - opt.d_min_limit_A = s.GetHighResolutionLimit_A().value_or(0.0); - opt.mosaicity_init_deg = s.GetDefaultMosaicity(); - // The joint problem is large (per-image G + mosaicity + all Itrue + wedge); give the solver - // enough room to converge rather than the per-image default. - opt.max_solver_time_s = 300.0; - opt.max_num_iterations = 50; - - double wedge = 0.0; - if (rotation) { - opt.partiality_model = GlobalScaleOptions::PartialityModel::Rotation; - wedge = experiment.GetRotationWedgeForScaling().value_or(0.0); - opt.wedge_deg = wedge; - opt.mosaicity_min_deg = s.GetMinMosaicity(); - opt.mosaicity_max_deg = s.GetMaxMosaicity(); - } else { - opt.partiality_model = GlobalScaleOptions::PartialityModel::Fixed; - } - - const auto result = GlobalScaleCeres(observations, opt); - - const double half_wedge = wedge / 2.0; - for (size_t i = 0; i < outcomes.size(); ++i) { - const double G = result.image_scale_g[i]; - const double mos = result.mosaicity_deg[i]; - for (auto &r: outcomes[i].reflections) { - if (rotation && std::isfinite(r.delta_phi_deg) && std::isfinite(r.zeta) && mos > 1e-6) { - const double c1 = r.zeta / std::sqrt(2.0); - r.partiality = static_cast( - (std::erf((r.delta_phi_deg + half_wedge) * c1 / mos) - - std::erf((r.delta_phi_deg - half_wedge) * c1 / mos)) / 2.0); - } - const double denom = r.partiality * G; - r.image_scale_corr = (std::isfinite(r.rlp) && std::isfinite(denom) && denom > 0.0) - ? static_cast(r.rlp / denom) : NAN; - } - if (std::isfinite(G)) - outcomes[i].image_scale_g = static_cast(G); - if (rotation && std::isfinite(mos)) - outcomes[i].mosaicity_deg = static_cast(mos); - } - logger.Info("Global scaling complete ({} images)", outcomes.size()); - } - // XDS-order scaling. The rot3d combine emits fulls with partiality == 1 (image_scale_corr == 1), // so they were only ever scaled as per-frame *partials* upstream - their per-frame scale is // entangled with the rocking-curve/partiality model. This refits a per-frame scale directly on @@ -457,8 +393,21 @@ ProcessResult JFJochProcess::Run(JFJochProcessObserver *observer) { // Scaling and merging (full analysis only). if (full && !cancelled_ && result.indexing_rate.has_value() && result.indexing_rate > 0 && (config_.run_scaling || !config_.reference_data.empty())) { - if (observer) - observer->OnPhase("Scaling and merging"); + // Scaling/merging is a long post-pass; report each sub-step as a phase so the GUI progress + // bar reflects what is happening instead of freezing on one label. Also time each phase + // (logged on transition) so the bottlenecks are visible. + auto t_phase = std::chrono::steady_clock::now(); + std::string prev_phase; + auto phase = [&](const std::string &p) { + const auto now = std::chrono::steady_clock::now(); + if (!prev_phase.empty()) + logger.Info("[timing] {}: {:.2f} s", prev_phase, + std::chrono::duration(now - t_phase).count()); + t_phase = now; + prev_phase = p; + if (observer) observer->OnPhase(p); + }; + phase("Scaling and merging"); const bool pixel_refine_path = experiment_.GetIndexingSettings().GetGeomRefinementAlgorithm() == GeomRefinementAlgorithmEnum::PixelRefine; @@ -466,31 +415,41 @@ ProcessResult JFJochProcess::Run(JFJochProcessObserver *observer) { // ScaleOnTheFly is only for the classical, no-reference path; with a reference (or // PixelRefine) each image is already scaled, so we merge directly. if (config_.reference_data.empty() && !pixel_refine_path) { - if (config_.global_scaling) { - logger.Info("Running global (joint) scaling ..."); - RunGlobalScaling(experiment_, indexer->GetIntegrationOutcome(), logger); - } else { - logger.Info("Running scaling ..."); - ScalingResult scale_result(0); - for (int i = 0; i < config_.scaling_iter; i++) { - auto merge_result = MergeAll(experiment_, indexer->GetIntegrationOutcome(), false); - scale_result = indexer->ScaleAllImages(merge_result); - } + logger.Info("Running scaling ..."); + ScalingResult scale_result(0); + double t_merge_all = 0.0, t_scale = 0.0; + for (int i = 0; i < config_.scaling_iter; i++) { + phase("Scaling images (iteration " + std::to_string(i + 1) + "/" + + std::to_string(config_.scaling_iter) + ")"); + const auto a = std::chrono::steady_clock::now(); + auto merge_result = MergeAll(experiment_, indexer->GetIntegrationOutcome(), false); + const auto b = std::chrono::steady_clock::now(); + scale_result = indexer->ScaleAllImages(merge_result); + const auto c = std::chrono::steady_clock::now(); + t_merge_all += std::chrono::duration(b - a).count(); + t_scale += std::chrono::duration(c - b).count(); } + logger.Info("[timing] scaling loop ({} iter): MergeAll(serial) {:.2f} s, ScaleAllImages(parallel) {:.2f} s", + config_.scaling_iter, t_merge_all, t_scale); } // -P rot3d: weight-sum each reflection's per-frame partials into one full before merging, so // the error model sees counting statistics (high ISa) instead of rocking-curve slicing scatter. const bool rot3d = experiment_.GetScalingSettings().GetCombine3D(); std::vector combined; - if (rot3d) + if (rot3d) { + phase("Combining 3D partials"); combined = CombineRotationObservations(indexer->GetIntegrationOutcome(), experiment_, &logger, config_.observation_dump_path); - if (rot3d && config_.scale_fulls) + } + if (rot3d && experiment_.GetScalingSettings().GetScaleFulls()) { + phase("Scaling fulls (XDS order)"); ScaleFulls(experiment_, combined, static_cast(config_.scaling_iter), config_.nthreads, logger); + } const std::vector &merge_input = rot3d ? combined : indexer->GetIntegrationOutcome(); + phase("Merging"); MergeOnTheFly merge_engine(experiment_); if (result.consensus_cell.has_value()) merge_engine.ReferenceCell(*result.consensus_cell); @@ -503,6 +462,7 @@ ProcessResult JFJochProcess::Run(JFJochProcessObserver *observer) { merge_engine.AddImage(outcome); auto merged_reflections = merge_engine.ExportReflections(); + phase("Computing statistics"); auto merged_statistics = merge_engine.MergeStats(merged_reflections, merge_input, config_.reference_data); logger.Info("Merge complete ({} unique reflections)", merged_reflections.size()); @@ -523,8 +483,16 @@ ProcessResult JFJochProcess::Run(JFJochProcessObserver *observer) { stats_text << merged_statistics; result.merge_statistics_text = stats_text.str(); - if (result.consensus_cell && write_files) + result.has_merge_statistics = true; + result.merge_statistics = merged_statistics; + result.error_model_isa = merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0; + result.has_reference = !config_.reference_data.empty(); + + if (result.consensus_cell && write_files) { + phase("Writing reflections"); WriteReflections(merged_reflections, *result.consensus_cell, experiment_, config_.output_prefix); + } + phase(""); // flush the last phase's timing } if (writer) { diff --git a/process/JFJochProcess.h b/process/JFJochProcess.h index 755d01ab..8cb89bee 100644 --- a/process/JFJochProcess.h +++ b/process/JFJochProcess.h @@ -15,6 +15,7 @@ #include "../common/JFJochMessages.h" // DataMessage #include "../common/JFJochReceiverPlots.h" // MeanProcessingTime #include "../image_analysis/spot_finding/SpotFindingSettings.h" +#include "../image_analysis/scale_merge/Merge.h" // MergeStatistics class JFJochHDF5Reader; @@ -53,13 +54,6 @@ struct ProcessConfig { // per-image live scaling path when present. bool run_scaling = false; int64_t scaling_iter = 3; - // Global (joint) scaling instead of the alternating merge/scale loop: refine all per-image - // scales, mosaicities and the shared per-HKL intensities together in one Ceres problem. - bool global_scaling = false; - // XDS-order scaling: after the rot3d combine, refit a per-frame scale on the combined fulls - // (Unity model, no partiality term) so the per-frame scale is decoupled from rocking-curve - // model error. Only meaningful with -P rot3d. - bool scale_fulls = false; std::vector reference_data; // Diagnostic: if set, the -P rot3d combine writes the unmerged fulls here (for comparison vs XDS). @@ -78,6 +72,13 @@ struct ProcessResult { MeanProcessingTime mean_processing_time{}; std::optional written_master_path; std::string merge_statistics_text; // populated when scaling/merging ran + + // Structured merge statistics (per-shell + overall), populated when scaling/merging ran. ISa is + // the error-model asymptotic I/sigma (1/b); has_reference is true when a reference MTZ drove CCref. + bool has_merge_statistics = false; + MergeStatistics merge_statistics; + double error_model_isa = 0.0; + bool has_reference = false; }; // Callbacks for progress and live results. Methods may be called from worker threads, so an diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index 56011515..dfc1310e 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -49,7 +49,6 @@ void print_usage() { std::cout << " Scaling and merging" << std::endl; std::cout << " -M, --scale-merge Scale and merge (refine mosaicity) and write scaled.hkl + image.dat" << std::endl; - std::cout << " --global-scale Joint global scaling (all images + Itrue in one fit) instead of the alternating loop; implies -M" << std::endl; std::cout << " --scale-fulls After -P rot3d combine, refit a per-frame scale on the fulls (XDS order, Unity model); implies -M" << std::endl; std::cout << " -P, --partiality Partiality model fixed|rot|rot3d|unity (default: fixed). rot3d = rot + 3D combine of per-frame partials" << std::endl; std::cout << " -A, --anomalous Anomalous mode (don't merge Friedel pairs)" << std::endl; @@ -94,7 +93,6 @@ enum { OPT_REFERENCE_COLUMN, OPT_DUMP_OBSERVATIONS, OPT_INTEGRATOR, - OPT_GLOBAL_SCALE, OPT_SCALE_FULLS }; @@ -116,7 +114,6 @@ static option long_options[] = { {"refine-bfactor", no_argument, nullptr, 'B'}, {"wedge", optional_argument, nullptr, 'w'}, {"scale-merge", no_argument, nullptr, 'M'}, - {"global-scale", no_argument, nullptr, OPT_GLOBAL_SCALE}, {"scale-fulls", no_argument, nullptr, OPT_SCALE_FULLS}, {"refine", required_argument, nullptr, 'r'}, @@ -277,7 +274,6 @@ int main(int argc, char **argv) { int rotation_indexing_image_count = 100; std::optional rotation_indexing_range; bool run_scaling = false; - bool global_scaling = false; bool scale_fulls = false; bool anomalous_mode = false; std::optional space_group_number; @@ -508,10 +504,6 @@ int main(int argc, char **argv) { case 'M': run_scaling = true; break; - case OPT_GLOBAL_SCALE: - run_scaling = true; - global_scaling = true; - break; case OPT_SCALE_FULLS: run_scaling = true; scale_fulls = true; @@ -722,6 +714,7 @@ int main(int argc, char **argv) { ScalingSettings scaling_settings; scaling_settings.SetPartialityModel(partiality_model); scaling_settings.Combine3D(combine_3d); + scaling_settings.ScaleFulls(scale_fulls); if (d_min_scale_merge) scaling_settings.HighResolutionLimit_A(d_min_scale_merge.value()); scaling_settings.MergeFriedel(!anomalous_mode); @@ -798,8 +791,6 @@ int main(int argc, char **argv) { config.rotation_indexing_image_count = rotation_indexing_image_count; config.forced_rotation_lattice = forced_rotation_lattice; config.run_scaling = run_scaling; - config.global_scaling = global_scaling; - config.scale_fulls = scale_fulls; config.scaling_iter = scaling_iter; config.reference_data = reference_data; config.observation_dump_path = dump_observations;