From abb61e27081e1972fdfecf45cc4c8d20dcc0af4d Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 26 Aug 2026 14:11:56 +0200 Subject: [PATCH] Unmerged export: carry the sort key with the part here too SumRockingEvents assembles rocking events for the unmerged MTZ through the same indirect comparator the anisotropy diagnostic had, on a larger array - it has neither the scaled-image test nor the usable() filter in front of it. It is latent, since --export-unmerged is off by default, but it is the same defect and it is a mechanical fix. Same treatment: the key travels with the part, plus an exact reserve. p_unmerged.mtz, p_unmerged_partials.mtz and p.hkl are md5-identical across eight runs on two crystals. With the export on it is worth about 0.9 s on a seventeen-million-partial run. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01CHMmeM1d489zvNFT7ZMN2P --- image_analysis/WriteReflections.cpp | 34 ++++++++++++++++++++--------- 1 file changed, 24 insertions(+), 10 deletions(-) diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index 6e701d74..df82fc39 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -617,28 +617,42 @@ std::vector SumRockingEvents(const std::vector & double min_partiality) { constexpr float MAX_FRAME_GAP = 2.0f; // == RotationScaleMerge's: what makes one rocking event - std::vector parts; + // The sort key travels with the part instead of being read back through the pointer, the way the + // merge's own ingest sort carries it (RotationScaleMerge's SortKey): there are millions of parts + // and an indirect compare is a cache miss on every one of them. The keys are the same values in + // the same order, so introsort makes the same comparisons and the same swaps and leaves the same + // order - which matters, because two parts can genuinely share (h,k,l) and image_number and the + // event sums below are floating point. + struct Part { + int32_t h, k, l; + float image_number; + const Reflection *r; + }; + std::vector parts; + size_t n_parts = 0; + for (const auto &outcome : outcomes) + n_parts += outcome.reflections.size(); + parts.reserve(n_parts); for (const auto &outcome : outcomes) for (const auto &r : outcome.reflections) - parts.push_back(&r); - std::sort(parts.begin(), parts.end(), [](const Reflection *a, const Reflection *b) { - return std::tie(a->h, a->k, a->l, a->image_number) - < std::tie(b->h, b->k, b->l, b->image_number); + parts.push_back({r.h, r.k, r.l, r.image_number, &r}); + std::sort(parts.begin(), parts.end(), [](const Part &a, const Part &b) { + return std::tie(a.h, a.k, a.l, a.image_number) < std::tie(b.h, b.k, b.l, b.image_number); }); std::vector fulls; for (size_t i = 0; i < parts.size(); ) { size_t j = i + 1; - while (j < parts.size() && parts[j]->h == parts[i]->h && parts[j]->k == parts[i]->k - && parts[j]->l == parts[i]->l - && parts[j]->image_number - parts[j - 1]->image_number <= MAX_FRAME_GAP) + while (j < parts.size() && parts[j].h == parts[i].h && parts[j].k == parts[i].k + && parts[j].l == parts[i].l + && parts[j].image_number - parts[j - 1].image_number <= MAX_FRAME_GAP) ++j; double sum_p = 0.0, sum_I = 0.0, sum_var = 0.0, sum_var_bkg = 0.0; double p_rlp = 0.0, p_frame = 0.0, p_x = 0.0, p_y = 0.0, p_delta_phi = 0.0, p_zeta = 0.0, p_bkg = 0.0; for (size_t m = i; m < j; ++m) { - const Reflection &r = *parts[m]; + const Reflection &r = *parts[m].r; const double p = r.partiality; sum_p += p; sum_I += static_cast(r.I) * r.rlp; @@ -652,7 +666,7 @@ std::vector SumRockingEvents(const std::vector & p_zeta += p * r.zeta; p_bkg += p * r.bkg; } - Reflection full = *parts[i]; + Reflection full = *parts[i].r; i = j; if (sum_p < min_partiality) continue;