Test PixelHistogramCUDA against PixelHistogram

Wire PixelHistogramCUDA into CMake and add a Catch2 test: 100 synthetic
5x5 frames histogrammed on CPU and GPU (batch_size 8, partial last
batch) must give identical edges, centers and values for uint8/16/32/64,
float and double storage. 8-byte storage uses uint64_t device counters.

PixelHistogramCUDA: values_pixel() throws std::out_of_range; add rows(),
cols(), n_bins() and batch_size().

CMakeLists.txt

@@ -423,6 +423,8 @@ if(AARE_CUDA)
     include/aare/cf/ClusterFinderCUDA.hpp
     include/aare/cf/clusterfinder_algo.cuh
     include/aare/cf/clusterfinder_kernel.cuh
+    include/aare/hist/PixelHistogramCUDA.hpp
+    include/aare/hist/pixelhistogram_kernel.cuh
     include/aare/RAII_device.hpp
     include/aare/utils/cuda_check.cuh)
 endif()
@@ -549,7 +551,9 @@ if(AARE_TESTS AND AARE_CUDA)
   target_sources(
     tests
     PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/src/cf/ClusterFinderCUDA.test.cpp
-            ${CMAKE_CURRENT_SOURCE_DIR}/src/cf/ClusterFinderCUDA_testimpl.cu)
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/cf/ClusterFinderCUDA_testimpl.cu
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/hist/PixelHistogramCUDA.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/hist/PixelHistogramCUDA_testimpl.cu)
   target_link_libraries(tests PRIVATE aare_cuda)
 endif()
 

include/aare/hist/PixelHistogramCUDA.hpp

@@ -261,7 +261,7 @@ class PixelHistogramCUDA {
     // Needs transpose form Bin-Major to Pixel-Major
     NDArray<StorageType, 1> values_pixel(size_t row, size_t col) const {
         if (row >= nrows || col >= ncols)
-            throw std::invalid_argument(
+            throw std::out_of_range(
                 "Incorrect pixel coordinate; out-of bound access");
         NDArray<DeviceStorageType, 1> pixel_hist(
             {static_cast<ssize_t>(m_nbins)});
@@ -310,6 +310,11 @@ class PixelHistogramCUDA {
                        });
         return centers;
     }
+
+    size_t rows() const noexcept { return nrows; }
+    size_t cols() const noexcept { return ncols; }
+    size_t n_bins() const noexcept { return m_nbins; }
+    size_t batch_size() const noexcept { return m_batch_size; }
 };
 
-} // namespace aare::cuda
+} // namespace aare::cuda

src/hist/PixelHistogramCUDA.test.cpp

@@ -0,0 +1,63 @@
+// SPDX-License-Identifier: MPL-2.0
+#include "PixelHistogramCUDA_testimpl.hpp"
+#include "aare/hist/PixelHistogram.hpp"
+#include <catch2/catch_template_test_macros.hpp>
+#include <cmath>
+#include <random>
+#include <vector>
+
+using namespace aare;
+
+namespace {
+
+NDArray<float, 3> make_synth(size_t nrows, size_t ncols, size_t nframes) {
+    std::mt19937 rng(12345);
+    std::uniform_real_distribution<float> base_d(4000.0f, 5000.0f);
+    std::uniform_real_distribution<float> sigma_d(5.0f, 15.0f);
+
+    size_t npix = nrows * ncols;
+    std::vector<float> base(npix);
+    std::vector<std::normal_distribution<float>> noise;
+    noise.reserve(npix);
+
+    for (size_t i = 0; i < npix; ++i) {
+        base[i] = base_d(rng);
+        noise.emplace_back(0.0f, sigma_d(rng));
+    }
+
+    NDArray<float, 3> stack({static_cast<ssize_t>(nframes),
+                             static_cast<ssize_t>(nrows),
+                             static_cast<ssize_t>(ncols)});
+
+    float *data_ = stack.data();
+
+    for (size_t f = 0; f < nframes; ++f) {
+        for (size_t i = 0; i < npix; ++i) {
+            data_[f * npix + i] =
+                std::round(base[i] + noise[i](rng)); // rounding is optional
+        }
+    }
+
+    return stack;
+}
+
+} // namespace
+
+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;
+    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)})));
+
+    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));
+}

src/hist/PixelHistogramCUDA_testimpl.cu

@@ -0,0 +1,32 @@
+// SPDX-License-Identifier: MPL-2.0
+#include "PixelHistogramCUDA_testimpl.hpp"
+#include "aare/hist/PixelHistogramCUDA.hpp"
+#include <type_traits>
+
+namespace aare::test_cuda {
+
+template <typename H>
+DeviceResult<H> run_gpu(NDView<float, 3> frames, size_t nbins, float xmin,
+                        float xmax, size_t batch_size) {
+    using D = std::conditional_t<sizeof(H) == 8, uint64_t, uint32_t>;
+    cuda::PixelHistogramCUDA<H, D> 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());
+    hist.fill_histogram(frames);
+    hist.unregister_input_buffer();
+    return {hist.bin_centers(), hist.bin_edges(), hist.values()};
+}
+
+#define INSTANTIATE(H)                                                         \
+    template DeviceResult<H> run_gpu<H>(NDView<float, 3>, size_t, float,       \
+                                        float, size_t);
+INSTANTIATE(uint8_t)
+INSTANTIATE(uint16_t)
+INSTANTIATE(uint32_t)
+INSTANTIATE(uint64_t)
+INSTANTIATE(float)
+INSTANTIATE(double)
+#undef INSTANTIATE
+
+} // namespace aare::test_cuda

src/hist/PixelHistogramCUDA_testimpl.hpp

@@ -0,0 +1,21 @@
+// SPDX-License-Identifier: MPL-2.0
+#pragma once
+#include "aare/NDArray.hpp"
+#include "aare/NDView.hpp"
+#include <cstddef>
+#include <cstdint>
+
+namespace aare::test_cuda {
+
+template <typename H> struct DeviceResult {
+    NDArray<float, 1> centers, edges;
+    NDArray<H, 3> values;
+};
+
+// Histogram `frames` (nframes, nrows, ncols) on the GPU. Explicitly
+// instantiated in the .cu for the storage types under test.
+template <typename H>
+DeviceResult<H> run_gpu(NDView<float, 3> frames, size_t nbins, float xmin,
+                        float xmax, size_t batch_size);
+
+} // namespace aare::test_cuda