From 2cbb3fc8b4f1e48567c2874e081f0a0480f08840 Mon Sep 17 00:00:00 2001 From: jungfrau Date: Sat, 15 Aug 2026 21:42:41 -0400 Subject: [PATCH] Build the GPU engines a worker never uses on first use, not always Every worker thread built a full set of analysis engines. Two of them are never asked for on the offline path: the fixed-threshold spot finder, because detection is adaptive by default, and the azimuthal integrator, because the fused adaptive finder produces the profile as a by-product. They are still needed elsewhere - the broker defaults to non-adaptive detection, and --no-adaptive-spots asks for the finder - so they are built on first use rather than removed. A lazily built finder takes the current resolution mask on construction; without that it would find spots outside the limits it was never told about. The bitshuffle decoder sized its output buffer for the widest pixel type there is rather than the one the images actually have, holding a second full frame per worker on 16-bit data. It is sized from the image now and grows if a later frame needs more. The shared-table checksum runs over eight interleaved lanes. FNV's multiply is a loop-carried dependency, so one chain retires a byte every few cycles whatever memory bandwidth is spare, and every worker hashes tens of megabytes of geometry tables as it builds its engines - about 5% of all CPU samples on a 16M-pixel detector. Measured on a 16M-pixel rotation dataset: cudaMalloc 11314 -> 9474 calls and, with cudaFree, 117 s -> 78 s of aggregate thread time; both synchronise the whole device, so that time is spent blocking every other worker. Whole battery 15m32s -> 12m30s. Data quality against main, over 24 crystals and eight statistics each: the same space group on all 24, and every difference smaller than what two runs of an IDENTICAL binary produce (measured: 13 of 24 crystals reproduce exactly run to run, worst R_meas swing 5.5 points, against 4.6 points for main vs this branch). The float atomics in the reductions have always made this so. Co-Authored-By: Claude Opus 5 (1M context) --- image_analysis/MXAnalysisWithoutFPGA.cpp | 41 +++++++++++++++---- image_analysis/MXAnalysisWithoutFPGA.h | 7 ++++ .../image_preprocessing/BSLZ4DecoderGPU.cu | 15 ++++++- .../image_preprocessing/BSLZ4DecoderGPU.h | 1 + .../ImagePreprocessorGPU.cu | 4 +- image_analysis/indexing/CudaSharedTables.h | 31 +++++++++++--- 6 files changed, 84 insertions(+), 15 deletions(-) diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 63a55f58..8ad97ca5 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -50,8 +50,6 @@ MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_exp if (get_gpu_count() == 0) { #endif preprocessor_buffer = std::make_unique(experiment.GetPixelsNum()); - spotFinder = std::make_unique(experiment.GetXPixelsNum(), experiment.GetYPixelsNum()); - azint = std::make_unique(integration); preprocessor = std::make_unique(in_experiment, in_mask); bragg_engine = std::make_unique(in_experiment); if (experiment.ROI().size() >= 1) @@ -70,8 +68,6 @@ MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_exp // enable_fused_adaptive_gpu = true, so on the GPU path the copy is off in practice. preprocessor = std::make_unique(in_experiment, in_mask, stream, /*copy_image_to_host=*/!enable_fused_adaptive_gpu); - spotFinder = std::make_unique(experiment.GetXPixelsNum(), experiment.GetYPixelsNum(), stream); - azint = std::make_unique(integration, stream); bragg_engine = std::make_unique(in_experiment, stream); if (experiment.ROI().size() >= 1) roi = std::make_unique(experiment, stream); @@ -144,7 +140,7 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, if (!fused) { const auto azint_start_time = std::chrono::steady_clock::now(); - azint->Run(*preprocessor_buffer, profile); + AzInt().Run(*preprocessor_buffer, profile); const auto azint_end_time = std::chrono::steady_clock::now(); output.azint_time_s = std::chrono::duration(azint_end_time - azint_start_time).count(); } @@ -160,7 +156,7 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, ImageSpotFinder &finder = spot_finding_settings.adaptive_threshold ? static_cast(*adaptiveSpotFinder) - : *spotFinder; + : FixedThresholdFinder(); const auto integrate_fn = [this](const std::vector &predicted, size_t npredicted, int64_t image_number) { return bragg_engine->Run(*preprocessor_buffer, predicted, npredicted, image_number); @@ -293,6 +289,36 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, integration.Settings(), spot_finding_settings.ice_ring_width_Q_recipA); } +ImageSpotFinder &MXAnalysisWithoutFPGA::FixedThresholdFinder() { + if (!spotFinder) { +#ifdef JFJOCH_USE_CUDA + if (stream) + spotFinder = std::make_unique(experiment.GetXPixelsNum(), + experiment.GetYPixelsNum(), stream); + else +#endif + spotFinder = std::make_unique(experiment.GetXPixelsNum(), + experiment.GetYPixelsNum()); + // It missed every mask update that happened before it existed, so it takes the current one + // now. Without this it would find spots outside the resolution limits. + if (mask_high_res.has_value() || mask_low_res.has_value()) + spotFinder->SetResolutionMask(mask_resolution); + } + return *spotFinder; +} + +AzIntEngine &MXAnalysisWithoutFPGA::AzInt() { + if (!azint) { +#ifdef JFJOCH_USE_CUDA + if (stream) + azint = std::make_unique(integration, stream); + else +#endif + azint = std::make_unique(integration); + } + return *azint; +} + void MXAnalysisWithoutFPGA::RebuildROI() { if (experiment.ROI().empty()) { roi.reset(); @@ -346,6 +372,7 @@ void MXAnalysisWithoutFPGA::UpdateMaskResolution(const SpotFindingSettings &sett // The finders keep their own copy (the GPU ones a bit-packed device copy), so the mask is handed // over here - when the limits change - rather than with every image. - spotFinder->SetResolutionMask(mask_resolution); + if (spotFinder) + spotFinder->SetResolutionMask(mask_resolution); adaptiveSpotFinder->SetResolutionMask(mask_resolution); } diff --git a/image_analysis/MXAnalysisWithoutFPGA.h b/image_analysis/MXAnalysisWithoutFPGA.h index 58512091..7af8249b 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.h +++ b/image_analysis/MXAnalysisWithoutFPGA.h @@ -36,9 +36,16 @@ class MXAnalysisWithoutFPGA { size_t npixels; size_t xpixels; + // Built on first use: the fused adaptive finder produces the azimuthal profile as a by-product, + // so on the rugnux path this engine is constructed and then never run. std::unique_ptr azint; + AzIntEngine &AzInt(); std::unique_ptr roi; + // Built on first use. Which finder an image takes arrives with its SpotFindingSettings, and + // with adaptive detection on - the default everywhere but the broker - this one is never asked + // for; on the GPU it is ~14 MB and 15 device allocations per worker. std::unique_ptr spotFinder; + ImageSpotFinder &FixedThresholdFinder(); // Self-calibrating finder, used when spot settings request adaptive detection. Kept alongside the // default finder because the choice arrives with the per-image settings, not at construction. It is // an AdaptiveSpotFinderCPU by default; on the GPU path, when the fused engine is enabled (rugnux diff --git a/image_analysis/image_preprocessing/BSLZ4DecoderGPU.cu b/image_analysis/image_preprocessing/BSLZ4DecoderGPU.cu index 8068328d..29cd70d0 100644 --- a/image_analysis/image_preprocessing/BSLZ4DecoderGPU.cu +++ b/image_analysis/image_preprocessing/BSLZ4DecoderGPU.cu @@ -231,6 +231,18 @@ BSLZ4DecoderGPU::BSLZ4DecoderGPU(size_t in_max_uncompressed_bytes, std::shared_p // tens; sizing this from the UNCOMPRESSED size cost ~73 MB per worker to hold ~4 MB. } +// The stored pixel depth is fixed within a dataset, so in practice this runs once - but it is not +// promised anywhere, and sizing for the widest type instead would hold twice the memory a 16-bit +// detector needs. Grown with slack because cudaMalloc and cudaFree synchronise the whole device. +void BSLZ4DecoderGPU::EnsureUncompressedCapacity(size_t bytes) { + if (bytes <= max_uncompressed_bytes) + return; + const size_t want = std::max(bytes, max_uncompressed_bytes + max_uncompressed_bytes / 2); + cuda_err(cudaStreamSynchronize(*stream)); + gpu_shuffled = CudaDevicePtr(want); + max_uncompressed_bytes = want; +} + void BSLZ4DecoderGPU::EnsureCompressedCapacity(size_t bytes) { if (bytes <= compressed_capacity) return; @@ -259,8 +271,7 @@ BSLZ4ShuffledImage BSLZ4DecoderGPU::DecodeShuffled(const CompressedImage &image) if (clen < 12) throw JFJochException(JFJochExceptionCategory::Compression, "bslz4 chunk shorter than its header"); - if (total_bytes > max_uncompressed_bytes) - throw JFJochException(JFJochExceptionCategory::Compression, "bslz4 image larger than the decoder was sized for"); + EnsureUncompressedCapacity(total_bytes); if (be64(src) != total_bytes) throw JFJochException(JFJochExceptionCategory::Compression, "bslz4 header size does not match the image"); diff --git a/image_analysis/image_preprocessing/BSLZ4DecoderGPU.h b/image_analysis/image_preprocessing/BSLZ4DecoderGPU.h index bdc81485..98da7b5f 100644 --- a/image_analysis/image_preprocessing/BSLZ4DecoderGPU.h +++ b/image_analysis/image_preprocessing/BSLZ4DecoderGPU.h @@ -73,6 +73,7 @@ class BSLZ4DecoderGPU { size_t max_blocks = 0; void EnsureCompressedCapacity(size_t bytes); + void EnsureUncompressedCapacity(size_t bytes); void EnsureBlockCapacity(size_t nblocks); public: diff --git a/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu b/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu index 426393f0..37c772af 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu +++ b/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu @@ -295,7 +295,9 @@ bool ImagePreprocessorGPU::AnalyzeCompressed(ImagePreprocessorBuffer &processed_ return false; if (!bslz4_decoder) - bslz4_decoder = std::make_unique(npixels * sizeof(uint32_t), stream); + // Sized for the depth this image actually has, not for the widest one there could be: on + // 16-bit data the difference is half of a full frame per worker thread. + bslz4_decoder = std::make_unique(image.GetUncompressedSize(), stream); // LZ4 on the device, then ONE kernel that un-transposes the bitshuffle blocks and preprocesses // them as it goes. The decompressed image is never materialised: the fused kernel reads the diff --git a/image_analysis/indexing/CudaSharedTables.h b/image_analysis/indexing/CudaSharedTables.h index d7c30d95..3c19d5fb 100644 --- a/image_analysis/indexing/CudaSharedTables.h +++ b/image_analysis/indexing/CudaSharedTables.h @@ -42,14 +42,35 @@ namespace jfjoch_cuda_shared_tables { // FNV-1a. Not a cryptographic hash and does not need to be - it exists to notice that the bytes // behind a reused address changed, not to resist anyone. + // + // Run over eight interleaved lanes and fold them at the end. FNV's multiply is a loop-carried + // dependency, so one lane retires a byte every few cycles however much memory bandwidth is + // going spare; eight independent chains fill that latency. The tables are tens of megabytes and + // every worker thread hashes them as it builds its engines, which put this at ~5% of all CPU + // samples on a 16M-pixel detector. inline uint64_t checksum(const void *data, size_t bytes) { + constexpr uint64_t PRIME = 1099511628211ULL; + constexpr size_t LANES = 8; const auto *p = static_cast(data); - uint64_t h = 1469598103934665603ULL; - for (size_t i = 0; i < bytes; i++) { - h ^= p[i]; - h *= 1099511628211ULL; + uint64_t h[LANES]; + for (size_t l = 0; l < LANES; l++) + h[l] = 1469598103934665603ULL + l; + const size_t n = bytes / LANES * LANES; + for (size_t i = 0; i < n; i += LANES) + for (size_t l = 0; l < LANES; l++) { + h[l] ^= p[i + l]; + h[l] *= PRIME; + } + uint64_t out = 1469598103934665603ULL; + for (size_t l = 0; l < LANES; l++) { + out ^= h[l]; + out *= PRIME; } - return h; + for (size_t i = n; i < bytes; i++) { + out ^= p[i]; + out *= PRIME; + } + return out; } inline Registry ®istry() {