From 2f54a1189deaed36a1065bfa707c080d4c5330ca Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 24 Aug 2026 18:50:50 +0200 Subject: [PATCH] Take the error-model split and the ASU grouping off one thread Three regions of the merge tail, measured with instrumented timers and confirmed against a cycle profile. On a tail-heavy dataset the scale and merge tail is 70% of the run's wall clock at six of thirty-two logical cores busy, with the GPU idle 88% of the time, so this is where the CPU headroom is. fit_error_model ran a serial four-level nth_element cascade over the whole sample pool, twelve times per dataset. The two halves either side of a partition are disjoint and their contents are already fixed by the parent's nth_element, so the recursion can descend both at once; it now does while a range is worth a thread. The bins are unchanged. ComputeAsuGroups sorted indices with an indirect comparator, taking a cache miss per comparison into an array far larger than the last-level cache. It now sorts packed key-and-run pairs. Tie order does not matter because the packed key encodes h, k, l and the hand exactly, so every run in a tie reduces to the same reflection. The per-thread histogram prefix walked thirty-two separate histograms column-wise on one thread. It becomes a parallel per-group total, one sequential scan over two flat arrays, and a parallel hand-out of the bases - the same sums in the same order. Faster on 21 of 23 matched pairs in an alternating A/B, and on 15 of 15 in the quieter of the two sessions: 0.6% to 2.3% of whole-run wall clock depending on the dataset, around 1.8% in aggregate, and 3 to 4% of the time spent outside the image loop. The reflection files are byte-identical on every dataset tested. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_016NNnL26LAvruQ9eLUUWvrJ --- docs/CHANGELOG.md | 1 + .../scale_merge/RotationScaleMerge.cpp | 70 +++++++++++++------ 2 files changed, 51 insertions(+), 20 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index c5bf88fe..026525a3 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -18,6 +18,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as * HDF5 and image stream: `mirror_y` records whether the assembled image is mirrored in Y relative to the detector's raw readout. * rugnux: an image integrated in pyFAI through the `.poni` file written by `--mode calibration` now comes out with the correct azimuth. Radial integration is unchanged. * rugnux: the `.poni` file declares pyFAI's `orientation`, which needs pyFAI 2024.01 or newer. +* rugnux: scaling and merging are faster, with identical output. * The per-image resolution estimate now predicts the resolution the merged data reach, rather than reporting the highest-resolution spot found; rugnux reports the run's value as `SPOT_RESOLUTION_ESTIMATE` in its report. * rugnux: fixing the space group with `-S` no longer prevents the lattice from being found; the group is applied to scaling and merging rather than to the indexing search. * rugnux: the detector geometry is also logged in XDS's convention (`ORGX`/`ORGY`, detector axis vectors, rotation axis), so it can be compared with an XDS refinement. diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 53c7f5a8..97df6bbb 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -927,24 +927,30 @@ int RotationScaleMerge::ComputeAsuGroups(const HKLKeyGenerator &keygen) { eligible[r] = 1; } }); - std::vector idx; - idx.reserve(n_run); + // The key travels with the run it belongs to, so the sort and the run-detection below read the + // key they are standing on instead of chasing it back through `key` - an indirect compare is a + // cache miss per comparison, and there are a few tens of millions of them. Which of two runs + // with the SAME key ends up first is left to the sort, and cannot matter: the packed key holds + // h, k, l and the hand exactly, so every run in a tie reduces to the same ASU reflection. + struct RunKey { uint64_t key; int32_t run; }; + std::vector sorted; + sorted.reserve(n_run); for (int r = 0; r < n_run; ++r) - if (eligible[r]) idx.push_back(r); - std::sort(idx.begin(), idx.end(), [&](int32_t a, int32_t b) { return key[a] < key[b]; }); + if (eligible[r]) sorted.push_back({key[r], r}); + std::sort(sorted.begin(), sorted.end(), [](const RunKey &a, const RunKey &b) { return a.key < b.key; }); group_h.clear(); group_k.clear(); group_l.clear(); int n_groups = 0; - for (size_t j = 0; j < idx.size(); ++j) { - if (j == 0 || key[idx[j]] != key[idx[j - 1]]) { - const int r = idx[j]; + for (size_t j = 0; j < sorted.size(); ++j) { + if (j == 0 || sorted[j].key != sorted[j - 1].key) { + const int r = sorted[j].run; const auto hkl = keygen(rawrun_h[r], rawrun_k[r], rawrun_l[r]); group_h.push_back(hkl.plus ? hkl.h : -hkl.h); group_k.push_back(hkl.plus ? hkl.k : -hkl.k); group_l.push_back(hkl.plus ? hkl.l : -hkl.l); ++n_groups; } - rawrun_group[idx[j]] = n_groups - 1; + rawrun_group[sorted[j].run] = n_groups - 1; } // Stamp the ASU-group id per obs from its raw hkl + the precomputed finiteness. For the GPU we build a @@ -1001,15 +1007,27 @@ int RotationScaleMerge::ComputeAsuGroups(const HKLKeyGenerator &keygen) { for (auto &f : f1) f.get(); // CSR starts + convert hist[t][g] into chunk t's write base for group g (exclusive prefix over t). + // Only the prefix ACROSS groups is sequential; the sum over the chunks inside one group, and + // handing that group's chunks their bases, are per-group and go on all threads. The single loop + // this replaces walked nt separate histograms one column at a time on one thread, which is nt + // cache misses per group. Same sums in the same order, so the same numbers. std::vector gstart(n_groups), gcount(n_groups); + const size_t hist_work = static_cast(n_groups) * static_cast(nt); + ParallelChunks(n_groups, ThreadsForWork(hist_work, nthreads), [&](int glo, int ghi) { + for (int g = glo; g < ghi; ++g) { + int c = 0; + for (int t = 0; t < nt; ++t) c += hist[t][g]; + gcount[g] = c; + } + }); int acc = 0; - for (int g = 0; g < n_groups; ++g) { - int base = acc; - gstart[g] = acc; - for (int t = 0; t < nt; ++t) { const int c = hist[t][g]; hist[t][g] = base; base += c; } - gcount[g] = base - acc; - acc = base; - } + for (int g = 0; g < n_groups; ++g) { gstart[g] = acc; acc += gcount[g]; } + ParallelChunks(n_groups, ThreadsForWork(hist_work, nthreads), [&](int glo, int ghi) { + for (int g = glo; g < ghi; ++g) { + int base = gstart[g]; + for (int t = 0; t < nt; ++t) { const int c = hist[t][g]; hist[t][g] = base; base += c; } + } + }); // Pass 2 (parallel): each chunk fills its obs into gperm at its per-group base (stable). std::vector gperm(acc); @@ -2594,16 +2612,28 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool if (a.dev2 != b.dev2) return a.dev2 < b.dev2; return a.d < b.d; }; - const std::function split = - [&](size_t lo, size_t hi, int b0, int b1) { + // The two halves a split leaves behind are disjoint ranges whose contents the nth_element + // above them has already fixed, so a child sees the same sequence whether it runs before, + // after or beside its sibling and the bins come out the same either way. Run them side by + // side while the range is still worth a thread; below that they stay on this one. + constexpr size_t SPLIT_MIN_PARALLEL = 1u << 17; + const std::function split = + [&](size_t lo, size_t hi, int b0, int b1, int depth) { if (b0 >= b1) return; const int mid = (b0 + b1) / 2; const size_t k = static_cast(mid) * per; std::nth_element(smp.begin() + lo, smp.begin() + k, smp.begin() + hi, by_I2); - split(lo, k, b0, mid); - split(k, hi, mid + 1, b1); + if (depth > 0 && hi - lo > SPLIT_MIN_PARALLEL) { + auto left = std::async(std::launch::async, + [&, lo, k, b0, mid, depth] { split(lo, k, b0, mid, depth - 1); }); + split(k, hi, mid + 1, b1, depth - 1); + left.get(); + } else { + split(lo, k, b0, mid, 0); + split(k, hi, mid + 1, b1, 0); + } }; - split(0, smp.size(), 1, n_bins); + split(0, smp.size(), 1, n_bins, nthreads > 1 ? 3 : 0); std::vector bs2(n_bins), bI2(n_bins), bd2(n_bins); ParallelChunks(n_bins, ThreadsForWork(smp.size(), nthreads, 8 * 32768), [&](int blo, int bhi) { for (int bin = blo; bin < bhi; ++bin) {