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

@@ -14,12 +14,22 @@
 #include <string>
 #include <type_traits>
 
+#ifdef AARE_TESTS
+#include "../tests/friend_test.hpp"
+#endif
+
 namespace aare::cuda {
 
 template <typename StorageType = uint16_t,
           typename DeviceStorageType = uint32_t, typename AxisType = float,
           int BLOCK_SIZE = 128>
 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<StorageType, 3> out) const {
         // Make sure the handle points to a well-shaped owning NDArray in
         // Pixel-Major layout

src/hist/PixelHistogramCUDA.test.cpp

@@ -2,7 +2,9 @@
 #include "PixelHistogramCUDA_testimpl.hpp"
 #include "aare/hist/PixelHistogram.hpp"
 #include <catch2/catch_template_test_macros.hpp>
+#include <catch2/generators/catch_generators.hpp>
 #include <cmath>
+#include <limits>
 #include <random>
 #include <vector>
 
@@ -41,23 +43,182 @@ NDArray<float, 3> make_synth(size_t nrows, size_t ncols, size_t nframes) {
     return stack;
 }
 
-} // namespace
+template <typename T>
+void fill_cpu(PixelHistogram<T> &cpu, NDView<const float, 3> 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<float, 2>(NDView<const float, 2>(
+            frames.data() + f * nrows * ncols,
+            {static_cast<ssize_t>(nrows), static_cast<ssize_t>(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<TestType> cpu(nrows, ncols, nbins, xmin, xmax);
-    for (size_t f = 0; f < nframes; ++f)
-        cpu.fill_async(NDArray<float, 2>(NDView<float, 2>(
-            frames.data() + f * nrows * ncols,
-            {static_cast<ssize_t>(nrows), static_cast<ssize_t>(ncols)})));
+    PixelHistogram<TestType> cpu(nrows, ncols, nbins, xmin, xmax,
+                                 /*n_threads=*/8);
+    fill_cpu(cpu, frames.view());
 
     auto gpu = test_cuda::run_gpu<TestType>(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<TestType>::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<float, 3> frames({static_cast<ssize_t>(nframes),
+                              static_cast<ssize_t>(nrows),
+                              static_cast<ssize_t>(ncols)},
+                             value);
+
+    PixelHistogram<TestType> cpu(nrows, ncols, nbins, xmin, xmax);
+    fill_cpu(cpu, frames.view());
+
+    auto gpu = test_cuda::run_gpu<TestType>(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<size_t>(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<float>::infinity();
+    // One pixel per case, the same value in every frame.
+    const std::vector<float> 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<float>::quiet_NaN(),
+        inf,
+        -inf,
+    };
+    const size_t ncols = cases.size(), nframes = 10, nbins = 120;
+
+    NDArray<float, 3> frames(
+        {static_cast<ssize_t>(nframes), 1, static_cast<ssize_t>(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<uint32_t, 3> expected(
+        {1, static_cast<ssize_t>(ncols), static_cast<ssize_t>(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<uint32_t>(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<uint32_t> cpu(nrows, ncols, nbins, xmin, xmax);
+    fill_cpu(cpu, frames.view());
+
+    auto gpu = test_cuda::run_gpu<uint32_t>(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<size_t> 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<uint32_t> 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<ssize_t>(begin),
+                                             static_cast<ssize_t>(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<std::pair<size_t, size_t>> 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<uint8_t, 1> slice({static_cast<ssize_t>(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());
 }

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 <algorithm>
+#include <limits>
+#include <numeric>
+#include <stdexcept>
+#include <string>
 #include <type_traits>
+#include <vector>
+
+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<StorageType> hist(nrows, ncols, nbins, 0.0f, 1.0f);
+        const size_t npix = nrows * ncols;
+
+        std::vector<uint32_t> 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<StorageType>)
+                    want = std::min<uint64_t>(
+                        want, std::numeric_limits<StorageType>::max());
+                const auto got = out(p / ncols, p % ncols, b);
+                if (got != static_cast<StorageType>(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<double>(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 <typename H>
 DeviceResult<H> run_gpu(NDView<float, 3> frames, size_t nbins, float xmin,
                         float xmax, size_t batch_size) {
@@ -29,4 +91,60 @@ INSTANTIATE(float)
 INSTANTIATE(double)
 #undef INSTANTIATE
 
+std::vector<NDArray<uint32_t, 3>>
+run_gpu_in_parts(NDView<float, 3> frames, size_t nbins, float xmin, float xmax,
+                 size_t batch_size, const std::vector<size_t> &ends) {
+    cuda::PixelHistogramCUDA<uint32_t> hist(
+        static_cast<size_t>(frames.shape(1)),
+        static_cast<size_t>(frames.shape(2)), nbins, xmin, xmax, batch_size);
+    hist.register_input_buffer(frames.data(), frames.total_bytes());
+
+    std::vector<NDArray<uint32_t, 3>> snapshots;
+    size_t begin = 0;
+    for (size_t end : ends) {
+        hist.fill_histogram(frames.sub_view(static_cast<ssize_t>(begin),
+                                            static_cast<ssize_t>(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 <typename Hist>
+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<float, 3> frames, size_t nbins, float xmin,
+                     float xmax, size_t batch_size,
+                     const std::vector<std::pair<size_t, size_t>> &pixels) {
+    const auto nrows = static_cast<size_t>(frames.shape(1));
+    const auto ncols = static_cast<size_t>(frames.shape(2));
+    cuda::PixelHistogramCUDA<uint8_t> 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

src/hist/PixelHistogramCUDA_testimpl.hpp

@@ -4,6 +4,8 @@
 #include "aare/NDView.hpp"
 #include <cstddef>
 #include <cstdint>
+#include <utility>
+#include <vector>
 
 namespace aare::test_cuda {
 
@@ -18,4 +20,29 @@ template <typename H>
 DeviceResult<H> run_gpu(NDView<float, 3> 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<NDArray<uint32_t, 3>>
+run_gpu_in_parts(NDView<float, 3> frames, size_t nbins, float xmin, float xmax,
+                 size_t batch_size, const std::vector<size_t> &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<uint8_t, 3> values;
+    std::vector<NDArray<uint8_t, 1>> 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<float, 3> frames, size_t nbins, float xmin,
+                     float xmax, size_t batch_size,
+                     const std::vector<std::pair<size_t, size_t>> &pixels);
+
 } // namespace aare::test_cuda