From 46bb3bdbab1e0d548ae0bf72f80ecdf6eecd151d Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 30 Jul 2026 18:22:39 +0200 Subject: [PATCH] Bragg prediction: say so when a frame overflows the prediction buffer Found while chasing a 12% run-to-run spread in the merged reflection count of one crystal. The GPU kernels claim output slots with an atomicAdd and, on overflow, undid the increment with an atomicSub - so the counter saturated at the capacity and the host could not tell a full buffer from an overflowing one. Which reflections survived was then decided by CUDA block scheduling and changed every run. Measured on that dataset: every frame predicts 23000-44000 against a 20000 buffer, and the spread reached the merged output (161591 / 165193 / 166110 / 166479 unique across four runs of the same command). Single-threaded runs diverge too - this is entirely GPU-side. Stop clamping the counter, so the true number predicted reaches the host, and warn once per predictor when it exceeds the buffer. Which reflections are kept is unchanged: making that reproducible means deciding what to keep when a frame predicts more than the pipeline carries, and the obvious answers are worse - the capacity is not the real limit, kPredictionOutput (10000, selected by smallest excitation error) is, and on this crystal both a bigger buffer and a strided selection collapse the merge, because the rotation combine rebuilds fulls from exactly the partials that a smallest-excitation-error cut throws away. Co-Authored-By: Claude Opus 5 (1M context) --- image_analysis/bragg_prediction/BraggPrediction.cpp | 12 ++++++++++++ image_analysis/bragg_prediction/BraggPrediction.h | 13 +++++++++++-- .../bragg_prediction/BraggPredictionGPU.cu | 9 ++++++--- .../bragg_prediction/BraggPredictionRotGPU.cu | 11 +++++++---- 4 files changed, 36 insertions(+), 9 deletions(-) diff --git a/image_analysis/bragg_prediction/BraggPrediction.cpp b/image_analysis/bragg_prediction/BraggPrediction.cpp index ac483c60..5e4460f3 100644 --- a/image_analysis/bragg_prediction/BraggPrediction.cpp +++ b/image_analysis/bragg_prediction/BraggPrediction.cpp @@ -4,9 +4,21 @@ #include #include "../../common/JFJochMath.h" +#include "../../common/Logger.h" #include "BraggPrediction.h" #include "../bragg_integration/SystematicAbsence.h" +void BraggPrediction::ReportOverflow(int predicted) { + if (overflow_reported) + return; + overflow_reported = true; + Logger("BraggPrediction").Warning( + "A frame predicted {} reflections but the buffer holds {} - the ones kept are whichever the GPU " + "wrote first, so this dataset does NOT process reproducibly. Results will differ between runs. " + "The cell is large enough that the prediction cap no longer fits it.", + predicted, kPredictionCapacity); +} + int BraggPrediction::TruncateToOutput(int count) { if (count <= kPredictionOutput) return count; diff --git a/image_analysis/bragg_prediction/BraggPrediction.h b/image_analysis/bragg_prediction/BraggPrediction.h index 60400b00..479dcf63 100644 --- a/image_analysis/bragg_prediction/BraggPrediction.h +++ b/image_analysis/bragg_prediction/BraggPrediction.h @@ -30,14 +30,23 @@ protected: const int max_reflections; std::vector reflections; + // A frame that predicts more than the buffer holds keeps an ARBITRARY subset of them: the GPU + // kernels claim slots with an atomicAdd, so which ones survive depends on block scheduling and + // changes from run to run. Say so, once per predictor, rather than let it pass silently - it is + // not a small effect (measured: a 3% run-to-run spread in the number of merged reflections, and + // every frame of that dataset overflowed). + void ReportOverflow(int predicted); + bool overflow_reported = false; + // Deterministically cap Calc's output at kPredictionOutput: if more were predicted, keep the ones // closest to the Ewald sphere (smallest excitation error), ties broken by hkl. Returns the kept // count. Below the cap it is a no-op. Call at the end of every Calc override. int TruncateToOutput(int count); public: // The prediction buffer holds up to kPredictionCapacity reflections so a strong lattice does not - // overflow it (the GPU kernels then fill it in a non-deterministic atomic order). Calc returns at - // most kPredictionOutput, the number that flows downstream and is serialized - kept low so the + // overflow it (the GPU kernels then fill it in a non-deterministic atomic order - ReportOverflow + // warns when that happens; a large unit cell, ~2.8e6 A^3, reaches ~40000 per frame and overflows on + // every one). Calc returns at most kPredictionOutput, the number that flows downstream and is serialized - kept low so the // per-image reflection list stays within the frame transport headroom. static constexpr int kPredictionCapacity = 20000; static constexpr int kPredictionOutput = 10000; diff --git a/image_analysis/bragg_prediction/BraggPredictionGPU.cu b/image_analysis/bragg_prediction/BraggPredictionGPU.cu index 6a606182..43ebefbf 100644 --- a/image_analysis/bragg_prediction/BraggPredictionGPU.cu +++ b/image_analysis/bragg_prediction/BraggPredictionGPU.cu @@ -148,9 +148,9 @@ namespace { int l = li - max_hkl; Reflection r{}; if (!compute_reflection(*kc, h, k, l, r)) return; - int pos = atomicAdd(counter, 1); + // See the rotation kernel: clamping the counter hides an overflow from the host. + const int pos = atomicAdd(counter, 1); if (pos < max_reflections) out[pos] = r; - else atomicSub(counter, 1); } inline KernelConsts BuildKernelConsts(const DiffractionExperiment &experiment, @@ -212,7 +212,10 @@ int BraggPredictionGPU::Calc(const DiffractionExperiment &experiment, cudaStreamSynchronize(stream); int count = *h_count.get(); - if (count > max_reflections) count = max_reflections; + if (count > max_reflections) { + ReportOverflow(count); + count = max_reflections; + } if (count == 0) return {}; diff --git a/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu b/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu index 5d317abe..c30d7b5f 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu +++ b/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu @@ -189,11 +189,11 @@ namespace { int n = compute_reflections_rot(*kc, h, k, l, r); for (int i = 0; i < n; ++i) { - int pos = atomicAdd(counter, 1); + // Do NOT clamp the counter back down on overflow: it then saturates at the capacity and the + // host cannot tell a full buffer from an overflowing one. Let it count the true total. + const int pos = atomicAdd(counter, 1); if (pos < max_reflections) out[pos] = r[i]; - else - atomicSub(counter, 1); } } @@ -289,7 +289,10 @@ int BraggPredictionRotGPU::Calc(const DiffractionExperiment &experiment, cudaStreamSynchronize(stream); int count = *h_count.get(); - if (count > max_reflections) count = max_reflections; + if (count > max_reflections) { + ReportOverflow(count); + count = max_reflections; + } if (count == 0) return 0; cudaMemcpyAsync(reflections.data(), d_out, sizeof(Reflection) * count, cudaMemcpyDeviceToHost, stream);