From 1d3ccdaa00a967d241be4cfcf1260b6bde24c461 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sun, 17 May 2026 20:00:32 +0200 Subject: [PATCH] Add possible reflection and completeness to MergeStatistics --- common/UnitCell.cpp | 3 + common/UnitCell.h | 7 +- gemmi_gph/gemmi/reciproc.hpp | 57 +++++++++++++++ image_analysis/scale_merge/Merge.cpp | 105 +++++++++++++++++++++------ image_analysis/scale_merge/Merge.h | 5 ++ tools/jfjoch_process.cpp | 2 +- tools/jfjoch_scale.cpp | 2 +- 7 files changed, 153 insertions(+), 28 deletions(-) create mode 100644 gemmi_gph/gemmi/reciproc.hpp diff --git a/common/UnitCell.cpp b/common/UnitCell.cpp index 30da2878..70cce7b9 100644 --- a/common/UnitCell.cpp +++ b/common/UnitCell.cpp @@ -69,3 +69,6 @@ std::optional MeanUnitCell(const std::vector &cells) { return ret; } +UnitCell::operator gemmi::UnitCell() const { + return {a, b, c, alpha, beta, gamma}; +} \ No newline at end of file diff --git a/common/UnitCell.h b/common/UnitCell.h index 3aa196dc..cec7b0fa 100644 --- a/common/UnitCell.h +++ b/common/UnitCell.h @@ -6,6 +6,7 @@ #include #include #include +#include struct UnitCell { float a; @@ -15,8 +16,10 @@ struct UnitCell { float beta; float gamma; - bool is_finite() const; - bool is_close(const UnitCell &ref, float dist_tolerance, float angle_tolerance_deg) const; + operator gemmi::UnitCell() const; + + [[nodiscard]] bool is_finite() const; + [[nodiscard]] bool is_close(const UnitCell &ref, float dist_tolerance, float angle_tolerance_deg) const; }; std::ostream &operator<<(std::ostream &output, const UnitCell &in); diff --git a/gemmi_gph/gemmi/reciproc.hpp b/gemmi_gph/gemmi/reciproc.hpp new file mode 100644 index 00000000..6d5caf09 --- /dev/null +++ b/gemmi_gph/gemmi/reciproc.hpp @@ -0,0 +1,57 @@ +// Copyright 2020 Global Phasing Ltd. +// +// Reciprocal space helper functions. + +#ifndef GEMMI_RECIPROC_HPP_ +#define GEMMI_RECIPROC_HPP_ + +#include "symmetry.hpp" // for SpaceGroup +#include "unitcell.hpp" // for UnitCell + +namespace gemmi { + +// dmin should include a tiny margin for numerical errors +template +void for_all_reflections(Func func, + const UnitCell& cell, const SpaceGroup* spacegroup, + double dmin, double dmax=0., bool unique=true) { + Miller lim = cell.get_hkl_limits(dmin); + double inv_dmin2 = 1. / sq(dmin); + double inv_dmax2 = 0.; + if (dmax > 0) + inv_dmax2 = dmax == INFINITY ? -1 : 1. / sq(dmax); + ReciprocalAsu asu(spacegroup); + GroupOps gops = spacegroup->operations(); + Miller hkl; + for (hkl[0] = -lim[0]; hkl[0] <= lim[0]; ++hkl[0]) + for (hkl[1] = -lim[1]; hkl[1] <= lim[1]; ++hkl[1]) + for (hkl[2] = -lim[2]; hkl[2] <= lim[2]; ++hkl[2]) + if (!unique || asu.is_in(hkl)) { + double inv_d2 = cell.calculate_1_d2(hkl); + if (inv_d2 <= inv_dmin2 && inv_d2 > inv_dmax2 && + !gops.is_systematically_absent(hkl)) + func(hkl); + } +} + +// dmin should include a tiny margin for numerical errors +inline int count_reflections(const UnitCell& cell, const SpaceGroup* spacegroup, + double dmin, double dmax=0., bool unique=true) { + int counter = 0; + for_all_reflections([&counter](const Miller&) { ++counter; }, + cell, spacegroup, dmin, dmax, unique); + return counter; +} + +inline std::vector +make_miller_vector(const UnitCell& cell, const SpaceGroup* spacegroup, + double dmin, double dmax=0., bool unique=true) { + std::vector hkls; + for_all_reflections([&hkls](const Miller& hkl) { hkls.push_back(hkl); }, + cell, spacegroup, dmin, dmax, unique); + return hkls; +} + + +} // namespace gemmi +#endif diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 0483ee61..3c5e2fc4 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -9,6 +9,8 @@ #include #include +#include + #include "../../common/ResolutionShells.h" #include "HKLKey.h" @@ -165,11 +167,58 @@ std::vector MergeAll(const DiffractionExperiment &x, return out; } +struct ShellAccum { + int total_obs = 0; + int unique = 0; + int possible = 0; + + double sum_i_over_sigma = 0.0; + int n_i_over_sigma = 0; + + double sum_x = 0.0; + double sum_y = 0.0; + double sum_x2 = 0.0; + double sum_y2 = 0.0; + double sum_xy = 0.0; + int n_cc_half = 0; +}; + +void CalcPossibleReflections(const DiffractionExperiment &x, + const UnitCell &cell, + double d_min, + double d_max, + const ResolutionShells &shells, + std::vector &acc) { + gemmi::UnitCell gemmi_cell = cell; + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(x.GetSpaceGroupNumber().value_or(1)); + + // Generate unique reflections + std::vector possible_hkls = gemmi::make_miller_vector(gemmi_cell, sg, d_min, d_max, true); + CrystalLattice lattice(cell); + const auto astar = lattice.Astar(); + const auto bstar = lattice.Bstar(); + const auto cstar = lattice.Cstar(); + + for (const auto& hkl: possible_hkls) { + const auto q = hkl[0] * astar + hkl[1] * bstar + hkl[2] * cstar; + const auto qlen = q.Length(); + if (qlen < 1e-6) + continue; + const auto d = 1.0 / qlen; + const auto shell = shells.GetShell(d); + if (!shell.has_value()) + continue; + const int s = *shell; + if (s >= 0 && s < acc.size()) + acc[s].possible++; + } +} + MergeStatistics MergeStats(const DiffractionExperiment &x, const std::vector &merged, const std::vector > &reflections, + const UnitCell &cell, const std::vector &merge_mask) { - if (!merge_mask.empty() && merge_mask.size() != reflections.size()) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Merge mask size mismatch"); @@ -204,22 +253,10 @@ MergeStatistics MergeStats(const DiffractionExperiment &x, const auto shell_mean_1_d2 = shells.GetShellMeanOneOverResSq(); const auto shell_min_res = shells.GetShellMinRes(); - struct ShellAccum { - int total_obs = 0; - int unique = 0; - double sum_i_over_sigma = 0.0; - int n_i_over_sigma = 0; - - double sum_x = 0.0; - double sum_y = 0.0; - double sum_x2 = 0.0; - double sum_y2 = 0.0; - double sum_xy = 0.0; - int n_cc_half = 0; - }; - std::vector acc(n_shells); + CalcPossibleReflections(x, cell, d_min_pad, d_max_pad, shells, acc); + for (const auto &m: merged) { const auto shell = shells.GetShell(m.d); if (!shell.has_value()) @@ -283,6 +320,7 @@ MergeStatistics MergeStats(const DiffractionExperiment &x, ss.d_max = s == 0 ? d_max_pad : shell_min_res[s - 1]; ss.total_observations = sa.total_obs; ss.unique_reflections = sa.unique; + ss.possible_unique_reflections = sa.possible; ss.mean_i_over_sigma = sa.n_i_over_sigma > 0 ? sa.sum_i_over_sigma / sa.n_i_over_sigma : 0.0; @@ -305,6 +343,7 @@ MergeStatistics MergeStats(const DiffractionExperiment &x, overall.d_min = d_min; overall.d_max = d_max; + int all_possible = 0; int all_unique = 0; double sum_i_over_sigma = 0.0; int n_i_over_sigma = 0; @@ -319,6 +358,7 @@ MergeStatistics MergeStats(const DiffractionExperiment &x, for (const auto &sa: acc) { overall.total_observations += sa.total_obs; all_unique += sa.unique; + all_possible += sa.possible; sum_i_over_sigma += sa.sum_i_over_sigma; n_i_over_sigma += sa.n_i_over_sigma; @@ -330,6 +370,7 @@ MergeStatistics MergeStats(const DiffractionExperiment &x, all_n_cc_half += sa.n_cc_half; } + overall.possible_unique_reflections = all_possible; overall.unique_reflections = all_unique; overall.mean_i_over_sigma = n_i_over_sigma > 0 ? sum_i_over_sigma / n_i_over_sigma : 0.0; @@ -351,21 +392,37 @@ MergeStatistics MergeStats(const DiffractionExperiment &x, void MergeStatistics::Print(Logger &logger) const { logger.Info(""); - logger.Info(" {:>8s} {:>8s} {:>8s} {:>8s} {:>8s}", "d_min", "N_obs", "N_uniq", "", "CC1/2"); - logger.Info(" {:->8s} {:->8s} {:->8s} {:->8s} {:->8s}", "", "", "", "", ""); + logger.Info(" {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s}", "d_min", "N_obs", "N_uniq", "N_possib", "Compl","", "CC1/2"); + logger.Info(" {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s}", "", "", "", "", "", "", ""); for (const auto &sh: shells) { if (sh.unique_reflections == 0) continue; - logger.Info(" {:8.2f} {:8d} {:8d} {:8.1f} {:8.3f}", - sh.d_min, sh.total_observations, sh.unique_reflections, - sh.mean_i_over_sigma, sh.cc_half*100.0); + double completeness = sh.possible_unique_reflections > 0 + ? static_cast(sh.unique_reflections) / sh.possible_unique_reflections * 100.0 : 0.0; + + logger.Info(" {:8.2f} {:8d} {:8d} {:8d} {:7.1f}% {:8.1f} {:7.1f}%", + sh.d_min, + sh.total_observations, + sh.unique_reflections, + sh.possible_unique_reflections, + completeness, + sh.mean_i_over_sigma, + sh.cc_half*100.0); } { const auto &ov = overall; - logger.Info(" {:->8s} {:->8s} {:->8s} {:->8s} {:->8s}", "", "", "", "", ""); - logger.Info(" {:>8s} {:8d} {:8d} {:8.1f} {:8.3f}", - "Overall", ov.total_observations, ov.unique_reflections, - ov.mean_i_over_sigma, ov.cc_half*100.0); + double completeness = ov.possible_unique_reflections > 0 + ? static_cast(ov.unique_reflections) / ov.possible_unique_reflections * 100.0 : 0.0; + + logger.Info(" {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s}", "", "", "", "", "", "", ""); + logger.Info(" {:>8s} {:8d} {:8d} {:8d} {:7.1f}% {:8.1f} {:7.1f}%", + "Overall", + ov.total_observations, + ov.unique_reflections, + ov.possible_unique_reflections, + completeness, + ov.mean_i_over_sigma, + ov.cc_half*100.0); } logger.Info(""); } diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index e5b26cdd..8ccb2dba 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -10,11 +10,15 @@ #include "../../common/Reflection.h" struct MergeStatisticsShell { + float d_min = 0.0f; float d_max = 0.0f; float mean_one_over_d2 = 0; + int total_observations = 0; int unique_reflections = 0; + int possible_unique_reflections = 0; + double mean_i_over_sigma = 0.0; double cc_half = 0.0f; @@ -47,4 +51,5 @@ std::vector MergeAll(const DiffractionExperiment &x, MergeStatistics MergeStats(const DiffractionExperiment &x, const std::vector &merged, const std::vector > &reflections, + const UnitCell &cell, const std::vector &merge_mask = {}); diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index df83107e..b6fa3a16 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -662,7 +662,7 @@ int main(int argc, char **argv) { logger.Info("Rejected {} images for merging due to low CC with reference", rejected_cc); auto merged_reflections = MergeAll(experiment, indexer.GetReflections(), merging_mask_uc); - auto merged_statistics = MergeStats(experiment, merged_reflections, indexer.GetReflections(), merging_mask_uc); + auto merged_statistics = MergeStats(experiment, merged_reflections, indexer.GetReflections(), *consensus_cell, merging_mask_uc); auto merge_end = std::chrono::steady_clock::now(); double merge_time = std::chrono::duration(merge_end - merge_start).count(); diff --git a/tools/jfjoch_scale.cpp b/tools/jfjoch_scale.cpp index 3d08c322..11f6df07 100644 --- a/tools/jfjoch_scale.cpp +++ b/tools/jfjoch_scale.cpp @@ -289,7 +289,7 @@ int main(int argc, char **argv) { logger.Info("Rejected {} images due to low CC with reference", rejected); auto merge_result = MergeAll(experiment, reflections, merge_mask); - auto merge_stats = MergeStats(experiment, merge_result, reflections, merge_mask); + auto merge_stats = MergeStats(experiment, merge_result, reflections, experiment.GetUnitCell().value(), merge_mask); auto merge_end = std::chrono::steady_clock::now(); double merge_time = std::chrono::duration(merge_end - merge_start).count();