From 50ebf1694ed66871ebca356cb46a94d5f5c406c6 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Thu, 25 Jun 2026 07:11:39 +0200 Subject: [PATCH] Add global (joint) scaling path: jfjoch_process --global-scale Refine all per-image scales, per-image mosaicities and the shared per-HKL true intensities (plus the global wedge) together in one Ceres problem (GlobalScale.{h,cpp}, GlobalScaleCeres), as an alternative to the existing alternating merge-then-scale-each-image loop. The result feeds back into each reflection's image_scale_corr exactly as ScaleOnTheFly does, so the rot3d combine and MergeOnTheFly run unchanged and every metric stays comparable. On HEWL crystal 2 the joint fit converges to the same place as the alternating loop (anomalous S-peak vs XDS 0.53x -> 0.54x, ISa 9.4 both), confirming the alternating path already reaches the joint optimum. Kept as a validated alternative and the home for future global corrections. A shared quadratic absorption surface in detector position was prototyped here and dropped: it fit large non-physical coefficients (radially degenerate with the per-HKL resolution structure), lowered the scaling-fit residual but raised the error-model b (ISa 9.4 -> 7.3) and did not improve anomalous accuracy. Co-Authored-By: Claude Opus 4.8 --- image_analysis/scale_merge/CMakeLists.txt | 2 + image_analysis/scale_merge/GlobalScale.cpp | 489 +++++++++++++++++++++ image_analysis/scale_merge/GlobalScale.h | 54 +++ process/JFJochProcess.cpp | 79 +++- process/JFJochProcess.h | 3 + tools/jfjoch_process.cpp | 11 +- 6 files changed, 632 insertions(+), 6 deletions(-) create mode 100644 image_analysis/scale_merge/GlobalScale.cpp create mode 100644 image_analysis/scale_merge/GlobalScale.h diff --git a/image_analysis/scale_merge/CMakeLists.txt b/image_analysis/scale_merge/CMakeLists.txt index 190040c7..6681329d 100644 --- a/image_analysis/scale_merge/CMakeLists.txt +++ b/image_analysis/scale_merge/CMakeLists.txt @@ -5,6 +5,8 @@ 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 new file mode 100644 index 00000000..c01c322b --- /dev/null +++ b/image_analysis/scale_merge/GlobalScale.cpp @@ -0,0 +1,489 @@ +// 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 new file mode 100644 index 00000000..85618efe --- /dev/null +++ b/image_analysis/scale_merge/GlobalScale.h @@ -0,0 +1,54 @@ +// 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 74e4004e..27a5da27 100644 --- a/process/JFJochProcess.cpp +++ b/process/JFJochProcess.cpp @@ -27,6 +27,7 @@ #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/SearchSpaceGroup.h" #include "../image_analysis/scale_merge/Combine3D.h" #include "../image_analysis/WriteReflections.h" @@ -54,6 +55,69 @@ namespace { ret.assign(unique_ordinals.begin(), unique_ordinals.end()); 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()); + } } JFJochProcess::JFJochProcess(JFJochHDF5Reader &reader, DiffractionExperiment experiment, @@ -380,11 +444,16 @@ 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) { - 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); + 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); + } } } diff --git a/process/JFJochProcess.h b/process/JFJochProcess.h index 100f480f..d0634421 100644 --- a/process/JFJochProcess.h +++ b/process/JFJochProcess.h @@ -53,6 +53,9 @@ 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; std::vector reference_data; // Diagnostic: if set, the -P rot3d combine writes the unmerged fulls here (for comparison vs XDS). diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index 682c76e5..d387c862 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -49,6 +49,7 @@ 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 << " -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; std::cout << " -B, --refine-bfactor Refine per image B-factor" << std::endl; @@ -91,7 +92,8 @@ enum { OPT_REJECT_DELTA_CCHALF, OPT_REFERENCE_COLUMN, OPT_DUMP_OBSERVATIONS, - OPT_INTEGRATOR + OPT_INTEGRATOR, + OPT_GLOBAL_SCALE }; static option long_options[] = { @@ -112,6 +114,7 @@ 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}, {"refine", required_argument, nullptr, 'r'}, {"two-pass-rotation", optional_argument, nullptr, 'R'}, @@ -271,6 +274,7 @@ int main(int argc, char **argv) { int rotation_indexing_image_count = 30; std::optional rotation_indexing_range; bool run_scaling = false; + bool global_scaling = false; bool anomalous_mode = false; std::optional space_group_number; std::optional fixed_reference_unit_cell; @@ -500,6 +504,10 @@ 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_MIN_PARTIALITY: min_partiality = std::stod(optarg); break; @@ -782,6 +790,7 @@ 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.scaling_iter = scaling_iter; config.reference_data = reference_data; config.observation_dump_path = dump_observations;