From 81b81ec68364065ffa64db7d17d8e8b7388dd849 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 15 Oct 2026 14:52:22 +0200 Subject: [PATCH] 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 --- .../BraggIntegrationEngineCPU.cpp | 43 ++++++++++++++----- .../BraggIntegrationEngineCPU.h | 5 ++- 2 files changed, 36 insertions(+), 12 deletions(-) diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index aa574ae17..46ba4fefe 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/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 FrameOrder(const std::vector &predicted, size_t npredicted) { + std::vector> keyed(npredicted); + for (size_t i = 0; i < npredicted; ++i) { + const auto band = static_cast(std::max(0.0f, predicted[i].predicted_y)) / 16; + const auto col = static_cast(std::max(0.0f, predicted[i].predicted_x)); + keyed[i] = {(band << 32) | col, static_cast(i)}; + } + std::sort(keyed.begin(), keyed.end()); + std::vector 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(xpixel), static_cast(ypixel), BRAGG_OWNER_NONE), static_mask(mask.GetPackedMask()) {} -void BraggIntegrationEngineCPU::MarkNeighbours(const std::vector &predicted, size_t npredicted) { +void BraggIntegrationEngineCPU::MarkNeighbours(const std::vector &predicted, + const std::vector &order) { const int W = static_cast(xpixel), H = static_cast(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 &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(xpixel), H = static_cast(ypixel); counts.predicted += npredicted; - MarkNeighbours(predicted, npredicted); - for (size_t i = 0; i < npredicted; ++i) { + const std::vector 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(std::floor(r.predicted_x - st.ex_out - 1.0f))); @@ -133,7 +153,9 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, const int W = static_cast(xpixel), H = static_cast(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 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 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(std::floor(r.predicted_x - claim))); const int x1 = std::min(W - 1, static_cast(std::ceil(r.predicted_x + claim))); @@ -190,13 +212,14 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, std::max(0, static_cast(std::floor(r.predicted_y - st.ey_out - 1.0f))), std::min(H - 1, static_cast(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(y) * W + b.x0), @@ -766,7 +789,7 @@ std::vector BraggIntegrationEngineCPU::ProfilePasses(const Sampler & if (!FitOnDevice(setups, shell_profiles, fits)) { std::vector cells_buf; std::vector 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); } diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h index e72497831..7ad908cf5 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h @@ -56,8 +56,9 @@ class BraggIntegrationEngineCPU : public BraggIntegrationEngine { TiledFrame refl_mask; TiledFrame owner; - // Marks every predicted reflection's r2 region in refl_mask (MASK_FLUX / MASK_TAIL). - void MarkNeighbours(const std::vector &predicted, size_t npredicted); + // Marks every predicted reflection's r2 region in refl_mask (MASK_FLUX / MASK_TAIL), visiting them in + // `order` (any order gives the same mask). + void MarkNeighbours(const std::vector &predicted, const std::vector &order); template std::vector RunImpl(const Sampler &img, const std::vector &predicted,