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) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011n8riB6X59oRjkrSHzNPAU

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;

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<OrderKey> order_keys;
+    std::vector<Reflection> 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

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);
 }
 

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);
 }