Merge: For anomalous signal, treat Friedel pairs separate, but add half-datasets

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;
 };
 

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<std::uint64_t>(h + bias) & mask;
     const std::uint64_t kk = static_cast<std::uint64_t>(k + bias) & mask;
     const std::uint64_t ll = static_cast<std::uint64_t>(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)

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 {

image_analysis/scale_merge/Merge.cpp

@@ -34,16 +34,15 @@ std::vector<MergedReflection> 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<uint64_t, Accum> acc;
@@ -65,25 +64,27 @@ std::vector<MergedReflection> 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;
+            const float wI = w * I_corr;
+
+            it->second.sum_wI += wI;
             it->second.sum_w += w;
-            it->second.sum_wI_anom[solution] += w * I_corr;
-            it->second.sum_w_anom[solution] += 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<MergedReflection> 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<float>(accum.sum_wI_anom[i] / accum.sum_w_anom[i]);
-                sigma_anom[i] = 1.0f / std::sqrt(static_cast<float>(accum.sum_w_anom[i]));
-            }
-        }
-        out.emplace_back(MergedReflection{
+        MergedReflection mr{
             .h = accum.h,
             .k = accum.k,
             .l = accum.l,
             .I = static_cast<float>(accum.sum_wI / accum.sum_w),
             .sigma = 1.0f / std::sqrt(static_cast<float>(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<float>(accum.sum_wI_half[i] / accum.sum_w_half[i]);
+                mr.sigma_half[i] = 1.0f / std::sqrt(static_cast<float>(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<uint64_t, HalfAccum> 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<double>(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);
 

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());
+}