mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 72 Wiki Activity

Reduce the anomalous split once per ASU group, not once per observation

Browse Source 
ComputeAsuGroups states the rule for itself - "one ASU reduction per distinct
raw hkl (not per observation)" - and the anomalous split then did a gemmi ASU
reduction and an unordered_map lookup for every one of the millions of fulls.

Both things it wants are properties of the observation's ASU GROUP rather than
of the observation: group_h/k/l is the group's SIGNED representative, so the
same reduction applied to it returns the Friedel-merged key and the hand
together. Reduce once per group into a dense accumulator indexed from there.
The hand only follows the group when the merge distinguishes the hands; a
Friedel-merged run holds both in one group and still has to ask per observation.

SigAno and the merged statistics are unchanged (2.96 over 53303 acentric pairs,
merge table byte-identical).

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-08-24 13:11:24 +02:00
co-authored by Claude Opus 5
parent 29eda90190
commit cf4fbcf96a
1 changed files with 30 additions and 7 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/RotationScaleMerge.cpp
+30 -7
View File
@@ -3146,22 +3146,45 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
const HKLKeyGenerator anom_keygen(/*merge_friedel=*/false, sg_num);
const gemmi::GroupOps gops = gemmi::find_spacegroup_by_number(sg_num)->operations();
struct AnomAcc { double swI[2] = {}; double sw[2] = {}; int32_t h = 0, k = 0, l = 0; float d = NAN; }; // [hand] 0=I(+) 1=I(-)
std::unordered_map<uint64_t, AnomAcc> anom;
anom.reserve(result.merged.size() * 2 + 1);
// The Friedel-merged key an observation lands on, and which hand it is, are properties of its
// ASU GROUP, not of the observation - group_h/k/l is that group's SIGNED representative, so the
// same reduction applied to it returns both. Reduce once per group (a few hundred thousand)
// instead of once per observation (millions), and index a dense accumulator from there, which
// is the rule ComputeAsuGroups states for itself. The hand only follows the group when the
// merge distinguishes the hands; a Friedel-merged run holds both in one group and must ask.
std::vector<int32_t> pair_of_group(n_groups, -1);
std::vector<uint8_t> hand_of_group(merge_friedel ? 0 : n_groups, 0);
std::vector<HKLKey> pair_hkl;
{
std::unordered_map<uint64_t, int32_t> pair_id;
pair_id.reserve(static_cast<size_t>(n_groups) + 1);
for (int g = 0; g < n_groups; ++g) {
const HKLKey gk = anom_keygen(group_h[g], group_k[g], group_l[g]);
if (!merge_friedel) hand_of_group[g] = gk.plus ? 0 : 1;
const uint64_t key = HKLKey{gk.h, gk.k, gk.l, true}.pack();
const auto [it, fresh] = pair_id.emplace(key, static_cast<int32_t>(pair_hkl.size()));
if (fresh) pair_hkl.push_back(HKLKey{gk.h, gk.k, gk.l, true});
pair_of_group[g] = it->second;
}
}
std::vector<AnomAcc> anom(pair_hkl.size());
for (const auto &o : fulls) {
if (!usable_merge(o)) continue;
if (rejected_obs[&o - fulls.data()]) continue; // outlier-rejected in the merge
const HKLKey ak = anom_keygen(o.h, o.k, o.l);
const int hand = ak.plus ? 0 : 1;
const int hand = merge_friedel ? (anom_keygen(o.h, o.k, o.l).plus ? 0 : 1)
: hand_of_group[o.group];
const float I_corr = o.I * o.corr;
const float sigma_corr = corrected_sigma(o, I_corr, o.sigma * o.corr);
if (!(sigma_corr > 0.0f) || !std::isfinite(sigma_corr)) continue;
const double w = 1.0 / (static_cast<double>(sigma_corr) * sigma_corr);
AnomAcc &a = anom[HKLKey{ak.h, ak.k, ak.l, true}.pack()];
a.h = ak.h; a.k = ak.k; a.l = ak.l; a.d = o.d;
AnomAcc &a = anom[pair_of_group[o.group]];
a.d = o.d;
a.swI[hand] += w * static_cast<double>(I_corr); a.sw[hand] += w;
}
for (const auto &[fkey, a] : anom) {
for (size_t p = 0; p < anom.size(); ++p) {
AnomAcc a = anom[p];
a.h = pair_hkl[p].h; a.k = pair_hkl[p].k; a.l = pair_hkl[p].l;
const uint64_t fkey = pair_hkl[p].pack();
// Centrics have I(+)=I(-) by symmetry; leave them without an anomalous split.
if (gops.is_reflection_centric(gemmi::Op::Miller{a.h, a.k, a.l})) continue;
AnomExport ex;