PixelHistogramCUDA: take frames as InputType, keep a floating-point axis

The third template parameter is now the frame type copied to the device,
so raw uint16 frames can be histogrammed without a host-side conversion.
The axis (range, edges, bin arithmetic) is AxisType: InputType for
floating-point input, float for integers of at most 16 bits; wider
integers are rejected at compile time. The kernel converts each value.

Test: uint16 frames on 1.2-ADU bins match PixelHistogram fed the same
values as float.

include/aare/hist/PixelHistogramCUDA.hpp

@@ -21,7 +21,7 @@
 namespace aare::cuda {
 
 template <typename StorageType = uint16_t,
-          typename DeviceStorageType = uint32_t, typename AxisType = float,
+          typename DeviceStorageType = uint32_t, typename InputType = float,
           int BLOCK_SIZE = 128>
 class PixelHistogramCUDA {
 #ifdef AARE_TESTS
@@ -30,6 +30,20 @@ class PixelHistogramCUDA {
     FRIEND_TEST(known_pattern_transpose_and_cap)
 #endif
 
+    // Frames are copied to the device as InputType, as the caller provides
+    // them. The histogram axis (range, bin edges, the bin arithmetic) must be
+    // floating point: an integral input gets a float axis, converted in the
+    // kernel. That conversion is exact up to 16-bit integers.
+    static_assert(std::is_floating_point_v<InputType> ||
+                      (std::is_integral_v<InputType> && sizeof(InputType) <= 2),
+                  "PixelHistogramCUDA: InputType must be floating point or an "
+                  "integer of at most 16 bits");
+
+  public:
+    using AxisType = std::conditional_t<std::is_floating_point_v<InputType>,
+                                        InputType, float>;
+
+  private:
     // nrows is declared first: its initializer validates all constructor
     // arguments, so it must run before anything below allocates.
     size_t nrows;
@@ -51,7 +65,7 @@ class PixelHistogramCUDA {
 
     // Device buffers for frames and histogram
     DeviceBuffer<DeviceStorageType> d_histogram;
-    std::array<DeviceBuffer<AxisType>, NUM_SLOTS> d_frames_batch;
+    std::array<DeviceBuffer<InputType>, NUM_SLOTS> d_frames_batch;
 
     // Events for synchronization
     std::array<Event, NUM_SLOTS> copy_done;
@@ -120,7 +134,7 @@ class PixelHistogramCUDA {
         // Initialize device buffers for frames
         for (int i = 0; i < NUM_SLOTS; ++i) {
             d_frames_batch[i] =
-                DeviceBuffer<AxisType>(nrows * ncols * m_batch_size);
+                DeviceBuffer<InputType>(nrows * ncols * m_batch_size);
         }
     }
 
@@ -141,7 +155,7 @@ class PixelHistogramCUDA {
         m_registered_input.unregister();
     }
 
-    void submit_batch(NDView<AxisType, 3> frames) {
+    void submit_batch(NDView<InputType, 3> frames) {
         auto frame_shape = frames.shape();
         if (static_cast<size_t>(frame_shape[0]) > m_batch_size ||
             static_cast<size_t>(frame_shape[1]) != nrows ||
@@ -156,7 +170,7 @@ class PixelHistogramCUDA {
         const int slot = m_next_batch++ % NUM_SLOTS;
 
         // pointer to the first frame in the batch
-        const AxisType *h_frames_batch = frames.data();
+        const InputType *h_frames_batch = frames.data();
 
         // Make sure `copy_stream` waits on any kernel (using `slot`) to finish
         // execution
@@ -166,7 +180,7 @@ class PixelHistogramCUDA {
         // Copy frames to device
         const size_t npix = nrows * ncols;
         CUDA_CHECK(cudaMemcpyAsync(d_frames_batch[slot].get(), h_frames_batch,
-                                   npix * batch_size * sizeof(AxisType),
+                                   npix * batch_size * sizeof(InputType),
                                    cudaMemcpyHostToDevice, copy_stream.get()));
 
         // Record the copy_done event (in the copy_stream) right after the copy
@@ -191,11 +205,11 @@ class PixelHistogramCUDA {
             cudaEventRecord(kernel_done[slot].get(), compute_stream.get()));
     }
 
-    void fill_histogram(NDView<AxisType, 3> frames) {
+    void fill_histogram(NDView<InputType, 3> frames) {
         auto N_frames = static_cast<size_t>(frames.shape(0));
         for (size_t i = 0; i < N_frames; i += m_batch_size) {
             size_t batch_end = std::min(i + m_batch_size, N_frames);
-            NDView<AxisType, 3> batch = frames.sub_view(i, batch_end);
+            NDView<InputType, 3> batch = frames.sub_view(i, batch_end);
             submit_batch(batch);
         }
     }

include/aare/hist/pixelhistogram_kernel.cuh

@@ -33,8 +33,10 @@ bin_to_pixel_major(const DeviceStorageType *__restrict__ d_histogram_binMajor,
         count < saturation_value ? count : saturation_value);
 }
 
-template <typename DeviceStorageType = uint32_t, typename AxisType = float>
-__global__ void fill_bin_major(const AxisType *__restrict__ d_frames_batch,
+// Frames arrive as InputType; each value is converted to the floating-point
+// AxisType before binning.
+template <typename DeviceStorageType, typename InputType, typename AxisType>
+__global__ void fill_bin_major(const InputType *__restrict__ d_frames_batch,
                                DeviceStorageType *__restrict__ d_histogram,
                                const size_t npix, const size_t nbins,
                                AxisType xmin, AxisType xmax, AxisType scale,
@@ -46,7 +48,8 @@ __global__ void fill_bin_major(const AxisType *__restrict__ d_frames_batch,
 
     for (size_t fid = 0; fid < batch_size; ++fid) {
 
-        AxisType val_pix_at_fid = d_frames_batch[fid * npix + pix_id];
+        const AxisType val_pix_at_fid =
+            static_cast<AxisType>(d_frames_batch[fid * npix + pix_id]);
         if (!(val_pix_at_fid >= xmin &&
               val_pix_at_fid < xmax)) // also skips NaN
             continue;

src/hist/PixelHistogramCUDA.test.cpp

@@ -160,6 +160,27 @@ TEST_CASE("PixelHistogramCUDA result does not depend on the batch size",
     REQUIRE((cpu.values() == gpu.values));
 }
 
+TEST_CASE("PixelHistogramCUDA bins uint16 frames on a float axis", "[cuda]") {
+    // 1000 bins over 1200 ADU: bins 1.2 wide, so an integral axis (scale
+    // 1000 / 1200 == 0 in integers) would put everything in bin 0.
+    const size_t nrows = 7, ncols = 13, nframes = 100, nbins = 1000;
+    const float xmin = 3900.0f, xmax = 5100.0f;
+    auto frames = make_synth(nrows, ncols, nframes); // integer-valued floats
+
+    NDArray<uint16_t, 3> raw(frames.shape());
+    for (ssize_t i = 0; i < frames.size(); ++i)
+        raw.data()[i] = static_cast<uint16_t>(frames.data()[i]);
+
+    // The CPU reference gets the same values as float.
+    PixelHistogram<uint32_t> cpu(nrows, ncols, nbins, xmin, xmax);
+    fill_cpu(cpu, frames.view());
+
+    auto gpu = test_cuda::run_gpu<uint32_t>(raw.view(), nbins, xmin, xmax,
+                                            /*batch_size=*/8);
+    REQUIRE((cpu.bin_edges() == gpu.edges));
+    REQUIRE((cpu.values() == gpu.values));
+}
+
 TEST_CASE("PixelHistogramCUDA accumulates over several fills with values() "
           "in between",
           "[cuda]") {

src/hist/PixelHistogramCUDA_testimpl.cu

@@ -67,28 +67,29 @@ 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,
+template <typename H, typename In>
+DeviceResult<H> run_gpu(NDView<In, 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);
+    cuda::PixelHistogramCUDA<H, D, In> 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)
+#define INSTANTIATE(H, In)                                                     \
+    template DeviceResult<H> run_gpu<H, In>(NDView<In, 3>, size_t, float,      \
+                                            float, size_t);
+INSTANTIATE(uint8_t, float)
+INSTANTIATE(uint16_t, float)
+INSTANTIATE(uint32_t, float)
+INSTANTIATE(uint64_t, float)
+INSTANTIATE(float, float)
+INSTANTIATE(double, float)
+INSTANTIATE(uint32_t, uint16_t)
 #undef INSTANTIATE
 
 std::vector<NDArray<uint32_t, 3>>

src/hist/PixelHistogramCUDA_testimpl.hpp

@@ -14,10 +14,11 @@ template <typename H> struct DeviceResult {
     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,
+// Histogram `frames` (nframes, nrows, ncols) on the GPU, with frames of type
+// In. Explicitly instantiated in the .cu for the storage and input types
+// under test; the axis is float for all of them.
+template <typename H, typename In>
+DeviceResult<H> run_gpu(NDView<In, 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