CPU Bragg integration: reflections walked in detector order, output unchanged

The prediction lists a frame's reflections by Miller index, which scatters consecutive ones over the
whole detector, so pass A's box walk, the reflection-mask and owner-map builds and the CPU profile fit
each waited on memory for every box's image rows and mask tiles (the pixel load alone was 18% of pass
A's samples). These loops now visit the reflections sorted by 16-row band, then column. Each of them
gives the same result in any order - integer sums, min/max, per-reflection outputs indexed by the
reflection - so the output is unchanged; the profile learning, whose double sums do depend on the
order, still walks the original one.

cytochrome C --no-gpu: user cycles 10.71e12 -> 10.20e12 (-4.8%); p.hkl and MTZ data identical.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp

@@ -53,6 +53,24 @@ struct ImageSampler {
     const void *address(size_t i) const { return p + i; }
 };
 
+// The reflections in detector order: by band of 16 rows (a mask tile's height), then by column. The
+// prediction lists them by Miller index, which scatters consecutive reflections over the whole frame, so
+// each box waited on memory for its image rows and mask tiles; in this order neighbours share them in
+// cache. Only for the loops whose results do not depend on the order the reflections are visited in.
+std::vector<uint32_t> FrameOrder(const std::vector<Reflection> &predicted, size_t npredicted) {
+    std::vector<std::pair<uint64_t, uint32_t>> keyed(npredicted);
+    for (size_t i = 0; i < npredicted; ++i) {
+        const auto band = static_cast<uint64_t>(std::max(0.0f, predicted[i].predicted_y)) / 16;
+        const auto col = static_cast<uint64_t>(std::max(0.0f, predicted[i].predicted_x));
+        keyed[i] = {(band << 32) | col, static_cast<uint32_t>(i)};
+    }
+    std::sort(keyed.begin(), keyed.end());
+    std::vector<uint32_t> order(npredicted);
+    for (size_t k = 0; k < npredicted; ++k)
+        order[k] = keyed[k].second;
+    return order;
+}
+
 }   // namespace
 
 BraggIntegrationEngineCPU::BraggIntegrationEngineCPU(const DiffractionExperiment &experiment,
@@ -62,7 +80,8 @@ BraggIntegrationEngineCPU::BraggIntegrationEngineCPU(const DiffractionExperiment
       owner(static_cast<int>(xpixel), static_cast<int>(ypixel), BRAGG_OWNER_NONE),
       static_mask(mask.GetPackedMask()) {}
 
-void BraggIntegrationEngineCPU::MarkNeighbours(const std::vector<Reflection> &predicted, size_t npredicted) {
+void BraggIntegrationEngineCPU::MarkNeighbours(const std::vector<Reflection> &predicted,
+                                               const std::vector<uint32_t> &order) {
     const int W = static_cast<int>(xpixel), H = static_cast<int>(ypixel);
     // --- Reflection mask: mark the r2 signal region of every predicted reflection, so a neighbour's
     //     signal is excluded from this reflection's r2..r3 background ring. MASK_FLUX where the
@@ -72,7 +91,7 @@ void BraggIntegrationEngineCPU::MarkNeighbours(const std::vector<Reflection> &pr
     //     neighbour's OWN frame - an elongated neighbour whose streak is still marked as a disk would
     //     leak its tails into this reflection's ring. ---
     refl_mask.Clear();
-    for (size_t i = 0; i < npredicted; ++i) {
+    for (const uint32_t i : order) {
         const auto &r = predicted[i];
         const uint8_t level = r.partiality >= NEIGHBOUR_MASK_MIN_PARTIALITY ? MASK_FLUX : MASK_TAIL;
         const BraggStencil st = MakeBraggStencil(r.predicted_x, r.predicted_y, stencil);
@@ -95,8 +114,9 @@ void BraggIntegrationEngineCPU::CountNeighbourCrowding(const std::vector<Reflect
                                                        size_t npredicted) {
     const int W = static_cast<int>(xpixel), H = static_cast<int>(ypixel);
     counts.predicted += npredicted;
-    MarkNeighbours(predicted, npredicted);
-    for (size_t i = 0; i < npredicted; ++i) {
+    const std::vector<uint32_t> order = FrameOrder(predicted, npredicted);
+    MarkNeighbours(predicted, order);
+    for (const uint32_t i : order) {
         const auto &r = predicted[i];
         const BraggStencil st = MakeBraggStencil(r.predicted_x, r.predicted_y, stencil);
         const int x0 = std::max(0, static_cast<int>(std::floor(r.predicted_x - st.ex_out - 1.0f)));
@@ -133,7 +153,9 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
     const int W = static_cast<int>(xpixel), H = static_cast<int>(ypixel);
     const bool do_clip = bkg_clip_nsigma > 0.0f && mode != IntegratorMode::BoxSum;
 
-    MarkNeighbours(predicted, npredicted);
+    // Every loop below whose results do not depend on the order walks the reflections in this one.
+    const std::vector<uint32_t> order = FrameOrder(predicted, npredicted);
+    MarkNeighbours(predicted, order);
     owner.Clear();
 
     // --- Signal-region ownership: a pixel inside two reflections' signal regions belongs to the
@@ -145,7 +167,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
     const bool use_owner = overlap != OverlapMode::Off;
     if (use_owner) {
         const float claim_sq = claim * claim;
-        for (size_t i = 0; i < npredicted; ++i) {
+        for (const uint32_t i : order) {
             const auto &r = predicted[i];
             const int x0 = std::max(0, static_cast<int>(std::floor(r.predicted_x - claim)));
             const int x1 = std::min(W - 1, static_cast<int>(std::ceil(r.predicted_x + claim)));
@@ -190,13 +212,14 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
                     std::max(0, static_cast<int>(std::floor(r.predicted_y - st.ey_out - 1.0f))),
                     std::min(H - 1, static_cast<int>(std::ceil(r.predicted_y + st.ey_out + 1.0f)))};
     };
-    for (size_t i = 0; i < npredicted; ++i) {
+    for (size_t k = 0; k < npredicted; ++k) {
+        const size_t i = order[k];
         // The next reflection's image rows, asked for while this one is worked on: the reflections
         // are spread over a frame far larger than the caches, so each box would otherwise wait on
         // memory row by row. (Written out here: GCC drops a call to a function that does nothing but
         // prefetch.)
-        if (i + 1 < npredicted) {
-            const Reflection &n = predicted[i + 1];
+        if (k + 1 < npredicted) {
+            const Reflection &n = predicted[order[k + 1]];
             const Box b = box_of(n, MakeBraggStencil(n.predicted_x, n.predicted_y, stencil));
             for (int y = b.y0; y <= b.y1 && b.x0 <= b.x1; ++y)
                 prefetch_row(img.address(static_cast<size_t>(y) * W + b.x0),
@@ -766,7 +789,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::ProfilePasses(const Sampler &
     if (!FitOnDevice(setups, shell_profiles, fits)) {
         std::vector<float> cells_buf;
         std::vector<FitPixel> fit_pixels;
-        for (size_t i = 0; i < npredicted; ++i)
+        for (const uint32_t i : FrameOrder(predicted, npredicted))
             if (setups[i].fit)
                 fits[i] = FitProfile(img, owner_map, setups[i], i, shell_profiles, cells_buf, fit_pixels);
     }