From 484a0e162ab5f7c04ccf8792c0974577cd82f8eb Mon Sep 17 00:00:00 2001 From: jungfrau Date: Sun, 23 Aug 2026 10:55:37 -0400 Subject: [PATCH] Give the predicted reflections an order of their own The GPU predictors claim their output slot with atomicAdd(counter, 1), so a reflection's position in the array is whatever order the blocks happened to finish in. That position is not private to the predictor. BraggOwnerKey packs it into the owner map as the tie-break between two centres equidistant from a shared pixel - the map's atomicMin is order-independent, but the number it compares is not - and the ingest and post-refine bucket sorts, whose comparators are deliberately not total, resolve their ties by the order they are handed. So two runs of the same binary on the same images integrated a different set of reflections. Measured on a large-cell rotation dataset: 63301112 observations against 63301139, and 89% of the merged intensities differing by more than 1% of themselves, median 1.8%. Single-threaded as well as at -N 48, which is what ruled out thread ordering and pointed here. Order the downloaded list by (h, k, l, delta_phi) before TruncateToOutput, whose own pick is then reproducible as well. hkl is a property of the reflection rather than of the schedule, and delta_phi separates the two rocking solutions one hkl can have. The CPU predictors already emit in hkl order, so the two paths now agree on it. Sorting a 20-byte key and gathering once, rather than sorting the 88-byte reflections in place, keeps this off the clock: on a crystal predicting some 35000 reflections a frame the run measures 52.2 s against 52.3 s before. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_011n8riB6X59oRjkrSHzNPAU --- .../bragg_prediction/BraggPrediction.cpp | 18 ++++++++++++++++++ .../bragg_prediction/BraggPrediction.h | 19 +++++++++++++++++++ .../bragg_prediction/BraggPredictionGPU.cu | 1 + .../bragg_prediction/BraggPredictionRotGPU.cu | 1 + 4 files changed, 39 insertions(+) diff --git a/image_analysis/bragg_prediction/BraggPrediction.cpp b/image_analysis/bragg_prediction/BraggPrediction.cpp index 87e042bc..7f939a1d 100644 --- a/image_analysis/bragg_prediction/BraggPrediction.cpp +++ b/image_analysis/bragg_prediction/BraggPrediction.cpp @@ -13,6 +13,24 @@ void BraggPrediction::GrowCapacity(int count) { max_reflections = count; } +void BraggPrediction::OrderOutput(int count) { + order_keys.resize(count); + for (int i = 0; i < count; i++) { + const Reflection &r = reflections[i]; + order_keys[i] = {r.h, r.k, r.l, r.delta_phi_deg, i}; + } + std::sort(order_keys.begin(), order_keys.end(), + [](const OrderKey &a, const OrderKey &b) { + if (a.h != b.h) return a.h < b.h; + if (a.k != b.k) return a.k < b.k; + if (a.l != b.l) return a.l < b.l; + return a.delta_phi_deg < b.delta_phi_deg; + }); + order_scratch.resize(count); + for (int i = 0; i < count; i++) order_scratch[i] = reflections[order_keys[i].index]; + std::copy(order_scratch.begin(), order_scratch.end(), reflections.begin()); +} + int BraggPrediction::TruncateToOutput(int count) { if (count <= output_limit) return count; diff --git a/image_analysis/bragg_prediction/BraggPrediction.h b/image_analysis/bragg_prediction/BraggPrediction.h index b70344f4..4d14a44b 100644 --- a/image_analysis/bragg_prediction/BraggPrediction.h +++ b/image_analysis/bragg_prediction/BraggPrediction.h @@ -58,6 +58,25 @@ protected: // excitation error on the still path - with hkl breaking what is left. Returns the kept count. Below // the cap it is a no-op. Call at the end of every Calc override. int TruncateToOutput(int count); + + // Put the first `count` predicted reflections in an order that depends only on the reflections + // themselves. The GPU predictors append at an atomic counter, so a reflection's POSITION in the + // array is decided by the order the blocks happened to finish - and that position is not private + // to the predictor: BraggOwnerKey packs it into the owner map as the tie-break between two centres + // equidistant from a shared pixel, and every downstream sort that is not a total order (the + // ingest and post-refine bucket sorts) resolves its ties by the order it receives. Two runs of the + // same binary on the same image therefore integrated a different set of reflections. hkl is a + // property of the reflection; delta_phi separates the two rocking solutions one hkl can have. + // Call before TruncateToOutput, whose own pick is then reproducible as well. + void OrderOutput(int count); + + // Scratch for OrderOutput. A Reflection is ~88 bytes and a large cell predicts tens of thousands + // of them per frame, so sorting the structs themselves moves several megabytes an image; sorting + // a 20-byte key and gathering once is the same order for a third of the traffic. Members rather + // than locals so the two allocations happen once per engine, not once per image. + struct OrderKey { int32_t h, k, l; float delta_phi_deg; int32_t index; }; + std::vector order_keys; + std::vector order_scratch; public: // The prediction buffer holds up to kPredictionCapacity reflections so a strong lattice does not // overflow it. Calc returns at most output_limit, the number that flows downstream and is serialized - kept low so the diff --git a/image_analysis/bragg_prediction/BraggPredictionGPU.cu b/image_analysis/bragg_prediction/BraggPredictionGPU.cu index 87a0b397..91a45ac7 100644 --- a/image_analysis/bragg_prediction/BraggPredictionGPU.cu +++ b/image_analysis/bragg_prediction/BraggPredictionGPU.cu @@ -232,6 +232,7 @@ int BraggPredictionGPU::Calc(const DiffractionExperiment &experiment, cudaMemcpyAsync(reflections.data(), d_out, sizeof(Reflection) * count, cudaMemcpyDeviceToHost, stream); cudaStreamSynchronize(stream); + OrderOutput(count); return TruncateToOutput(count); } diff --git a/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu b/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu index ecf71266..c67e4413 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu +++ b/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu @@ -328,6 +328,7 @@ int BraggPredictionRotGPU::Calc(const DiffractionExperiment &experiment, cudaMemcpyAsync(reflections.data(), d_out, sizeof(Reflection) * count, cudaMemcpyDeviceToHost, stream); cudaStreamSynchronize(stream); + OrderOutput(count); return TruncateToOutput(count); }