More uniform binding of Numpy arrays (#381)

- Introduce a standard way to bind numpy arrays: py::array_t<T> and
.noconvert() on the argument
- Take py::array_t by value since it's just a small wrapper
- Update bindings

python/src/bind_Cluster.hpp

@@ -24,14 +24,21 @@ void define_Cluster(py::module &m, const std::string &typestr) {
     py::class_<Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>>(
         m, class_name.c_str(), py::buffer_protocol())
 
+        // conversion ok: small input, accept lists and any numeric dtype
         .def(py::init([](CoordType x, CoordType y,
                          py::array_t<Type, py::array::forcecast> data) {
-            py::buffer_info buf_info = data.request();
+            constexpr py::ssize_t n_pixels = ClusterSizeX * ClusterSizeY;
+            if (data.ndim() != 1 || data.size() != n_pixels) {
+                throw py::value_error(fmt::format(
+                    "Cluster data must be a 1D array of {} values, got {}D "
+                    "array with {} values",
+                    n_pixels, data.ndim(), data.size()));
+            }
             Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType> cluster;
             cluster.x = x;
             cluster.y = y;
             auto r = data.template unchecked<1>(); // no bounds checks
-            for (py::ssize_t i = 0; i < data.size(); ++i) {
+            for (py::ssize_t i = 0; i < n_pixels; ++i) {
                 cluster.data[i] = r(i);
             }
             return cluster;

python/src/bind_ClusterFile.hpp

@@ -138,9 +138,10 @@ void define_ClusterFile(py::module &m, const std::string &typestr) {
              "Estimate the number of clusters from the file size. Frame "
              "headers can make this larger than the actual count.")
         .def(
-            "set_noise_map",
+            "set_noise_map", // conversions ok
             [](ClusterFile<ClusterType> &self,
-               py::array_t<int32_t, py::array::c_style> noise_map) {
+               py::array_t<int32_t, py::array::c_style | py::array::forcecast>
+                   noise_map) {
                 auto view = make_view_2d(noise_map);
                 self.set_noise_map(view);
             },
@@ -149,9 +150,10 @@ void define_ClusterFile(py::module &m, const std::string &typestr) {
             "[y, x]. The map must cover every cluster center coordinate.")
 
         .def(
-            "set_gain_map",
+            "set_gain_map", // conversions ok
             [](ClusterFile<ClusterType> &self,
-               py::array_t<double, py::array::c_style> gain_map) {
+               py::array_t<double, py::array::c_style | py::array::forcecast>
+                   gain_map) {
                 auto view = make_view_2d(gain_map);
                 self.set_gain_map(view);
             },

python/src/bind_ClusterFinder.hpp

@@ -43,12 +43,14 @@ void define_ClusterFinder(py::module &m, const std::string &typestr) {
             &ClusterFinder<ClusterType, uint16_t, pd_type>::set_nSigma,
             R"(number of sigma above the pedestal to consider a photon during cluster finding.)")
 
-        .def("push_pedestal_frame",
-             [](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
-                py::array_t<uint16_t> frame) {
-                 auto view = make_view_2d(frame);
-                 self.push_pedestal_frame(view);
-             })
+        .def(
+            "push_pedestal_frame",
+            [](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
+               py::array_t<uint16_t> frame) {
+                auto view = make_view_2d(frame);
+                self.push_pedestal_frame(view);
+            },
+            py::arg("frame").noconvert())
         .def("clear_pedestal",
              &ClusterFinder<ClusterType, uint16_t, pd_type>::clear_pedestal)
         .def("update_threshold",
@@ -84,7 +86,7 @@ void define_ClusterFinder(py::module &m, const std::string &typestr) {
                 self.find_clusters(view, frame_number);
                 return;
             },
-            py::arg(), py::arg("frame_number") = 0);
+            py::arg("frame").noconvert(), py::arg("frame_number") = 0);
 }
 
 #pragma GCC diagnostic pop

python/src/bind_ClusterFinderMT.hpp

@@ -37,12 +37,14 @@ void define_ClusterFinderMT(py::module &m, const std::string &typestr) {
              py::arg("capacity") = 2048, py::arg("n_threads") = 3,
              py::arg("queue_depth") = 16,
              py::arg("min_pedestal_samples") = 1000)
-        .def("push_pedestal_frame",
-             [](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
-                py::array_t<uint16_t> frame) {
-                 auto view = make_view_2d(frame);
-                 self.push_pedestal_frame(view);
-             })
+        .def(
+            "push_pedestal_frame",
+            [](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
+               py::array_t<uint16_t> frame) {
+                auto view = make_view_2d(frame);
+                self.push_pedestal_frame(view);
+            },
+            py::arg("frame").noconvert())
         .def(
             "find_clusters",
             [](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
@@ -50,7 +52,7 @@ void define_ClusterFinderMT(py::module &m, const std::string &typestr) {
                 auto view = make_view_2d(frame);
                 self.find_clusters(view, frame_number);
             },
-            py::arg(), py::arg("frame_number") = 0)
+            py::arg("frame").noconvert(), py::arg("frame_number") = 0)
         .def_property_readonly(
             "cluster_size",
             [](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self) {

python/src/bind_ClusterVector.hpp

@@ -46,8 +46,7 @@ void define_ClusterVector(py::module &m, const std::string &typestr) {
 
         .def(
             "__call__",
-            [](ClusterVector<ClusterType> &self,
-               py::array_t<bool, py::array::c_style> mask) {
+            [](ClusterVector<ClusterType> &self, py::array_t<bool> mask) {
                 return self(make_view_1d(mask));
             },
             py::arg("mask").noconvert(), R"doc(

python/src/bind_FastPedestal.hpp

@@ -106,8 +106,7 @@ void define_fast_pedestal_bindings(py::module &m, const std::string &name) {
             "Return an independent copy of the pedestal and its state.")
         .def(
             "push_ema",
-            [](FastPedestal<SUM_TYPE> &pedestal,
-               py::array_t<uint16_t, py::array::c_style> &frame) {
+            [](FastPedestal<SUM_TYPE> &pedestal, py::array_t<uint16_t> frame) {
                 pedestal.push_ema(make_view_2d(frame));
             },
             py::arg("frame").noconvert(),
@@ -115,8 +114,7 @@ void define_fast_pedestal_bindings(py::module &m, const std::string &name) {
             "already be ready for this update.")
         .def(
             "add_init_frame",
-            [](FastPedestal<SUM_TYPE> &pedestal,
-               py::array_t<uint16_t, py::array::c_style> &frame) {
+            [](FastPedestal<SUM_TYPE> &pedestal, py::array_t<uint16_t> frame) {
                 pedestal.add_init_frame(make_view_2d(frame));
             },
             py::arg("frame").noconvert(),

python/src/bind_Interpolator.hpp

@@ -101,9 +101,7 @@ void define_interpolation_bindings(py::module &m) {
     auto interpolator =
         py::class_<aare::Interpolator>(m, "Interpolator")
             .def(py::init(
-                [](py::array_t<double,
-                               py::array::c_style | py::array::forcecast>
-                       etacube,
+                [](py::array_t<double> etacube,
                    py::array_t<double> xbins, py::array_t<double> ybins,
                    py::array_t<double> ebins) {
                     return Interpolator(
@@ -124,9 +122,9 @@ void define_interpolation_bindings(py::module &m) {
                 ebins: 
                     bin edges of photon energy
                 )doc",
-                py::arg("etacube"),
-                py::arg("xbins"), py::arg("ybins"),
-                py::arg("ebins"))
+                py::arg("etacube").noconvert(),
+                py::arg("xbins").noconvert(), py::arg("ybins").noconvert(),
+                py::arg("ebins").noconvert())
 
             .def(py::init(
                 [](py::array_t<double> xbins, py::array_t<double> ybins,
@@ -146,13 +144,12 @@ void define_interpolation_bindings(py::module &m) {
                     bin edges of etay
                 ebins: 
                     bin edges of photon energy
-                )", py::arg("xbins"),
-                py::arg("ybins"), py::arg("ebins"))
+                )", py::arg("xbins").noconvert(),
+                py::arg("ybins").noconvert(), py::arg("ebins").noconvert())
             .def(
                 "rosenblatttransform",
                 [](Interpolator &self,
-                   py::array_t<double,
-                               py::array::c_style | py::array::forcecast>
+                   py::array_t<double>
                        etacube) {
                     return self.rosenblatttransform(make_view_3d(etacube));
                 },
@@ -162,7 +159,7 @@ void define_interpolation_bindings(py::module &m) {
                 etacube: 
                     joint distribution of eta_x, eta_y and photon energy (**Note:** for the joint distribution first dimension is eta_x, second: eta_y, third: energy bins.)
                 )",
-                py::arg("etacube"))
+                py::arg("etacube").noconvert())
             .def("get_ietax",
                  [](Interpolator &self) {
                      auto *ptr = new NDArray<double, 3>{};

python/src/bind_Pedestal.hpp

@@ -85,11 +85,7 @@ void define_pedestal_bindings(py::module &m, const std::string &name) {
         // TODO! add push for other data types
         .def(
             "push",
-            [](Pedestal<SUM_TYPE> &pedestal,
-               py::array_t<uint16_t, py::array::c_style> &f) {
-                if (f.ndim() != 2) {
-                    throw py::value_error("Frame must be 2-dimensional");
-                }
+            [](Pedestal<SUM_TYPE> &pedestal, py::array_t<uint16_t> f) {
                 auto v = make_view_2d(f);
                 pedestal.push(v);
             },
@@ -99,15 +95,8 @@ void define_pedestal_bindings(py::module &m, const std::string &name) {
             "n_samples values per pixel, new values have weight 1 / n_samples.")
         .def(
             "push_with_threshold",
-            [](Pedestal<SUM_TYPE> &pedestal,
-               py::array_t<uint16_t, py::array::c_style> &f,
-               py::array_t<SUM_TYPE, py::array::c_style> &threshold) {
-                if (f.ndim() != 2) {
-                    throw py::value_error("Frame must be 2-dimensional");
-                }
-                if (threshold.ndim() != 2) {
-                    throw py::value_error("Threshold must be 2-dimensional");
-                }
+            [](Pedestal<SUM_TYPE> &pedestal, py::array_t<uint16_t> f,
+               py::array_t<SUM_TYPE> threshold) {
                 auto frame_view = make_view_2d(f);
                 auto threshold_view = make_view_2d(threshold);
                 pedestal.push_with_threshold(frame_view, threshold_view);

python/src/bind_PedestalTrackingPixelHistogram.hpp

@@ -50,8 +50,7 @@ void define_pedestal_tracking_pixel_histogram_bindings(py::module &m) {
         .def(
             "push_pedestal_no_update",
             [](PedestalTrackingPixelHistogram &self,
-               py::array_t<PedestalTrackingPixelHistogram::FrameType, 0>
-                   frame) {
+               py::array_t<PedestalTrackingPixelHistogram::FrameType> frame) {
                 auto view = make_view_2d(frame);
                 self.push_pedestal_no_update(view);
             },
@@ -104,8 +103,7 @@ void define_pedestal_tracking_pixel_histogram_bindings(py::module &m) {
         .def(
             "fill_async",
             [](PedestalTrackingPixelHistogram &self,
-               py::array_t<PedestalTrackingPixelHistogram::FrameType, 0>
-                   image) {
+               py::array_t<PedestalTrackingPixelHistogram::FrameType> image) {
                 // Copy the numpy buffer into an owned NDArray while we
                 // still hold the GIL so we don't depend on the array's
                 // backing storage outliving this call.

python/src/bind_PixelHistogram.hpp

@@ -46,7 +46,7 @@ void define_pixel_histogram_binding(py::module &m, const char *class_name,
 
         .def(
             "fill_async",
-            [](Hist &self, py::array_t<double, 0> image) {
+            [](Hist &self, py::array_t<double> image) {
                 // Copy the numpy buffer into an owned NDArray while we
                 // still hold the GIL so we don't depend on the array's
                 // backing storage outliving this call.

python/src/bind_calibration.hpp

@@ -11,38 +11,34 @@
 namespace py = pybind11;
 
 template <typename DataType>
-py::array_t<DataType> pybind_apply_calibration(
-    py::array_t<uint16_t, py::array::c_style | py::array::forcecast> data,
-    py::array_t<DataType, py::array::c_style | py::array::forcecast> pedestal,
-    py::array_t<DataType, py::array::c_style | py::array::forcecast>
-        calibration,
-    int n_threads = 4) {
+py::array_t<DataType>
+pybind_apply_calibration(py::array_t<uint16_t> data,
+                         py::array_t<DataType> pedestal,
+                         py::array_t<DataType> calibration, int n_threads = 4) {
 
     auto data_span = make_view_3d(data); // data is always 3D
     /* No pointer is passed, so NumPy will allocate the buffer */
     auto result = py::array_t<DataType>(data_span.shape());
     auto res = make_view_3d(result);
-    if (data.ndim() == 3 && pedestal.ndim() == 3 && calibration.ndim() == 3) {
+    if (pedestal.ndim() == 3 && calibration.ndim() == 3) {
         auto ped = make_view_3d(pedestal);
         auto cal = make_view_3d(calibration);
         aare::apply_calibration<DataType, 3>(res, data_span, ped, cal,
                                              n_threads);
-    } else if (data.ndim() == 3 && pedestal.ndim() == 2 &&
-               calibration.ndim() == 2) {
+    } else if (pedestal.ndim() == 2 && calibration.ndim() == 2) {
         auto ped = make_view_2d(pedestal);
         auto cal = make_view_2d(calibration);
         aare::apply_calibration<DataType, 2>(res, data_span, ped, cal,
                                              n_threads);
     } else {
         throw std::runtime_error(
-            "Invalid number of dimensions for data, pedestal or calibration");
+            "Invalid number of dimensions for pedestal or calibration");
     }
     return result;
 }
 
-py::array_t<int> pybind_count_switching_pixels(
-    py::array_t<uint16_t, py::array::c_style | py::array::forcecast> data,
-    ssize_t n_threads = 4) {
+py::array_t<int> pybind_count_switching_pixels(py::array_t<uint16_t> data,
+                                               ssize_t n_threads = 4) {
 
     auto data_span = make_view_3d(data);
     auto arr = new NDArray<int, 2>{};
@@ -51,9 +47,8 @@ py::array_t<int> pybind_count_switching_pixels(
 }
 
 template <typename T>
-py::array_t<T> pybind_calculate_pedestal(
-    py::array_t<uint16_t, py::array::c_style | py::array::forcecast> data,
-    ssize_t n_threads) {
+py::array_t<T> pybind_calculate_pedestal(py::array_t<uint16_t> data,
+                                         ssize_t n_threads) {
 
     auto data_span = make_view_3d(data);
     auto arr = new NDArray<T, 3>{};
@@ -62,9 +57,8 @@ py::array_t<T> pybind_calculate_pedestal(
 }
 
 template <typename T>
-py::array_t<T> pybind_calculate_pedestal_g0(
-    py::array_t<uint16_t, py::array::c_style | py::array::forcecast> data,
-    ssize_t n_threads) {
+py::array_t<T> pybind_calculate_pedestal_g0(py::array_t<uint16_t> data,
+                                            ssize_t n_threads) {
 
     auto data_span = make_view_3d(data);
     auto arr = new NDArray<T, 2>{};

python/src/ctb_raw_file.hpp

@@ -24,177 +24,119 @@
 namespace py = pybind11;
 using namespace ::aare;
 
+// One binding for the three ADC SAR decoders. Each is registered under
+// its original Python name so existing transforms keep working.
+template <void (*decode)(NDView<const uint8_t, 2>, NDView<uint16_t, 2>)>
+py::array_t<uint16_t> decode64to16(py::array_t<uint8_t> input) {
+    auto bytes = make_const_view_2d(input);
+    constexpr ssize_t bytes_per_word = sizeof(uint64_t);
+    auto output = py::array_t<uint16_t>(
+        {bytes.shape(0), bytes.shape(1) / bytes_per_word});
+    decode(bytes, make_view_2d(output));
+    return output;
+}
+
+constexpr auto decode64to16_doc =
+    "Decode packed 64-bit ADC samples. Takes a two-dimensional, "
+    "C-contiguous uint8 array whose row length is a multiple of 8 bytes "
+    "and returns a uint16 array with one value per 64-bit word.";
+
 void define_ctb_raw_file_io_bindings(py::module &m) {
 
-    m.def("adc_sar_05_06_07_08decode64to16", [](py::array_t<uint8_t> input) {
-        if (input.ndim() != 2) {
-            throw std::runtime_error(
-                "Only 2D arrays are supported at this moment");
-        }
+    m.def("adc_sar_05_06_07_08decode64to16",
+          &decode64to16<adc_sar_05_06_07_08decode64to16>,
+          py::arg("input").noconvert(), decode64to16_doc);
 
-        // Create a 2D output array with the same shape as the input
-        std::vector<ssize_t> shape{input.shape(0),
-                                   input.shape(1) /
-                                       static_cast<ssize_t>(bits_per_byte)};
-        py::array_t<uint16_t> output(shape);
+    m.def("adc_sar_05_decode64to16", &decode64to16<adc_sar_05_decode64to16>,
+          py::arg("input").noconvert(), decode64to16_doc);
 
-        // Create a view of the input and output arrays
-        NDView<uint64_t, 2> input_view(
-            reinterpret_cast<uint64_t *>(input.mutable_data()),
-            {output.shape(0), output.shape(1)});
-        NDView<uint16_t, 2> output_view(output.mutable_data(),
-                                        {output.shape(0), output.shape(1)});
+    m.def("adc_sar_04_decode64to16", &decode64to16<adc_sar_04_decode64to16>,
+          py::arg("input").noconvert(), decode64to16_doc);
 
-        adc_sar_05_06_07_08decode64to16(input_view, output_view);
+    m.def(
+        "apply_custom_weights",
+        [](py::array_t<uint16_t> input, py::array_t<double> weights) {
+            auto input_view = make_view_1d(input);
+            auto weights_view = make_view_1d(weights);
+            auto output = py::array_t<double>(input_view.size());
+            apply_custom_weights(input_view, make_view_1d(output),
+                                 weights_view);
+            return output;
+        },
+        py::arg("input").noconvert(), py::arg("weights").noconvert(),
+        "Apply per-bit weights to every value of a one-dimensional, "
+        "C-contiguous uint16 array and return a float64 array.");
 
-        return output;
-    });
+    m.def(
+        "expand24to32bit",
+        [](py::array_t<uint8_t> input, uint32_t offset) {
+            constexpr ssize_t bytes_per_channel = 3; // 24 bit
+            auto input_view = make_const_view_1d(input);
+            auto output =
+                py::array_t<uint32_t>(input_view.size() / bytes_per_channel);
+            aare::expand24to32bit(input_view, make_view_1d(output),
+                                  aare::BitOffset(offset));
+            return output;
+        },
+        py::arg("input").noconvert(), py::arg("offset"),
+        "Expand packed 24-bit values from a one-dimensional, C-contiguous "
+        "uint8 array into a uint32 array.");
 
-    m.def("adc_sar_05_decode64to16", [](py::array_t<uint8_t> input) {
-        if (input.ndim() != 2) {
-            throw std::runtime_error(
-                "Only 2D arrays are supported at this moment");
-        }
+    m.def(
+        "expand4to8bit",
+        [](py::array_t<uint8_t> input) {
+            auto input_view = make_const_view_1d(input);
+            auto output = py::array_t<uint8_t>(input_view.size() * 2);
+            aare::expand4to8bit(input_view, make_view_1d(output));
+            return output;
+        },
+        py::arg("input").noconvert(),
+        "Expand two 4-bit values per byte of a one-dimensional, C-contiguous "
+        "uint8 array into a uint8 array of twice the length.");
 
-        // Create a 2D output array with the same shape as the input
-        std::vector<ssize_t> shape{input.shape(0),
-                                   input.shape(1) /
-                                       static_cast<ssize_t>(bits_per_byte)};
-        py::array_t<uint16_t> output(shape);
+    m.def(
+        "decode_my302",
+        [](py::array_t<uint8_t> input, uint32_t offset) {
+            // Physical layout of the chip
+            constexpr ssize_t channels = 64;
+            constexpr ssize_t counters = 3;
+            constexpr ssize_t bytes_per_channel = 3; // 24 bit
+            constexpr ssize_t n_outputs = 2;
 
-        // Create a view of the input and output arrays
-        NDView<uint64_t, 2> input_view(
-            reinterpret_cast<uint64_t *>(input.mutable_data()),
-            {output.shape(0), output.shape(1)});
-        NDView<uint16_t, 2> output_view(output.mutable_data(),
-                                        {output.shape(0), output.shape(1)});
+            auto input_view = make_const_view_1d(input);
 
-        adc_sar_05_decode64to16(input_view, output_view);
+            ssize_t expected_size = channels * counters * bytes_per_channel;
 
-        return output;
-    });
+            // If we have an offset we need one extra byte per output
+            aare::BitOffset bitoff(offset);
+            if (bitoff.value())
+                expected_size += n_outputs;
 
-    m.def("adc_sar_04_decode64to16", [](py::array_t<uint8_t> input) {
-        if (input.ndim() != 2) {
-            throw std::runtime_error(
-                "Only 2D arrays are supported at this moment");
-        }
+            if (input_view.size() != expected_size) {
+                throw py::value_error(
+                    fmt::format("{} Expected an input size of {} bytes. Called "
+                                "with input size of {}",
+                                LOCATION, expected_size, input_view.size()));
+            }
 
-        // Create a 2D output array with the same shape as the input
-        std::vector<ssize_t> shape{input.shape(0),
-                                   input.shape(1) /
-                                       static_cast<ssize_t>(bits_per_byte)};
-        py::array_t<uint16_t> output(shape);
+            auto output = py::array_t<uint32_t>(channels * counters);
+            auto output_view = make_view_1d(output);
 
-        // Create a view of the input and output arrays
-        NDView<uint64_t, 2> input_view(
-            reinterpret_cast<uint64_t *>(input.mutable_data()),
-            {output.shape(0), output.shape(1)});
-        NDView<uint16_t, 2> output_view(output.mutable_data(),
-                                        {output.shape(0), output.shape(1)});
+            const ssize_t in_step = input_view.size() / n_outputs;
+            const ssize_t out_step = output_view.size() / n_outputs;
+            for (ssize_t i = 0; i != n_outputs; ++i) {
+                NDView<const uint8_t, 1> in_part(
+                    input_view.data() + in_step * i, {in_step});
+                NDView<uint32_t, 1> out_part(output_view.data() + out_step * i,
+                                             {out_step});
+                aare::expand24to32bit(in_part, out_part, bitoff);
+            }
 
-        adc_sar_04_decode64to16(input_view, output_view);
-
-        return output;
-    });
-
-    m.def("apply_custom_weights",
-          [](py::array_t<uint16_t, py::array::c_style | py::array::forcecast>
-                 &input,
-             py::array_t<double, py::array::c_style | py::array::forcecast>
-                 &weights) {
-              // Create new array with same shape as the input array
-              // (uninitialized values)
-              py::buffer_info buf = input.request();
-              py::array_t<double> output(buf.shape);
-
-              // Use NDViews to call into the C++ library
-              auto weights_view = make_view_1d(weights);
-              NDView<uint16_t, 1> input_view(input.mutable_data(),
-                                             {input.size()});
-              NDView<double, 1> output_view(output.mutable_data(),
-                                            {output.size()});
-
-              apply_custom_weights(input_view, output_view, weights_view);
-              return output;
-          });
-
-    m.def("expand24to32bit",
-          [](py::array_t<uint8_t, py::array::c_style | py::array::forcecast>
-                 &input,
-             uint32_t offset) {
-              aare::BitOffset bitoff(offset);
-              py::buffer_info buf = input.request();
-
-              constexpr uint32_t bytes_per_channel = 3; // 24 bit
-              py::array_t<uint32_t> output(buf.size / bytes_per_channel);
-
-              NDView<uint8_t, 1> input_view(input.mutable_data(),
-                                            {input.size()});
-              NDView<uint32_t, 1> output_view(output.mutable_data(),
-                                              {output.size()});
-
-              aare::expand24to32bit(input_view, output_view, bitoff);
-              return output;
-          });
-
-    m.def("expand4to8bit",
-          [](py::array_t<uint8_t, py::array::c_style | py::array::forcecast>
-                 &input) {
-              py::buffer_info buf = input.request();
-
-              py::array_t<uint8_t> output(buf.size * 2);
-
-              NDView<uint8_t, 1> input_view(input.mutable_data(),
-                                            {input.size()});
-              NDView<uint8_t, 1> output_view(output.mutable_data(),
-                                             {output.size()});
-
-              aare::expand4to8bit(input_view, output_view);
-              return output;
-          });
-
-    m.def("decode_my302",
-          [](py::array_t<uint8_t, py::array::c_style | py::array::forcecast>
-                 &input,
-             uint32_t offset) {
-              // Physical layout of the chip
-              constexpr size_t channels = 64;
-              constexpr size_t counters = 3;
-              constexpr size_t bytes_per_channel = 3; // 24 bit
-              constexpr int n_outputs = 2;
-
-              ssize_t expected_size = channels * counters * bytes_per_channel;
-
-              // If whe have an offset we need one extra byte per output
-              aare::BitOffset bitoff(offset);
-              if (bitoff.value())
-                  expected_size += n_outputs;
-
-              if (input.size() != expected_size) {
-                  throw std::runtime_error(fmt::format(
-                      "{} Expected an input size of {} bytes. Called "
-                      "with input size of {}",
-                      LOCATION, expected_size, input.size()));
-              }
-
-              py::buffer_info buf = input.request();
-              py::array_t<uint32_t> output(channels * counters);
-
-              for (int i = 0; i != n_outputs; ++i) {
-                  auto step = input.size() / n_outputs;
-                  auto out_step = output.size() / n_outputs;
-                  NDView<uint8_t, 1> input_view(input.mutable_data() + step * i,
-                                                {input.size() / n_outputs});
-                  NDView<uint32_t, 1> output_view(output.mutable_data() +
-                                                      out_step * i,
-                                                  {output.size() / n_outputs});
-
-                  aare::expand24to32bit(input_view, output_view, bitoff);
-              }
-
-              return output;
-          });
+            return output;
+        },
+        py::arg("input").noconvert(), py::arg("offset"),
+        "Decode a Mythen 302 readout from a one-dimensional, C-contiguous "
+        "uint8 array into a uint32 array of 64 channels times 3 counters.");
 
     py::class_<CtbRawFile>(m, "CtbRawFile")
         .def(py::init<const std::filesystem::path &>())

python/src/fit.hpp

@@ -60,7 +60,8 @@ template <typename Model> void bind_fit_model(py::module &m, const char *name) {
         .def_property("compute_errors", &FM::compute_errors,
                       &FM::SetComputeErrors)
         .def(
-            "__call__",
+            "__call__", // conversion ok, we want to be able to call with any
+                        // dtype
             [](const FM & /*self*/,
                py::array_t<double, py::array::c_style | py::array::forcecast> x,
                py::array_t<double, py::array::c_style | py::array::forcecast>
@@ -78,7 +79,7 @@ template <typename Model> void bind_fit_model(py::module &m, const char *name) {
             },
             py::arg("x"), py::arg("par"))
         .def(
-            "fit",
+            "fit", // conversion ok
             [](const FM &self,
                py::array_t<double, py::array::c_style | py::array::forcecast> x,
                py::array_t<double, py::array::c_style | py::array::forcecast> y,
@@ -163,11 +164,6 @@ fit_dispatch(const aare::FitModel<Model> &model,
                 py::array_t<double, py::array::c_style | py::array::forcecast>>(
                 y_err_obj);
 
-            if (y_err.ndim() != 3) {
-                throw std::runtime_error(
-                    "For 3D input y, y_err must also be 3D.");
-            }
-
             auto err_out =
                 new NDArray<double, 3>({y.shape(0), y.shape(1), npar}, 0.0);
             auto y_view_err = make_view_3d(y_err);
@@ -208,11 +204,6 @@ fit_dispatch(const aare::FitModel<Model> &model,
                 py::array_t<double, py::array::c_style | py::array::forcecast>>(
                 y_err_obj);
 
-            if (y_err.ndim() != 1) {
-                throw std::runtime_error(
-                    "For 1D input y, y_err must also be 1D.");
-            }
-
             auto y_view_err = make_view_1d(y_err);
             result = aare::fit_pixel<Model>(model, x_view, y_view, y_view_err);
         } else {
@@ -241,7 +232,8 @@ void define_fit_bindings(py::module &m) {
 
     m.def(
         "fit",
-        [](py::object model_obj,
+        [](py::object model_obj, // conversion ok, we want to be able to call
+                                 // with any dtype
            py::array_t<double, py::array::c_style | py::array::forcecast> x,
            py::array_t<double, py::array::c_style | py::array::forcecast> y,
            py::object y_err_obj, int n_threads) -> py::object {

python/src/var_cluster.hpp

@@ -21,40 +21,42 @@ void define_var_cluster_finder_bindings(py::module &m) {
                  auto *ptr = new NDArray<int, 2>(self.labeled());
                  return return_image_data(ptr);
              })
-        .def("set_noiseMap",
-             [](VarClusterFinder<double> &self,
-                py::array_t<double, py::array::c_style | py::array::forcecast>
-                    noise_map) {
-                 auto noise_map_span = make_view_2d(noise_map);
-                 self.set_noiseMap(noise_map_span);
-             })
+        .def(
+            "set_noiseMap", // conversions ok
+            [](VarClusterFinder<double> &self,
+               py::array_t<double, py::array::c_style | py::array::forcecast>
+                   noise_map) {
+                auto noise_map_span = make_view_2d(noise_map);
+                self.set_noiseMap(noise_map_span);
+            },
+            py::arg("noise_map"))
         .def("set_numberOfNeighbours",
              &VarClusterFinder<double>::set_numberOfNeighbours)
         .def("set_empty_surroundingPixels",
              &VarClusterFinder<double>::set_empty_surroundingPixels)
         .def("set_peripheralThresholdFactor",
              &VarClusterFinder<double>::set_peripheralThresholdFactor)
-        .def("find_clusters",
-             [](VarClusterFinder<double> &self,
-                py::array_t<double, py::array::c_style | py::array::forcecast>
-                    img) {
-                 auto view = make_view_2d(img);
-                 self.find_clusters(view);
-             })
-        .def("find_clusters_X",
-             [](VarClusterFinder<double> &self,
-                py::array_t<double, py::array::c_style | py::array::forcecast>
-                    img) {
-                 auto img_span = make_view_2d(img);
-                 self.find_clusters_X(img_span);
-             })
-        .def("single_pass",
-             [](VarClusterFinder<double> &self,
-                py::array_t<double, py::array::c_style | py::array::forcecast>
-                    img) {
-                 auto img_span = make_view_2d(img);
-                 self.single_pass(img_span);
-             })
+        .def(
+            "find_clusters",
+            [](VarClusterFinder<double> &self, py::array_t<double> img) {
+                auto view = make_view_2d(img);
+                self.find_clusters(view);
+            },
+            py::arg("img").noconvert())
+        .def(
+            "find_clusters_X",
+            [](VarClusterFinder<double> &self, py::array_t<double> img) {
+                auto img_span = make_view_2d(img);
+                self.find_clusters_X(img_span);
+            },
+            py::arg("img").noconvert())
+        .def(
+            "single_pass",
+            [](VarClusterFinder<double> &self, py::array_t<double> img) {
+                auto img_span = make_view_2d(img);
+                self.single_pass(img_span);
+            },
+            py::arg("img").noconvert())
         .def("hits",
              [](VarClusterFinder<double> &self) {
                  auto ptr = new std::vector<VarClusterFinder<double>::Hit>(

python/tests/test_Cluster.py

@@ -162,3 +162,21 @@ def test_3x3_reduction():
 
 
 
+
+
+@pytest.mark.parametrize(
+    "cluster_type, n_pixels",
+    [(_aare.Cluster3x3i, 9), (_aare.Cluster2x2i, 4)],
+)
+def test_cluster_constructor_requires_exact_data_length(cluster_type, n_pixels):
+    """Data must hold exactly one value per pixel"""
+    for bad in (
+        np.ones(n_pixels - 1, dtype=np.int32),
+        np.ones(n_pixels + 1, dtype=np.int32),
+        np.ones((n_pixels, 1), dtype=np.int32),
+    ):
+        with pytest.raises(ValueError, match=f"{n_pixels} values"):
+            cluster_type(0, 0, bad)
+
+    cluster = cluster_type(0, 0, np.arange(n_pixels, dtype=np.int32))
+    assert (cluster.data.ravel() == np.arange(n_pixels)).all()

python/tests/test_ClusterVector.py

@@ -146,7 +146,7 @@ def test_masking_requires_c_contiguous_array():
     mask = np.array([True, False, True, False], dtype=bool)[::2]
     assert not mask.flags.c_contiguous
 
-    with pytest.raises(TypeError):
+    with pytest.raises(ValueError, match="C-contiguous"):
         cv(mask)
 
 

python/tests/test_Pedestal.py

@@ -172,7 +172,7 @@ def test_pedestal_rejects_noncontiguous_inputs(dtype, input_name, layout):
         invalid = np.ones((4, 6), dtype=input_dtype)[:2, :3][:, ::-1]
     assert not invalid.flags.c_contiguous
 
-    with pytest.raises(TypeError):
+    with pytest.raises(ValueError, match="C-contiguous"):
         if input_name == "push_frame":
             pedestal.push(invalid)
         elif input_name == "threshold_frame":
@@ -192,9 +192,8 @@ def test_pedestal_rejects_inputs_with_wrong_ndim(dtype, input_name, shape):
     threshold = np.full((2, 3), 10, dtype=dtype)
     input_dtype = dtype if input_name == "threshold" else np.uint16
     invalid = np.ones(shape, dtype=input_dtype)
-    name = "Threshold" if input_name == "threshold" else "Frame"
 
-    with pytest.raises(ValueError, match=f"{name} must be 2-dimensional"):
+    with pytest.raises(ValueError):
         if input_name == "push_frame":
             pedestal.push(invalid)
         elif input_name == "threshold_frame":