From b53f0d6474a33e6fffe337554a3b5e5e9cc72d74 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Mon, 11 May 2026 10:37:39 +0200 Subject: [PATCH] jfjoch_process: Print statistics moved to Merge.h/Merge.cpp and ScalingResult.cpp/ScalingResult.h --- common/ScalingSettings.cpp | 4 ++ common/ScalingSettings.h | 3 +- image_analysis/scale_merge/CMakeLists.txt | 3 +- image_analysis/scale_merge/Merge.cpp | 21 +++++++ image_analysis/scale_merge/Merge.h | 2 + image_analysis/scale_merge/ScaleOnTheFly.cpp | 15 +++-- image_analysis/scale_merge/ScalingResult.cpp | 38 ++++++++++++ image_analysis/scale_merge/ScalingResult.h | 9 +-- tools/jfjoch_process.cpp | 64 +++----------------- 9 files changed, 89 insertions(+), 70 deletions(-) create mode 100644 image_analysis/scale_merge/ScalingResult.cpp diff --git a/common/ScalingSettings.cpp b/common/ScalingSettings.cpp index 7969c331..04f0f8b1 100644 --- a/common/ScalingSettings.cpp +++ b/common/ScalingSettings.cpp @@ -77,6 +77,10 @@ double ScalingSettings::GetMaxWedge() const { return 10.0; } +double ScalingSettings::GetDefaultMosaicity() const { + return 0.1; +} + ScalingSettings &ScalingSettings::RotationWedgeForScaling(std::optional input) { if (input) { // TODO: Use fmt diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h index fae31415..941b6dea 100644 --- a/common/ScalingSettings.h +++ b/common/ScalingSettings.h @@ -6,7 +6,7 @@ #include #include "JFJochException.h" -enum class PartialityModel { Fixed, Rotation, Unity, Still }; +enum class PartialityModel { Fixed, Rotation, Unity }; class ScalingSettings { std::optional partiality_mode; @@ -35,6 +35,7 @@ public: [[nodiscard]] double GetMaxB() const; [[nodiscard]] double GetMinMosaicity() const; + [[nodiscard]] double GetDefaultMosaicity() const; [[nodiscard]] double GetMaxMosaicity() const; [[nodiscard]] double GetMinWedge() const; diff --git a/image_analysis/scale_merge/CMakeLists.txt b/image_analysis/scale_merge/CMakeLists.txt index a4b0a686..2ad08bcf 100644 --- a/image_analysis/scale_merge/CMakeLists.txt +++ b/image_analysis/scale_merge/CMakeLists.txt @@ -7,5 +7,6 @@ ADD_LIBRARY(JFJochScaleMerge FrenchWilson.cpp FrenchWilson.h ScaleOnTheFly.h HKLKey.cpp HKLKey.h - ScalingResult.h) + ScalingResult.h + ScalingResult.cpp) TARGET_LINK_LIBRARIES(JFJochScaleMerge Ceres::ceres Eigen3::Eigen JFJochCommon) \ No newline at end of file diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index cb68acd1..f57d8b0d 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -293,3 +293,24 @@ MergeResult MergeReflections(const std::vector > &observ return out; } + +void MergeStatistics::Print(Logger &logger) const { + logger.Info(""); + logger.Info(" {:>8s} {:>8s} {:>8s} {:>8s}", "d_min", "N_obs", "N_uniq", ""); + logger.Info(" {:->8s} {:->8s} {:->8s} {:->8s}", "", "", "", "", "", ""); + for (const auto &sh: shells) { + if (sh.unique_reflections == 0) + continue; + logger.Info(" {:8.2f} {:8d} {:8d} {:8.1f}", + sh.d_min, sh.total_observations, sh.unique_reflections, + sh.mean_i_over_sigma); + } + { + const auto &ov = overall; + logger.Info(" {:->8s} {:->8s} {:->8s} {:->8s}", "", "", "", ""); + logger.Info(" {:>8s} {:8d} {:8d} {:8.1f}", + "Overall", ov.total_observations, ov.unique_reflections, + ov.mean_i_over_sigma); + } + logger.Info(""); +} diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 7d249fcf..8af88753 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -7,6 +7,7 @@ #include #include +#include "../../common/Logger.h" #include "../../common/DiffractionExperiment.h" #include "../../common/Reflection.h" #include "gemmi/symmetry.hpp" @@ -26,6 +27,7 @@ struct MergeStatisticsShell { struct MergeStatistics { std::vector shells; MergeStatisticsShell overall; + void Print(Logger &logger) const; }; struct MergeResult { diff --git a/image_analysis/scale_merge/ScaleOnTheFly.cpp b/image_analysis/scale_merge/ScaleOnTheFly.cpp index 79676491..45f96f6a 100644 --- a/image_analysis/scale_merge/ScaleOnTheFly.cpp +++ b/image_analysis/scale_merge/ScaleOnTheFly.cpp @@ -96,8 +96,6 @@ bool ScaleOnTheFly::Accept(const Reflection &r) { switch (model) { case PartialityModel::Rotation: return std::isfinite(r.zeta) && r.zeta > 0.0f; - case PartialityModel::Still: - return std::isfinite(r.dist_ewald); case PartialityModel::Fixed: case PartialityModel::Unity: return true; @@ -113,11 +111,17 @@ ScaleOnTheFlyResult ScaleOnTheFly::Scale(std::vector &reflections, s ScaleOnTheFlyResult result{ .B = 0.0, - .G = 1.0, - .mos = mosaicity_deg.value_or(0.0), - .wedge = rot_wedge_deg.value_or(0.0) + .G = 1.0 }; + if (model == PartialityModel::Rotation) { + result.mos = mosaicity_deg.value_or(s.GetDefaultMosaicity()); + result.wedge = rot_wedge_deg.value_or(0.0); + } else { + result.mos = NAN; + result.wedge = NAN; + } + size_t n_reflections = 0; HKLKeyGenerator key_generator(s.GetMergeFriedel(), sg); for (const auto &r: reflections) { @@ -207,6 +211,7 @@ ScaleOnTheFlyResult ScaleOnTheFly::Scale(std::vector &reflections, s break; } default: + // For fixed partiality there is no need to change anything break; } r.scaling_correction = static_cast(r.rlp / (B_term * r.partiality * result.G)); diff --git a/image_analysis/scale_merge/ScalingResult.cpp b/image_analysis/scale_merge/ScalingResult.cpp new file mode 100644 index 00000000..8c683da9 --- /dev/null +++ b/image_analysis/scale_merge/ScalingResult.cpp @@ -0,0 +1,38 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "ScalingResult.h" +#include +#include + +#include "../../common/JFJochException.h" + +ScalingResult::ScalingResult(size_t n) + : image_scale_g(n, NAN), + mosaicity_deg(n, NAN), + image_bfactor_Ang2(n, NAN), + rotation_wedge_deg(n, NAN) { +} + +void ScalingResult::SaveToFile(const std::string &filename) { + const std::string img_path = filename + "_image.dat"; + std::ofstream img_file(img_path); + if (!img_file) { + throw JFJochException(JFJochExceptionCategory::FileWriteError + , "Cannot open {} for writing"); + } + + img_file << "# image_id G B mosaicity_deg wedge_deg\n"; + + for (size_t i = 0; i < image_scale_g.size(); ++i) { + img_file << i + << " " << image_scale_g[i] + << " " << image_bfactor_Ang2[i] + << " " << mosaicity_deg[i] + << " " << rotation_wedge_deg[i] + << "\n"; + } + + img_file.close(); +} + diff --git a/image_analysis/scale_merge/ScalingResult.h b/image_analysis/scale_merge/ScalingResult.h index 6b97304e..880f5c55 100644 --- a/image_analysis/scale_merge/ScalingResult.h +++ b/image_analysis/scale_merge/ScalingResult.h @@ -4,6 +4,7 @@ #pragma once #include +#include struct ScalingResult { std::vector image_scale_g; @@ -11,10 +12,6 @@ struct ScalingResult { std::vector image_bfactor_Ang2; std::vector rotation_wedge_deg; - explicit ScalingResult(size_t n) - : image_scale_g(n, NAN), - mosaicity_deg(n, NAN), - image_bfactor_Ang2(n, NAN), - rotation_wedge_deg(n, NAN) { - } + explicit ScalingResult(size_t n); + void SaveToFile(const std::string &filename); }; diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index 8e93cbfe..20e67be5 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -581,7 +581,7 @@ int main(int argc, char **argv) { auto merge_result = indexer.Merge(); auto scale_end = std::chrono::steady_clock::now(); double scale_time = std::chrono::duration(scale_end - scale_start).count(); - /* + if (!fixed_space_group) { logger.Info("Searching for space group from P1-merged reflections ..."); @@ -608,7 +608,7 @@ int main(int argc, char **argv) { } } logger.Info(""); - + /* if (sg_search.best_space_group.has_value()) { logger.Info("Re-running scaling in detected space group {}", sg_search.best_space_group->short_name()); @@ -626,8 +626,8 @@ int main(int argc, char **argv) { } } else { logger.Warning("No space group accepted; keeping P1-merged result"); - } - } */ + } */ + } end_msg.image_scale_factor = scale_result.image_scale_g; @@ -635,59 +635,9 @@ int main(int argc, char **argv) { scale_time, merge_result.merged.size()); // Print resolution-shell statistics table - { - const auto &stats = merge_result.statistics; - logger.Info(""); - logger.Info(" {:>8s} {:>8s} {:>8s} {:>8s}", - "d_min", "N_obs", "N_uniq", ""); - logger.Info(" {:->8s} {:->8s} {:->8s} {:->8s}", - "", "", "", "", "", ""); - for (const auto &sh: stats.shells) { - if (sh.unique_reflections == 0) - continue; - logger.Info(" {:8.2f} {:8d} {:8d} {:8.1f}", - sh.d_min, sh.total_observations, sh.unique_reflections, - sh.mean_i_over_sigma); - } - { - const auto &ov = stats.overall; - logger.Info(" {:->8s} {:->8s} {:->8s} {:->8s}", - "", "", "", ""); - logger.Info(" {:>8s} {:8d} {:8d} {:8.1f}", - "Overall", ov.total_observations, ov.unique_reflections, - ov.mean_i_over_sigma); - } - logger.Info(""); - } - - { - const std::string img_path = output_prefix + "_image.dat"; - std::ofstream img_file(img_path); - if (!img_file) { - logger.Error("Cannot open {} for writing", img_path); - } else { - if (experiment.GetPartialityModel() == PartialityModel::Rotation) { - img_file << "# image_id G B mosaicity_deg wedge_deg\n"; - for (size_t i = 0; i < scale_result.image_scale_g.size(); ++i) { - img_file << i - << " " << scale_result.image_scale_g[i] - << " " << scale_result.image_bfactor_Ang2[i] - << " " << scale_result.mosaicity_deg[i] - << " " << scale_result.rotation_wedge_deg[i] - << "\n"; - } - } else { - img_file << "# image_id G B\n"; - for (size_t i = 0; i < scale_result.image_scale_g.size(); ++i) { - img_file << i - << " " << scale_result.image_scale_g[i] - << " " << scale_result.image_bfactor_Ang2[i] - << "\n"; - } - } - img_file.close(); - } - } + const auto &stats = merge_result.statistics; + stats.Print(logger); + scale_result.SaveToFile(output_prefix + "_scale.dat"); { FrenchWilsonOptions fw_opts;