From f923c71bf5f89dd78234b804b949fe11e2e3e4ee Mon Sep 17 00:00:00 2001 From: kferjaoui Date: Wed, 21 Oct 2026 12:21:31 +0200 Subject: [PATCH] Add PixelHistogramCUDA tests for edge cases and the readback path - Saturation: uint8/uint16 storage capped at the type's max, as on CPU. - Range: xmin and the last float below xmax are kept; xmax, below xmin, NaN and +-inf are dropped (checked against expected counts, since PixelHistogram still counts NaN in bin 0). - Batch size 1, 7, 25, 100 and 150 give the CPU result. - Several fill_histogram() calls with values() snapshots in between. - values_pixel() equals the pixel's slice of values(), including capped bins, and throws std::out_of_range for bad coordinates. - Known pattern written into the accumulator through a FRIEND_TEST hook (AARE_TESTS only) checks the transpose and cap without the fill kernel. - Comparison tests use non-square 7x13 frames. --- include/aare/hist/PixelHistogramCUDA.hpp | 11 +- src/hist/PixelHistogramCUDA.test.cpp | 179 +++++++++++++++++++++-- src/hist/PixelHistogramCUDA_testimpl.cu | 118 +++++++++++++++ src/hist/PixelHistogramCUDA_testimpl.hpp | 27 ++++ 4 files changed, 325 insertions(+), 10 deletions(-) diff --git a/include/aare/hist/PixelHistogramCUDA.hpp b/include/aare/hist/PixelHistogramCUDA.hpp index 6a9837ca..b332196b 100644 --- a/include/aare/hist/PixelHistogramCUDA.hpp +++ b/include/aare/hist/PixelHistogramCUDA.hpp @@ -14,12 +14,22 @@ #include #include +#ifdef AARE_TESTS +#include "../tests/friend_test.hpp" +#endif + namespace aare::cuda { template class PixelHistogramCUDA { +#ifdef AARE_TESTS + // Writes a known pattern into d_histogram to test values() without the + // fill kernel (PixelHistogramCUDA_testimpl.cu). + FRIEND_TEST(known_pattern_transpose_and_cap) +#endif + // nrows is declared first: its initializer validates all constructor // arguments, so it must run before anything below allocates. size_t nrows; @@ -199,7 +209,6 @@ class PixelHistogramCUDA { // In theory, this fct captures the histogram built from batches processed // so far i.e. possibly runs between 2 consecutive batches and before all // frames where processed... - // TODO: Needs transpose form Bin-Major to Pixel-Major void values(NDView out) const { // Make sure the handle points to a well-shaped owning NDArray in // Pixel-Major layout diff --git a/src/hist/PixelHistogramCUDA.test.cpp b/src/hist/PixelHistogramCUDA.test.cpp index c166ee8c..1e53ad78 100644 --- a/src/hist/PixelHistogramCUDA.test.cpp +++ b/src/hist/PixelHistogramCUDA.test.cpp @@ -2,7 +2,9 @@ #include "PixelHistogramCUDA_testimpl.hpp" #include "aare/hist/PixelHistogram.hpp" #include +#include #include +#include #include #include @@ -41,23 +43,182 @@ NDArray make_synth(size_t nrows, size_t ncols, size_t nframes) { return stack; } -} // namespace +template +void fill_cpu(PixelHistogram &cpu, NDView frames) { + const size_t nframes = frames.shape(0); + const size_t nrows = frames.shape(1); + const size_t ncols = frames.shape(2); + for (size_t f = 0; f < nframes; ++f) + cpu.fill_async(NDArray(NDView( + frames.data() + f * nrows * ncols, + {static_cast(nrows), static_cast(ncols)}))); +} -TEMPLATE_TEST_CASE("PixelHistogramCUDA matches PixelHistogram", "[cuda]", - uint8_t, uint16_t, uint32_t, uint64_t, float, double) { - const size_t nrows = 5, ncols = 5, nframes = 100, nbins = 120; +} // anonymous namespace + +TEMPLATE_TEST_CASE("PixelHistogramCUDA matches PixelHistogram on noisy frames", + "[cuda]", uint8_t, uint16_t, uint32_t, uint64_t, float, + double) { + // Non-square, so swapped rows and columns or a wrong stride show up. + const size_t nrows = 7, ncols = 13, nframes = 100, nbins = 120; const float xmin = 3900.0f, xmax = 5100.0f; auto frames = make_synth(nrows, ncols, nframes); - PixelHistogram cpu(nrows, ncols, nbins, xmin, xmax); - for (size_t f = 0; f < nframes; ++f) - cpu.fill_async(NDArray(NDView( - frames.data() + f * nrows * ncols, - {static_cast(nrows), static_cast(ncols)}))); + PixelHistogram cpu(nrows, ncols, nbins, xmin, xmax, + /*n_threads=*/8); + fill_cpu(cpu, frames.view()); auto gpu = test_cuda::run_gpu(frames.view(), nbins, xmin, xmax, /*batch_size=*/8); REQUIRE((cpu.bin_edges() == gpu.edges)); REQUIRE((cpu.bin_centers() == gpu.centers)); REQUIRE((cpu.values() == gpu.values)); +} + +TEMPLATE_TEST_CASE( + "PixelHistogramCUDA saturates at the storage maximum like PixelHistogram", + "[cuda]", uint8_t, uint16_t) { + const size_t cap = std::numeric_limits::max(); + const size_t nrows = 1, ncols = 1, nframes = 2 * cap, nbins = 120; + const float xmin = 3900.0f, xmax = 5100.0f; + + // Centre of bin 60: 4505.0f + const size_t bin_idx = 60; + const float value = xmin + (bin_idx + 0.5f) * (xmax - xmin) / nbins; + + NDArray frames({static_cast(nframes), + static_cast(nrows), + static_cast(ncols)}, + value); + + PixelHistogram cpu(nrows, ncols, nbins, xmin, xmax); + fill_cpu(cpu, frames.view()); + + auto gpu = test_cuda::run_gpu(frames.view(), nbins, xmin, xmax, + /*batch_size=*/1024); + + REQUIRE((cpu.values() == gpu.values)); + + auto gpu_count = gpu.values(0, 0, bin_idx); + + REQUIRE(static_cast(gpu_count) == cap); +} + +TEST_CASE("PixelHistogramCUDA keeps [xmin, xmax) and drops NaN and infinities", + "[cuda]") { + const float xmin = 3900.0f, xmax = 5100.0f; + const float inf = std::numeric_limits::infinity(); + // One pixel per case, the same value in every frame. + const std::vector cases = { + xmin, // first bin + std::nextafter(xmax, xmin), // last bin + xmax, // dropped: the range is half-open + std::nextafter(xmin, -inf), // dropped + std::numeric_limits::quiet_NaN(), + inf, + -inf, + }; + const size_t ncols = cases.size(), nframes = 10, nbins = 120; + + NDArray frames( + {static_cast(nframes), 1, static_cast(ncols)}); + for (size_t f = 0; f < nframes; ++f) + for (size_t c = 0; c < ncols; ++c) + frames(f, 0, c) = cases[c]; + + // Pixel-major. Pixel (0,0) counts in the first bin, pixel (0,1) in the + // last; all other pixels stay empty. + NDArray expected( + {1, static_cast(ncols), static_cast(nbins)}, 0); + expected(0, 0, 0) = nframes; + expected(0, 1, nbins - 1) = nframes; + + // Compared with the expected counts rather than with PixelHistogram: + // the CPU class still counts NaN in bin 0. + auto gpu = test_cuda::run_gpu(frames.view(), nbins, xmin, xmax, + /*batch_size=*/4); + REQUIRE((gpu.values == expected)); +} + +TEST_CASE("PixelHistogramCUDA result does not depend on the batch size", + "[cuda]") { + // 1: one frame per batch, both slots reused many times. 7: partial last + // batch. 25: divides nframes. 100: exactly one batch. 150: one batch + // smaller than its slot. + const size_t batch_size = GENERATE(1, 7, 25, 100, 150); + CAPTURE(batch_size); + + const size_t nrows = 7, ncols = 13, nframes = 100, nbins = 120; + const float xmin = 3900.0f, xmax = 5100.0f; + auto frames = make_synth(nrows, ncols, nframes); + + PixelHistogram cpu(nrows, ncols, nbins, xmin, xmax); + fill_cpu(cpu, frames.view()); + + auto gpu = test_cuda::run_gpu(frames.view(), nbins, xmin, xmax, + batch_size); + REQUIRE((cpu.values() == gpu.values)); +} + +TEST_CASE("PixelHistogramCUDA accumulates over several fills with values() " + "in between", + "[cuda]") { + const size_t nrows = 7, ncols = 13, nframes = 100, nbins = 120; + const float xmin = 3900.0f, xmax = 5100.0f; + auto frames = make_synth(nrows, ncols, nframes); + + // With batch_size 8: part 1 ends in a partial batch, part 2 is a single + // frame, part 3 fills the rest after two values() calls. + const std::vector ends = {30, 31, nframes}; + auto snapshots = test_cuda::run_gpu_in_parts(frames.view(), nbins, xmin, + xmax, /*batch_size=*/8, ends); + REQUIRE(snapshots.size() == ends.size()); + + PixelHistogram cpu(nrows, ncols, nbins, xmin, xmax); + size_t begin = 0; + for (size_t i = 0; i < ends.size(); ++i) { + CAPTURE(ends[i]); + fill_cpu(cpu, frames.view().sub_view(static_cast(begin), + static_cast(ends[i]))); + REQUIRE((cpu.values() == snapshots[i])); + begin = ends[i]; + } +} + +TEST_CASE("PixelHistogramCUDA values_pixel matches the pixel's slice of " + "values()", + "[cuda]") { + // 2000 frames push peak bins past 255, so the uint8_t cap in values_pixel + // is exercised next to uncapped tail bins (saw_cap checks it happens). + const size_t nrows = 7, ncols = 13, nframes = 2000, nbins = 120; + const float xmin = 3900.0f, xmax = 5100.0f; + auto frames = make_synth(nrows, ncols, nframes); + + const std::vector> pixels = { + {0, 0}, {0, ncols - 1}, {nrows - 1, 0}, {nrows - 1, ncols - 1}, {3, 6}}; + auto gpu = test_cuda::run_gpu_values_pixel(frames.view(), nbins, xmin, xmax, + /*batch_size=*/64, pixels); + REQUIRE(gpu.pixels.size() == pixels.size()); + + bool saw_cap = false; + for (size_t i = 0; i < pixels.size(); ++i) { + const auto [row, col] = pixels[i]; + CAPTURE(row, col); + NDArray slice({static_cast(nbins)}); + for (size_t b = 0; b < nbins; ++b) { + slice(b) = gpu.values(row, col, b); + saw_cap = saw_cap || slice(b) == 255; + } + REQUIRE((gpu.pixels[i] == slice)); + } + REQUIRE(saw_cap); + + REQUIRE(gpu.row_out_of_range_throws); + REQUIRE(gpu.col_out_of_range_throws); +} + +TEST_CASE("PixelHistogramCUDA values() transposes and caps a known " + "accumulator pattern", + "[cuda]") { + REQUIRE_NOTHROW(test_cuda::check_known_pattern()); } \ No newline at end of file diff --git a/src/hist/PixelHistogramCUDA_testimpl.cu b/src/hist/PixelHistogramCUDA_testimpl.cu index fb540799..0194cb14 100644 --- a/src/hist/PixelHistogramCUDA_testimpl.cu +++ b/src/hist/PixelHistogramCUDA_testimpl.cu @@ -1,10 +1,72 @@ // SPDX-License-Identifier: MPL-2.0 #include "PixelHistogramCUDA_testimpl.hpp" #include "aare/hist/PixelHistogramCUDA.hpp" +#include +#include +#include +#include +#include #include +#include + +namespace aare::cuda { + +// Writes count = b * npix + p, unique per cell, into the bin-major +// accumulator, then checks that values() returns it pixel-major and capped +// at StorageType's max. Throws on the first mismatch. The check is a lambda +// because only this friend function may touch the private members. +void known_pattern_transpose_and_cap_impl() { + auto check = [](auto storage_tag, size_t nrows, size_t ncols, size_t nbins, + const char *type_name) { + using StorageType = decltype(storage_tag); + PixelHistogramCUDA hist(nrows, ncols, nbins, 0.0f, 1.0f); + const size_t npix = nrows * ncols; + + std::vector pattern(npix * nbins); + // Bin-major counts 0 .. npix * nbins - 1: within a bin, successive + // pixels hold successive values; one pixel's bins are npix apart. + std::iota(pattern.begin(), pattern.end(), 0u); + // Copied on the compute stream: after the constructor's memset, + // before values()'s transpose. + CUDA_CHECK(cudaMemcpyAsync( + hist.d_histogram.get(), pattern.data(), hist.d_histogram.bytes(), + cudaMemcpyHostToDevice, hist.compute_stream.get())); + + const auto out = hist.values(); + for (size_t p = 0; p < npix; ++p) { + for (size_t b = 0; b < nbins; ++b) { + uint64_t want = b * npix + p; + if constexpr (std::is_integral_v) + want = std::min( + want, std::numeric_limits::max()); + const auto got = out(p / ncols, p % ncols, b); + if (got != static_cast(want)) + throw std::runtime_error( + std::string(type_name) + ", " + std::to_string(nrows) + + "x" + std::to_string(ncols) + " px, " + + std::to_string(nbins) + " bins: pixel " + + std::to_string(p) + ", bin " + std::to_string(b) + + ": expected " + std::to_string(want) + ", got " + + std::to_string(static_cast(got))); + } + } + }; + + // Odd sizes leave partial blocks in the transpose. Counts reach ~88k: + // capped for uint8_t and uint16_t, exact for uint32_t and float. + check(uint8_t{}, 37, 53, 45, "uint8_t"); + check(uint16_t{}, 37, 53, 45, "uint16_t"); + check(uint32_t{}, 37, 53, 45, "uint32_t"); + check(float{}, 37, 53, 45, "float"); + check(uint32_t{}, 1, 1, 1, "uint32_t"); +} + +} // namespace aare::cuda namespace aare::test_cuda { +void check_known_pattern() { cuda::known_pattern_transpose_and_cap_impl(); } + template DeviceResult run_gpu(NDView frames, size_t nbins, float xmin, float xmax, size_t batch_size) { @@ -29,4 +91,60 @@ INSTANTIATE(float) INSTANTIATE(double) #undef INSTANTIATE +std::vector> +run_gpu_in_parts(NDView frames, size_t nbins, float xmin, float xmax, + size_t batch_size, const std::vector &ends) { + cuda::PixelHistogramCUDA hist( + static_cast(frames.shape(1)), + static_cast(frames.shape(2)), nbins, xmin, xmax, batch_size); + hist.register_input_buffer(frames.data(), frames.total_bytes()); + + std::vector> snapshots; + size_t begin = 0; + for (size_t end : ends) { + hist.fill_histogram(frames.sub_view(static_cast(begin), + static_cast(end))); + snapshots.push_back(hist.values()); + begin = end; + } + hist.unregister_input_buffer(); + return snapshots; +} + +namespace { + +// The exception can't cross into the Catch2 file as a REQUIRE_THROWS, so +// report whether values_pixel(row, col) threw std::out_of_range. +template +bool throws_out_of_range(const Hist &hist, size_t row, size_t col) { + try { + hist.values_pixel(row, col); + } catch (const std::out_of_range &) { + return true; + } + return false; +} + +} // namespace + +PixelQueries +run_gpu_values_pixel(NDView frames, size_t nbins, float xmin, + float xmax, size_t batch_size, + const std::vector> &pixels) { + const auto nrows = static_cast(frames.shape(1)); + const auto ncols = static_cast(frames.shape(2)); + cuda::PixelHistogramCUDA hist(nrows, ncols, nbins, xmin, xmax, + batch_size); + hist.register_input_buffer(frames.data(), frames.total_bytes()); + hist.fill_histogram(frames); + hist.unregister_input_buffer(); + + PixelQueries out{hist.values(), {}, false, false}; + for (const auto &[row, col] : pixels) + out.pixels.push_back(hist.values_pixel(row, col)); + out.row_out_of_range_throws = throws_out_of_range(hist, nrows, 0); + out.col_out_of_range_throws = throws_out_of_range(hist, 0, ncols); + return out; +} + } // namespace aare::test_cuda diff --git a/src/hist/PixelHistogramCUDA_testimpl.hpp b/src/hist/PixelHistogramCUDA_testimpl.hpp index 18e8a0c3..6f29fb6c 100644 --- a/src/hist/PixelHistogramCUDA_testimpl.hpp +++ b/src/hist/PixelHistogramCUDA_testimpl.hpp @@ -4,6 +4,8 @@ #include "aare/NDView.hpp" #include #include +#include +#include namespace aare::test_cuda { @@ -18,4 +20,29 @@ template DeviceResult run_gpu(NDView frames, size_t nbins, float xmin, float xmax, size_t batch_size); +// Fill one histogram from consecutive parts of `frames`, part i ending at +// frame ends[i], and take a values() snapshot after each part. +std::vector> +run_gpu_in_parts(NDView frames, size_t nbins, float xmin, float xmax, + size_t batch_size, const std::vector &ends); + +// Known pattern written straight into the device accumulator, read back +// through values(); tests the transpose and the cap without the fill kernel. +// Throws std::runtime_error describing the first mismatch. +void check_known_pattern(); + +struct PixelQueries { + NDArray values; + std::vector> pixels; // values_pixel(), one per query + bool row_out_of_range_throws; // values_pixel(nrows, 0) + bool col_out_of_range_throws; // values_pixel(0, ncols) +}; + +// Fill a uint8_t histogram from `frames`, then read values() and +// values_pixel() for each (row, col) in `pixels`. +PixelQueries +run_gpu_values_pixel(NDView frames, size_t nbins, float xmin, + float xmax, size_t batch_size, + const std::vector> &pixels); + } // namespace aare::test_cuda