From 674ac4afa24cd1f6450c6cc96f081361185ac6b2 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sun, 17 May 2026 10:27:19 +0200 Subject: [PATCH] Merge: For anomalous signal, treat Friedel pairs separate, but add half-datasets --- common/Reflection.h | 4 +- image_analysis/scale_merge/HKLKey.cpp | 6 +- image_analysis/scale_merge/HKLKey.h | 2 +- image_analysis/scale_merge/Merge.cpp | 126 ++++++++------------------ tests/MergeScaleTest.cpp | 15 ++- 5 files changed, 60 insertions(+), 93 deletions(-) diff --git a/common/Reflection.h b/common/Reflection.h index 96d053dd..12f7b23d 100644 --- a/common/Reflection.h +++ b/common/Reflection.h @@ -36,8 +36,8 @@ struct MergedReflection { int32_t l = 0; float I = NAN; float sigma = NAN; - float I_anom[2] = {NAN, NAN}; - float sigma_anom[2] = {NAN, NAN}; + float I_half[2] = {NAN, NAN}; + float sigma_half[2] = {NAN, NAN}; float d = 0.0; }; diff --git a/image_analysis/scale_merge/HKLKey.cpp b/image_analysis/scale_merge/HKLKey.cpp index 0f72d0ce..3bdbfe61 100644 --- a/image_analysis/scale_merge/HKLKey.cpp +++ b/image_analysis/scale_merge/HKLKey.cpp @@ -6,7 +6,7 @@ #include "HKLKey.h" #include "gemmi/symmetry.hpp" -uint64_t HKLKey::pack_no_anom() const { +uint64_t HKLKey::pack() const { constexpr int bits = 21; constexpr int bias = 1 << (bits - 1); // 1,048,576 constexpr int max_value = bias - 1; @@ -17,13 +17,13 @@ uint64_t HKLKey::pack_no_anom() const { k < min_value || k > max_value || l < min_value || l > max_value) { throw std::out_of_range("HKL index outside packable range"); - } + } const std::uint64_t hh = static_cast(h + bias) & mask; const std::uint64_t kk = static_cast(k + bias) & mask; const std::uint64_t ll = static_cast(l + bias) & mask; - return hh | (kk << bits) | (ll << (2 * bits)); + return (hh << 1) | (kk << (bits + 1)) | (ll << (2 * bits + 1)) | (plus ? 1ULL : 0ULL); } HKLKeyGenerator::HKLKeyGenerator(bool merge_friedel, int32_t space_group_number) diff --git a/image_analysis/scale_merge/HKLKey.h b/image_analysis/scale_merge/HKLKey.h index 74af3aff..51fbab48 100644 --- a/image_analysis/scale_merge/HKLKey.h +++ b/image_analysis/scale_merge/HKLKey.h @@ -24,7 +24,7 @@ struct HKLKey { return h == o.h && k == o.k && l == o.l && plus == o.plus; } - [[nodiscard]] uint64_t pack_no_anom() const; + [[nodiscard]] uint64_t pack() const; }; class HKLKeyGenerator { diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index c950f25c..51702499 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -34,16 +34,15 @@ std::vector MergeAll(const DiffractionExperiment &x, auto min_partiality = scaling_settings.GetMinPartiality(); struct Accum { - // Keep anomalous + / - together, but separate int32_t h = 0; int32_t k = 0; int32_t l = 0; float d = NAN; double sum_wI = 0.0; double sum_w = 0.0; - double sum_wI_anom[2] = {0.0, 0.0}; - double sum_w_anom[2] = {0.0, 0.0}; - bool present[2] = {false, false}; + double sum_wI_half[2] = {0.0, 0.0}; + double sum_w_half[2] = {0.0, 0.0}; + size_t n_half[2] = {0, 0}; }; std::unordered_map acc; @@ -65,25 +64,27 @@ std::vector MergeAll(const DiffractionExperiment &x, continue; auto hkl = key_generator(r); - auto hkl_key = hkl.pack_no_anom(); + auto hkl_key = hkl.pack(); auto it = acc.find(hkl_key); if (it == acc.end()) it = acc.emplace(hkl_key, Accum{ - .h = hkl.h, - .k = hkl.k, - .l = hkl.l - }).first; - - int solution = hkl.plus ? 0 : 1; + .h = hkl.plus ? hkl.h : -hkl.h, + .k = hkl.plus ? hkl.k : -hkl.k, + .l = hkl.plus ? hkl.l : -hkl.l, + }).first; const float w = 1.0f / (sigma_corr * sigma_corr); - it->second.sum_wI += w * I_corr; - it->second.sum_w += w; - it->second.sum_wI_anom[solution] += w * I_corr; - it->second.sum_w_anom[solution] += w; + const float wI = w * I_corr; + + it->second.sum_wI += wI; + it->second.sum_w += w; + + const int half = ((it->second.n_half[0] + it->second.n_half[1]) % 2 == 0) ? 0 : 1; + it->second.sum_wI_half[half] += wI; + it->second.sum_w_half[half] += w; + it->second.n_half[half]++; - it->second.present[solution] = true; if (!std::isfinite(it->second.d) && std::isfinite(r.d) && r.d > 0.0f) it->second.d = r.d; } @@ -95,25 +96,26 @@ std::vector MergeAll(const DiffractionExperiment &x, if (accum.sum_w <= 0.0) continue; - float I_anom[2] = {NAN, NAN}; - float sigma_anom[2] = {NAN, NAN}; - - for (int i = 0; i < 2; ++i) { - if (accum.present[i] && accum.sum_w_anom[i] > 0.0) { - I_anom[i] = static_cast(accum.sum_wI_anom[i] / accum.sum_w_anom[i]); - sigma_anom[i] = 1.0f / std::sqrt(static_cast(accum.sum_w_anom[i])); - } - } - out.emplace_back(MergedReflection{ + MergedReflection mr{ .h = accum.h, .k = accum.k, .l = accum.l, .I = static_cast(accum.sum_wI / accum.sum_w), .sigma = 1.0f / std::sqrt(static_cast(accum.sum_w)), - .I_anom = {I_anom[0], I_anom[1]}, - .sigma_anom = {sigma_anom[0], sigma_anom[1]}, + .I_half = {NAN, NAN}, + .sigma_half = {NAN, NAN}, .d = accum.d - }); + }; + + if (accum.n_half[0] + accum.n_half[1] > 0 + && accum.sum_w_half[0] > 0.0 + && accum.sum_w_half[1] > 0.0) { + for (int i = 0; i < 2; ++i) { + mr.I_half[i] = static_cast(accum.sum_wI_half[i] / accum.sum_w_half[i]); + mr.sigma_half[i] = 1.0f / std::sqrt(static_cast(accum.sum_w_half[i])); + } + } + out.emplace_back(mr); } return out; } @@ -136,15 +138,6 @@ MergeStatistics MergeStats(const DiffractionExperiment &x, auto scaling_settings = x.GetScalingSettings(); HKLKeyGenerator key_generator(scaling_settings.GetMergeFriedel(), x.GetSpaceGroupNumber().value_or(1)); - struct HalfAccum { - float d = NAN; - double sum_wI[2] = {0.0, 0.0}; - double sum_w[2] = {0.0, 0.0}; - int n[2] = {0, 0}; - }; - - std::unordered_map half_acc; - for (const auto &m: merged) { if (!std::isfinite(m.d) || m.d <= 0.0f) continue; @@ -192,8 +185,16 @@ MergeStatistics MergeStats(const DiffractionExperiment &x, acc[s].unique++; acc[s].sum_i_over_sigma += m.I / m.sigma; ++acc[s].n_i_over_sigma; - } + if (std::isfinite(m.I_half[0]) && std::isfinite(m.I_half[1])) { + acc[s].n_cc_half += 1; + acc[s].sum_x += m.I_half[0]; + acc[s].sum_y += m.I_half[1]; + acc[s].sum_x2 += m.I_half[0] * m.I_half[0]; + acc[s].sum_y2 += m.I_half[1] * m.I_half[1]; + acc[s].sum_xy += m.I_half[0] * m.I_half[1]; + } + } } } @@ -221,56 +222,9 @@ MergeStatistics MergeStats(const DiffractionExperiment &x, const int s = *shell; if (s >= 0 && s < n_shells) acc[s].total_obs++; - - const auto hkl = key_generator(r); - const auto hkl_key = hkl.pack_no_anom(); - - auto it = half_acc.find(hkl_key); - if (it == half_acc.end()) - it = half_acc.emplace(hkl_key, HalfAccum{}).first; - - if (!std::isfinite(it->second.d) && std::isfinite(r.d) && r.d > 0.0f) - it->second.d = r.d; - - const int half = (it->second.n[0] + it->second.n[1]) % 2; - const double w = 1.0 / static_cast(sigma_corr * sigma_corr); - - it->second.sum_wI[half] += w * I_corr; - it->second.sum_w[half] += w; - it->second.n[half]++; } } - for (const auto &[key, ha]: half_acc) { - if (ha.n[0] == 0 || ha.n[1] == 0) - continue; - if (ha.sum_w[0] <= 0.0 || ha.sum_w[1] <= 0.0) - continue; - if (!std::isfinite(ha.d) || ha.d <= 0.0f) - continue; - - const auto shell = shells.GetShell(ha.d); - if (!shell.has_value()) - continue; - - const int s = *shell; - if (s < 0 || s >= n_shells) - continue; - - const double x = ha.sum_wI[0] / ha.sum_w[0]; - const double y = ha.sum_wI[1] / ha.sum_w[1]; - - if (!std::isfinite(x) || !std::isfinite(y)) - continue; - - acc[s].sum_x += x; - acc[s].sum_y += y; - acc[s].sum_x2 += x * x; - acc[s].sum_y2 += y * y; - acc[s].sum_xy += x * y; - acc[s].n_cc_half++; - } - MergeStatistics out; out.shells.resize(n_shells); diff --git a/tests/MergeScaleTest.cpp b/tests/MergeScaleTest.cpp index 534871d7..71d76c38 100644 --- a/tests/MergeScaleTest.cpp +++ b/tests/MergeScaleTest.cpp @@ -48,7 +48,6 @@ TEST_CASE("HKLKey_SG96_MergeFriedel") { CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(2, 1, 3)); } - TEST_CASE("HKLKey_SG96_NoMergeFriedel") { HKLKeyGenerator hkl_key_gen(false, *gemmi::find_spacegroup_by_number(96)); CHECK(hkl_key_gen(-1, -2, -3) != hkl_key_gen(1,2,3)); @@ -64,3 +63,17 @@ TEST_CASE("HKLKey_SG96_NoMergeFriedel") { CHECK(hkl_key_gen(1,2,3) != hkl_key_gen(-2,-1,3)); CHECK(hkl_key_gen(1,2,3) != hkl_key_gen(2, 1, 3)); } + +TEST_CASE("HKLKey_pack_friedel") { + HKLKeyGenerator hkl_key_gen(false, *gemmi::find_spacegroup_by_number(1)); + CHECK(hkl_key_gen(-1, -2, -3).pack() == hkl_key_gen(1,2,3).pack()); + CHECK(hkl_key_gen(-1,-2,-3).pack() == hkl_key_gen(-1,-2,-3).pack()); + CHECK(hkl_key_gen(-1,-2,-3).pack() != hkl_key_gen(1,-2,-3).pack()); +} + +TEST_CASE("HKLKey_pack_no_friedel") { + HKLKeyGenerator hkl_key_gen(true, *gemmi::find_spacegroup_by_number(1)); + CHECK(hkl_key_gen(-1, -2, -3).pack() == hkl_key_gen(1,2,3).pack()); + CHECK(hkl_key_gen(-1,-2,-3).pack() == hkl_key_gen(-1,-2,-3).pack()); + CHECK(hkl_key_gen(-1,-2,-3).pack() == hkl_key_gen(1,-2,-3).pack()); +}