Merge branch 'main' into dev/strixels/remap_simple

python/aare/Cluster.py

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 from . import _aare 
 import numpy as np
 from .factory import _type_to_char

python/aare/__init__.py

@@ -22,9 +22,10 @@ from ._aare import (
     RawSubFile,
 )
 from ._aare import Pedestal_d, Pedestal_f, ClusterFinder_Cluster3x3i, VarClusterFinder
-from ._aare import DetectorType, ReadoutMode 
+from ._aare import DetectorType, FrameDiscardPolicy, ReadoutMode
 from ._aare import hitmap
 from ._aare import ROI
+from ._aare import xy
 from ._aare import corner 
 
 from ._aare import UDPPortPosition
@@ -40,6 +41,7 @@ from .Cluster import Cluster
 
 from ._aare import Gaussian, RisingScurve, FallingScurve, Pol1, Pol2, GaussianErfcPlateau, GaussianChargeSharing, GaussianChargeSharingKb
 from ._aare import fit
+from ._aare import Minimizer
 from ._aare import Interpolator
 from ._aare import calculate_eta2, calculate_eta3, calculate_cross_eta3, calculate_full_eta2
 from ._aare import reduce_to_2x2, reduce_to_3x3

python/aare/_version.py

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 # This file is used to get the version of the package from the VERSION file
 from importlib.metadata import PackageNotFoundError, version
 from pathlib import Path

python/aare/calibration.py

@@ -5,9 +5,12 @@ def load_calibration(fname, hg0=False):
     """
     Load calibration data from a file.
     
-    Parameters:
-    fname (str): Path to the calibration file.
-    hg0 (bool): If True, load HG0 calibration data instead of G0.
+    Parameters
+    ----------
+    fname : str
+        Path to the calibration file.
+    hg0 : bool
+        If True, load HG0 calibration data instead of G0.
 
     """
     gains = 3

python/src/bind_Cluster.hpp

@@ -13,9 +13,6 @@ namespace py = pybind11;
 
 using namespace aare;
 
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wunused-parameter"
-
 template <typename Type, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
           typename CoordType>
 void define_Cluster(py::module &m, const std::string &typestr) {
@@ -24,14 +21,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;
@@ -133,5 +137,3 @@ void reduce_to_2x2(py::module &m) {
         row-major order.
         )doc");
 }
-
-#pragma GCC diagnostic pop

python/src/bind_ClusterCollector.hpp

@@ -20,9 +20,6 @@ namespace py = pybind11;
 
 using namespace aare;
 
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wunused-parameter"
-
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
           typename CoordType = uint16_t>
 void define_ClusterCollector(py::module &m, const std::string &typestr) {
@@ -42,5 +39,3 @@ void define_ClusterCollector(py::module &m, const std::string &typestr) {
             },
             py::return_value_policy::take_ownership);
 }
-
-#pragma GCC diagnostic pop

python/src/bind_ClusterFile.hpp

@@ -13,11 +13,6 @@
 #include <string>
 #include <utility>
 
-// Disable warnings for unused parameters, as we ignore some
-// in the __exit__ method
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wunused-parameter"
-
 namespace py = pybind11;
 using namespace ::aare;
 
@@ -138,9 +133,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 +145,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);
             },
@@ -169,14 +166,12 @@ void define_ClusterFile(py::module &m, const std::string &typestr) {
         .def("__enter__", [](ClusterFile<ClusterType> &self) { return &self; })
         .def("__exit__",
              [](ClusterFile<ClusterType> &self,
-                const std::optional<pybind11::type> &exc_type,
-                const std::optional<pybind11::object> &exc_value,
-                const std::optional<pybind11::object> &traceback) {
-                 self.close();
-             })
+                [[maybe_unused]] const std::optional<pybind11::type> &exc_type,
+                [[maybe_unused]] const std::optional<pybind11::object>
+                    &exc_value,
+                [[maybe_unused]] const std::optional<pybind11::object>
+                    &traceback) { self.close(); })
         .def("__iter__", [](ClusterFile<ClusterType> &self) { return &self; })
         .def("__next__",
              [](File &self) { return ChunkIterator(self.chunks()).next(); });
 }
-
-#pragma GCC diagnostic pop

python/src/bind_ClusterFileSink.hpp

@@ -17,9 +17,6 @@
 namespace py = pybind11;
 using namespace aare;
 
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wunused-parameter"
-
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
           typename CoordType = uint16_t>
 void define_ClusterFileSink(py::module &m, const std::string &typestr) {
@@ -34,5 +31,3 @@ void define_ClusterFileSink(py::module &m, const std::string &typestr) {
                       const std::filesystem::path &>())
         .def("stop", &ClusterFileSink<ClusterType>::stop);
 }
-
-#pragma GCC diagnostic pop

python/src/bind_ClusterFinder.hpp

@@ -20,9 +20,6 @@ namespace py = pybind11;
 
 using namespace aare;
 
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wunused-parameter"
-
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
           typename CoordType = uint16_t>
 void define_ClusterFinder(py::module &m, const std::string &typestr) {
@@ -43,16 +40,22 @@ 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",
              &ClusterFinder<ClusterType, uint16_t, pd_type>::update_threshold)
+        .def_property_readonly(
+            "pedestal_ready",
+            &ClusterFinder<ClusterType, uint16_t, pd_type>::pedestal_ready,
+            R"(True once enough pedestal frames have been pushed to find clusters.)")
         .def_property_readonly(
             "pedestal",
             [](ClusterFinder<ClusterType, uint16_t, pd_type> &self) {
@@ -84,7 +87,5 @@ 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

@@ -20,9 +20,6 @@ namespace py = pybind11;
 
 using namespace aare;
 
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wunused-parameter"
-
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
           typename CoordType = uint16_t>
 void define_ClusterFinderMT(py::module &m, const std::string &typestr) {
@@ -37,12 +34,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,16 +49,21 @@ 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) {
+            []([[maybe_unused]] ClusterFinderMT<ClusterType, uint16_t, pd_type>
+                   &self) {
                 return py::make_tuple(ClusterSizeX, ClusterSizeY);
             })
         .def("clear_pedestal",
              &ClusterFinderMT<ClusterType, uint16_t, pd_type>::clear_pedestal)
         .def("update_threshold",
              &ClusterFinderMT<ClusterType, uint16_t, pd_type>::update_threshold)
+        .def(
+            "pedestal_ready",
+            &ClusterFinderMT<ClusterType, uint16_t, pd_type>::pedestal_ready,
+            R"(waits for the queued frames, then returns True if all cluster finders have enough pedestal frames to find clusters.)")
         .def("sync", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::sync)
         .def("stop", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::stop)
         .def("start", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::start)
@@ -87,5 +91,3 @@ void define_ClusterFinderMT(py::module &m, const std::string &typestr) {
              py::arg("nSigma"),
              R"(sets the number of sigma for all cluster finders.)");
 }
-
-#pragma GCC diagnostic pop

python/src/bind_ClusterVector.hpp

@@ -18,9 +18,6 @@ namespace py = pybind11;
 
 using namespace aare;
 
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wunused-parameter"
-
 template <typename Type, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
           typename CoordType = uint16_t>
 void define_ClusterVector(py::module &m, const std::string &typestr) {
@@ -46,11 +43,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) {
-                if (mask.ndim() != 1) {
-                    throw py::value_error("Mask must be one-dimensional");
-                }
+            [](ClusterVector<ClusterType> &self, py::array_t<bool> mask) {
                 return self(make_view_1d(mask));
             },
             py::arg("mask").noconvert(), R"doc(
@@ -127,8 +120,8 @@ void define_ClusterVector(py::module &m, const std::string &typestr) {
              "padding.")
         .def_property_readonly(
             "fmt",
-            [typestr](ClusterVector<ClusterType> &self) {
-                return fmt_format<ClusterType>;
+            [typestr]([[maybe_unused]] ClusterVector<ClusterType> &self) {
+                return fmt_format<ClusterType>();
             },
             "PEP 3118 format string for one stored cluster.")
 
@@ -149,11 +142,11 @@ void define_ClusterVector(py::module &m, const std::string &typestr) {
         .def_buffer(
             [typestr](ClusterVector<ClusterType> &self) -> py::buffer_info {
                 return py::buffer_info(
-                    self.data(),             /* Pointer to buffer */
-                    self.item_size(),        /* Size of one scalar */
-                    fmt_format<ClusterType>, /* Format descriptor */
-                    1,                       /* Number of dimensions */
-                    {self.size()},           /* Buffer dimensions */
+                    self.data(),               /* Pointer to buffer */
+                    self.item_size(),          /* Size of one scalar */
+                    fmt_format<ClusterType>(), /* Format descriptor */
+                    1,                         /* Number of dimensions */
+                    {self.size()},             /* Buffer dimensions */
                     {self.item_size()} /* Strides (in bytes) for each index */
                 );
             });
@@ -232,5 +225,3 @@ void define_3x3_reduction(py::module &m) {
         )doc",
         py::arg("clustervector"));
 }
-
-#pragma GCC diagnostic pop

python/src/bind_Defs.hpp

@@ -1,3 +1,5 @@
+// SPDX-License-Identifier: MPL-2.0
+#include <pybind11/operators.h>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 
@@ -89,6 +91,21 @@ void define_defs_bindings(py::module &m) {
                    self.ymin == other.ymin && self.ymax == other.ymax;
         });
 
+    py::class_<xy>(m, "xy")
+        .def(py::init<uint32_t, uint32_t>(), py::arg("row"), py::arg("col"))
+        .def_readwrite("row", &xy::row)
+        .def_readwrite("col", &xy::col)
+        .def("__repr__",
+             [](const xy &self) {
+                 return fmt::format("xy(row={}, col={})", self.row, self.col);
+             })
+        .def("__iter__",
+             [](const xy &self) {
+                 return py::iter(py::make_tuple(self.row, self.col));
+             })
+        .def(py::self == py::self)
+        .def(py::self != py::self);
+
     py::enum_<DetectorType>(m, "DetectorType")
         .value("Jungfrau", DetectorType::Jungfrau)
         .value("Eiger", DetectorType::Eiger)
@@ -111,4 +128,9 @@ void define_defs_bindings(py::module &m) {
         .value("GATED", TimingMode::GATED)
         .value("BURST_TRIGGER", TimingMode::BURST_TRIGGER)
         .value("TRIGGER_GATED", TimingMode::TRIGGER_GATED);
+
+    py::enum_<FrameDiscardPolicy>(m, "FrameDiscardPolicy")
+        .value("NoDiscard", FrameDiscardPolicy::NoDiscard)
+        .value("Discard", FrameDiscardPolicy::Discard)
+        .value("DiscardPartial", FrameDiscardPolicy::DiscardPartial);
 }

python/src/bind_Eta.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #include "aare/CalculateEta.hpp"
 #include "aare/to_string.hpp"
 

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_MultiThreadedFileReader.hpp

@@ -16,6 +16,7 @@
 #include <pybind11/stl/filesystem.h>
 
 namespace py = pybind11;
+using aare::ssize_t;
 
 inline py::dtype multi_threaded_reader_numpy_dtype(const aare::Dtype &dtype) {
     using aare::Dtype;

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_RawFile.hpp

@@ -318,7 +318,14 @@ void define_raw_file_io_bindings(py::module &m) {
             Cols for the given ROI.
             )")
         .def_property_readonly("bitdepth", &RawFile::bitdepth)
-        .def_property_readonly("geometry", &RawFile::geometry)
+        .def_property_readonly("geometry", &RawFile::geometry, R"(
+            Number of modules along each axis of the detector.
+
+            Returns
+            -------
+            xy
+                Module rows and columns.
+            )")
         .def_property_readonly("detector_type", &RawFile::detector_type)
         .def_property_readonly("master", &RawFile::master)
         .def_property_readonly("n_modules", &RawFile::n_modules)

python/src/bind_calibration.hpp

@@ -9,40 +9,37 @@
 #include <pybind11/stl_bind.h>
 
 namespace py = pybind11;
+using aare::ssize_t;
 
 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 +48,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 +58,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>{};
@@ -73,10 +68,27 @@ py::array_t<T> pybind_calculate_pedestal_g0(
 }
 
 void bind_calibration(py::module &m) {
+    py::options options;
+    options.disable_function_signatures();
+
     m.def("apply_calibration", &pybind_apply_calibration<double>,
           py::arg("raw_data").noconvert(), py::kw_only(),
           py::arg("pd").noconvert(), py::arg("cal").noconvert(),
-          py::arg("n_threads") = 4);
+          py::arg("n_threads") = 4,
+          R"(
+        Apply the calibration to the raw data and return the result as a 3D array of doubles.
+        
+        Parameters
+        ----------
+        raw_data : array_like
+            3D array of shape ``(frames, rows, cols)``.
+        pd : array_like
+            Pedestal array.
+        cal : array_like
+            Calibration array.
+        n_threads : int, optional
+            Number of threads to use. Defaults to 4.
+        )");
 
     m.def("apply_calibration", &pybind_apply_calibration<float>,
           py::arg("raw_data").noconvert(), py::kw_only(),

python/src/bind_testing.hpp

@@ -0,0 +1,65 @@
+// SPDX-License-Identifier: MPL-2.0
+// Probes used by the python test-suite to exercise the helpers in
+// np_helper.hpp. Internal, no API guarantees.
+#pragma once
+
+#include "aare/NDArray.hpp"
+#include "aare/NDView.hpp"
+#include "np_helper.hpp"
+
+#include <algorithm>
+#include <cstdint>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <type_traits>
+
+namespace py = pybind11;
+using namespace ::aare;
+
+// Arrays are taken with the default flags (no c_style) so that nothing is
+// made contiguous by pybind11 before it reaches make_view
+template <typename T, ssize_t Ndim>
+void define_make_view_probes(py::module &m, const std::string &suffix) {
+
+    // Copy what C++ sees through the view back to python
+    m.def(("view_roundtrip_" + suffix).c_str(),
+          [](py::array_t<T> arr) {
+              auto view = make_view<Ndim>(arr);
+              return return_image_data(new NDArray<T, Ndim>(view));
+          },
+          py::arg("arr").noconvert());
+
+    // Same as above but through a read only view
+    m.def(("const_view_roundtrip_" + suffix).c_str(),
+          [](py::array_t<T> arr) {
+              auto view = make_const_view<Ndim>(arr);
+              static_assert(
+                  std::is_same_v<decltype(view), NDView<const T, Ndim>>);
+              auto out = new NDArray<T, Ndim>(view.shape());
+              std::copy(view.begin(), view.end(), out->begin());
+              return return_image_data(out);
+          },
+          py::arg("arr").noconvert());
+
+    // Write through the view
+    m.def(("view_fill_" + suffix).c_str(),
+          [](py::array_t<T> arr, T value) {
+              auto view = make_view<Ndim>(arr);
+              view = value;
+          },
+          py::arg("arr").noconvert(), py::arg("value"));
+}
+
+void define_testing_bindings(py::module &m) {
+    auto testing = m.def_submodule(
+        "testing", "Internal helpers for the test-suite, no API guarantees");
+
+    // noconvert + overloads gives dispatch on dtype
+    define_make_view_probes<double, 1>(testing, "1d");
+    define_make_view_probes<double, 2>(testing, "2d");
+    define_make_view_probes<double, 3>(testing, "3d");
+    define_make_view_probes<uint16_t, 1>(testing, "1d");
+    define_make_view_probes<uint16_t, 2>(testing, "2d");
+    define_make_view_probes<uint16_t, 3>(testing, "3d");
+}

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/file.hpp

@@ -21,11 +21,6 @@
 namespace py = pybind11;
 using namespace ::aare;
 
-// Disable warnings for unused parameters, as we ignore some
-// in the __exit__ method
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wunused-parameter"
-
 void define_file_io_bindings(py::module &m) {
 
     PYBIND11_NUMPY_DTYPE(DetectorHeader, frameNumber, expLength, packetNumber,
@@ -122,9 +117,12 @@ void define_file_io_bindings(py::module &m) {
              })
         .def("__enter__", [](File &self) { return &self; })
         .def("__exit__",
-             [](File &self, const std::optional<pybind11::type> &exc_type,
-                const std::optional<pybind11::object> &exc_value,
-                const std::optional<pybind11::object> &traceback) {
+             []([[maybe_unused]] File &self,
+                [[maybe_unused]] const std::optional<pybind11::type> &exc_type,
+                [[maybe_unused]] const std::optional<pybind11::object>
+                    &exc_value,
+                [[maybe_unused]] const std::optional<pybind11::object>
+                    &traceback) {
                  //  self.close();
              })
         .def("__iter__", [](File &self) { return &self; })
@@ -146,7 +144,7 @@ void define_file_io_bindings(py::module &m) {
                 self.read_into(
                     reinterpret_cast<std::byte *>(image.mutable_data()));
                 return image;
-            } catch (std::runtime_error &e) {
+            } catch (std::runtime_error &) {
                 throw py::stop_iteration();
             }
         });
@@ -160,6 +158,4 @@ void define_file_io_bindings(py::module &m) {
         .def_property_readonly("start", &ScanParameters::start)
         .def_property_readonly("stop", &ScanParameters::stop)
         .def_property_readonly("step", &ScanParameters::step);
-
-#pragma GCC diagnostic pop
 }

python/src/fit.hpp

@@ -11,6 +11,7 @@
 
 namespace py = pybind11;
 using namespace pybind11::literals;
+using aare::ssize_t;
 
 template <typename Model>
 py::object
@@ -22,9 +23,40 @@ fit_dispatch(const aare::FitModel<Model> &model,
 template <typename Model> void bind_fit_model(py::module &m, const char *name) {
     using FM = aare::FitModel<Model>;
     py::class_<FM>(m, name)
-        .def(py::init<unsigned int, unsigned int, double, bool>(),
+        .def(py::init<unsigned int, unsigned int, double, bool,
+                      aare::Minimizer>(),
              py::arg("strategy") = 0, py::arg("max_calls") = 100,
-             py::arg("tolerance") = 0.5, py::arg("compute_errors") = false)
+             py::arg("tolerance") = 0.5, py::arg("compute_errors") = false,
+             py::arg("minimizer") = aare::Minimizer::Migrad,
+             R"doc(
+            Fit configuration for this model.
+
+            Parameters
+            ----------
+            strategy : int
+                Minuit2 strategy (0 = fast, 1 = Minuit2's default). Ignored
+                by ``Minimizer.LevenbergMarquardt`` and
+                ``Minimizer.VarPro``.
+            max_calls : int
+                Work budget per pixel: Minuit2 function calls, or model
+                evaluations for ``Minimizer.LevenbergMarquardt`` and
+                ``Minimizer.VarPro`` (one per iteration, two
+                when a step crosses a limit). 0 selects Minuit's default
+                budget.
+            tolerance : float
+                EDM tolerance in Minuit's convention. Migrad,
+                LevenbergMarquardt and VarPro stop below
+                ``0.002 * tolerance``, Fumili below ``1e-4 * tolerance``.
+            compute_errors : bool
+                Also return parameter errors: from Hesse with Migrad, from
+                the linearised covariance with Fumili, LevenbergMarquardt
+                and VarPro. Fixed parameters and parameters
+                ending on a limit report 0.
+            minimizer : Minimizer
+                ``Minimizer.Migrad`` (default), ``Minimizer.Fumili``,
+                ``Minimizer.LevenbergMarquardt`` or
+                ``Minimizer.VarPro``.
+            )doc")
         .def("SetParLimits",
              py::overload_cast<unsigned int, double, double>(&FM::SetParLimits),
              py::arg("idx"), py::arg("lo"), py::arg("hi"))
@@ -59,8 +91,10 @@ template <typename Model> void bind_fit_model(py::module &m, const char *name) {
         .def_property("tolerance", &FM::tolerance, &FM::SetTolerance)
         .def_property("compute_errors", &FM::compute_errors,
                       &FM::SetComputeErrors)
+        .def_property("minimizer", &FM::minimizer, &FM::SetMinimizer)
         .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 +112,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,
@@ -88,6 +122,11 @@ template <typename Model> void bind_fit_model(py::module &m, const char *name) {
             R"doc(
             Fit this model to 1D or 3D data using Minuit2.
 
+            The minimizer is selected by the ``minimizer`` property
+            (``Minimizer.Migrad`` by default, ``Minimizer.Fumili``,
+            ``Minimizer.LevenbergMarquardt`` or
+            ``Minimizer.VarPro``).
+
             Parameters
             ----------
             x : array_like, shape (n_scan,)
@@ -163,11 +202,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);
@@ -186,8 +220,23 @@ fit_dispatch(const aare::FitModel<Model> &model,
                                 "chi2"_a = return_image_data(chi2_out));
             }
         } else {
-
+            // Unweighted fit; parameter errors are still available when
+            // requested, as in the 1D path.
             NDView<double, 3> dummy_err{};
+
+            if (model.compute_errors()) {
+                auto err_out =
+                    new NDArray<double, 3>({y.shape(0), y.shape(1), npar}, 0.0);
+
+                aare::fit_3d<Model>(model, x_view, y_view, dummy_err,
+                                    par_out->view(), err_out->view(),
+                                    chi2_out->view(), n_threads);
+
+                return py::dict("par"_a = return_image_data(par_out),
+                                "par_err"_a = return_image_data(err_out),
+                                "chi2"_a = return_image_data(chi2_out));
+            }
+
             NDView<double, 3> dummy_err_out{};
 
             aare::fit_3d<Model>(model, x_view, y_view, dummy_err,
@@ -208,11 +257,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 {
@@ -226,7 +270,82 @@ fit_dispatch(const aare::FitModel<Model> &model,
     }
 }
 
+// Resolve the model type behind a Python object and run the fit.
+py::object dispatch_any_model(
+    py::object model_obj,
+    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) {
+    using namespace aare::model;
+
+#define AARE_DISPATCH_MODEL(Model)                                             \
+    if (py::isinstance<aare::FitModel<Model>>(model_obj)) {                    \
+        const auto &mdl = model_obj.cast<const aare::FitModel<Model> &>();     \
+        return fit_dispatch<Model>(mdl, x, y, y_err_obj, n_threads);           \
+    }
+
+    AARE_DISPATCH_MODEL(Pol1)
+    AARE_DISPATCH_MODEL(Pol2)
+    AARE_DISPATCH_MODEL(Gaussian)
+    AARE_DISPATCH_MODEL(GaussianErfcPlateau)
+    AARE_DISPATCH_MODEL(GaussianChargeSharing)
+    AARE_DISPATCH_MODEL(GaussianChargeSharingKb)
+    AARE_DISPATCH_MODEL(RisingScurve)
+    AARE_DISPATCH_MODEL(FallingScurve)
+
+#undef AARE_DISPATCH_MODEL
+
+    throw std::runtime_error(
+        "Unknown model type. Expected Pol1, Pol2, Gaussian, "
+        "GaussianErfcPlateau, GaussianChargeSharing, GaussianChargeSharingKb, "
+        "RisingScurve or FallingScurve.");
+}
+
 void define_fit_bindings(py::module &m) {
+    // The enum must exist before it is used as a constructor default below.
+    py::enum_<aare::Minimizer>(m, "Minimizer", R"doc(
+        Minimizer used by the fit models.
+
+        ``Migrad`` (default) is Minuit2's variable-metric minimizer. With
+        ``compute_errors`` its parameter errors come from Hesse.
+
+        ``Fumili`` is Minuit2's Gauss-Newton minimizer with a trust region.
+        It takes the gradient and a linearised Hessian from the analytic
+        derivatives of the model and usually needs far fewer function
+        evaluations than Migrad. Its parameter errors come from the
+        covariance of the linearised Hessian, without an extra Hesse step.
+        Minuit2's Fumili cannot converge once a two-sided limit becomes
+        active; a pixel without a valid Fumili minimum is refitted with
+        Migrad.
+
+        ``LevenbergMarquardt`` is the built-in damped Gauss-Newton solver. It
+        uses the analytic Jacobian, reflects steps at parameter limits and
+        reuses its buffers between pixels, so data cubes are fitted without
+        per-pixel allocations. Its parameter errors are
+        ``sqrt(diag((J^T J)^-1))``; ``max_calls`` counts model evaluations
+        and ``strategy`` is ignored.
+
+        ``VarPro`` is the built-in variable projection solver.
+        Every bundled model is linear in some of its parameters (the
+        amplitude of a Gaussian, the baseline and amplitude of an S-curve,
+        every parameter of a polynomial); each trial point solves those
+        exactly and the Levenberg-Marquardt iteration runs over the
+        remaining nonlinear parameters only. It needs fewer evaluations than
+        ``LevenbergMarquardt`` and cannot be misled by poor start values of
+        the linear parameters, which it ignores. A pixel whose linear
+        solution violates a limit on a linear parameter, or that does not
+        converge, is refitted with ``LevenbergMarquardt``. Its parameter
+        errors are the same Gauss-Newton estimates; ``max_calls`` counts
+        evaluations and ``strategy`` is ignored.
+
+        With every minimizer, fixed parameters and parameters that end on a
+        limit report an error of 0 and a failed fit returns zeros.
+        )doc")
+        .value("Migrad", aare::Minimizer::Migrad)
+        .value("Fumili", aare::Minimizer::Fumili)
+        .value("LevenbergMarquardt", aare::Minimizer::LevenbergMarquardt)
+        .value("VarPro", aare::Minimizer::VarPro);
+
     // ── Bind model classes ──────────────────────────────────────────
     bind_fit_model<aare::model::Gaussian>(m, "Gaussian");
     bind_fit_model<aare::model::GaussianErfcPlateau>(m, "GaussianErfcPlateau");
@@ -241,90 +360,24 @@ 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 {
-            using namespace aare::model;
-
-            // ── Polynomial of degree 1 ───────
-            if (py::isinstance<aare::FitModel<Pol1>>(model_obj)) {
-                const auto &mdl =
-                    model_obj.cast<const aare::FitModel<Pol1> &>();
-                return fit_dispatch<Pol1>(mdl, x, y, y_err_obj, n_threads);
-            }
-
-            // ── Polynomial of degree 2 ───────
-            if (py::isinstance<aare::FitModel<Pol2>>(model_obj)) {
-                const auto &mdl =
-                    model_obj.cast<const aare::FitModel<Pol2> &>();
-                return fit_dispatch<Pol2>(mdl, x, y, y_err_obj, n_threads);
-            }
-
-            // ── Gaussian ───────
-            if (py::isinstance<aare::FitModel<Gaussian>>(model_obj)) {
-                const auto &mdl =
-                    model_obj.cast<const aare::FitModel<Gaussian> &>();
-                return fit_dispatch<Gaussian>(mdl, x, y, y_err_obj, n_threads);
-            }
-
-            // ── GaussianErfcPlateau ───────
-            if (py::isinstance<aare::FitModel<GaussianErfcPlateau>>(
-                    model_obj)) {
-                const auto &mdl =
-                    model_obj
-                        .cast<const aare::FitModel<GaussianErfcPlateau> &>();
-                return fit_dispatch<GaussianErfcPlateau>(mdl, x, y, y_err_obj,
-                                                         n_threads);
-            }
-
-            // ── GaussianChargeSharing ───────
-            if (py::isinstance<aare::FitModel<GaussianChargeSharing>>(
-                    model_obj)) {
-                const auto &mdl =
-                    model_obj
-                        .cast<const aare::FitModel<GaussianChargeSharing> &>();
-                return fit_dispatch<GaussianChargeSharing>(mdl, x, y, y_err_obj,
-                                                           n_threads);
-            }
-
-            // ── GaussianChargeSharingKb ───────
-            if (py::isinstance<aare::FitModel<GaussianChargeSharingKb>>(
-                    model_obj)) {
-                const auto &mdl = model_obj.cast<
-                    const aare::FitModel<GaussianChargeSharingKb> &>();
-                return fit_dispatch<GaussianChargeSharingKb>(
-                    mdl, x, y, y_err_obj, n_threads);
-            }
-
-            // ── Rising Scurve ───────
-            if (py::isinstance<aare::FitModel<RisingScurve>>(model_obj)) {
-                const auto &mdl =
-                    model_obj.cast<const aare::FitModel<RisingScurve> &>();
-                return fit_dispatch<RisingScurve>(mdl, x, y, y_err_obj,
-                                                  n_threads);
-            }
-
-            // ── Falling Scurve ───────
-            if (py::isinstance<aare::FitModel<FallingScurve>>(model_obj)) {
-                const auto &mdl =
-                    model_obj.cast<const aare::FitModel<FallingScurve> &>();
-                return fit_dispatch<FallingScurve>(mdl, x, y, y_err_obj,
-                                                   n_threads);
-            }
-
-            throw std::runtime_error(
-                "Unknown model type. Expected Pol1, Pol2, Gaussian, "
-                "RisingScurve or FallingScurve.");
+            return dispatch_any_model(model_obj, x, y, y_err_obj, n_threads);
         },
         R"(
-        Fit a model to 1D or 3D data using Minuit2.
+        Fit a model to 1D or 3D data with the minimizer selected by the
+        model's ``minimizer`` property.
  
         Parameters
         ----------
-        model : Pol1, Pol2, Gaussian, RisingScurve, or FallingScurve
-            Configured model object.  User-set limits, fixed parameters,
-            and start values take precedence over automatic estimates.
+        model : object
+            Configured model object: Pol1, Pol2, Gaussian, GaussianErfcPlateau,
+            GaussianChargeSharing, GaussianChargeSharingKb, RisingScurve or
+            FallingScurve.  User-set limits, fixed parameters, and start values
+            take precedence over automatic estimates.
         x : array_like, shape (n_scan,)
             Scan points (e.g. energy or threshold values).
         y : array_like, shape (n_scan,) or (rows, cols, n_scan)
@@ -336,15 +389,12 @@ void define_fit_bindings(py::module &m) {
  
         Returns
         -------
-        For 1D input:
-            numpy array of shape (2*npar+1,) if compute_errors else (npar+1,).
-            Layout: [params..., (errors...,) chi2].
- 
-        For 3D input:
-            dict with keys:
-              "par"     : (rows, cols, npar) fitted parameters.
-              "par_err" : (rows, cols, npar) parameter errors (if compute_errors).
-              "chi2"    : (rows, cols)       chi-squared per pixel.
+        dict
+            ``par`` holds the fitted parameters, shape (npar,) for 1D input or
+            (rows, cols, npar) for 3D input.  ``chi2`` holds the chi-squared,
+            shape (1,) or (rows, cols).  ``par_err`` holds the parameter
+            errors when ``compute_errors`` is enabled; fixed parameters and
+            parameters on a limit report 0.  A failed fit returns zeros.
         )",
         py::arg("model"), py::arg("x"), py::arg("y"),
         py::arg("y_err") = py::none(), py::arg("n_threads") = 4);

python/src/jungfrau_data_file.hpp

@@ -15,11 +15,6 @@
 namespace py = pybind11;
 using namespace ::aare;
 
-// Disable warnings for unused parameters, as we ignore some
-// in the __exit__ method
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wunused-parameter"
-
 auto read_dat_frame(JungfrauDataFile &self) {
     py::array_t<JungfrauDataHeader> header(1);
     py::array_t<uint16_t> image({self.rows(), self.cols()});
@@ -94,20 +89,20 @@ void define_jungfrau_data_file_io_bindings(py::module &m) {
               )")
         .def("__enter__", [](JungfrauDataFile &self) { return &self; })
         .def("__exit__",
-             [](JungfrauDataFile &self,
-                const std::optional<pybind11::type> &exc_type,
-                const std::optional<pybind11::object> &exc_value,
-                const std::optional<pybind11::object> &traceback) {
+             []([[maybe_unused]] JungfrauDataFile &self,
+                [[maybe_unused]] const std::optional<pybind11::type> &exc_type,
+                [[maybe_unused]] const std::optional<pybind11::object>
+                    &exc_value,
+                [[maybe_unused]] const std::optional<pybind11::object>
+                    &traceback) {
                  //  self.close();
              })
         .def("__iter__", [](JungfrauDataFile &self) { return &self; })
         .def("__next__", [](JungfrauDataFile &self) {
             try {
                 return read_dat_frame(self);
-            } catch (std::runtime_error &e) {
+            } catch (std::runtime_error &) {
                 throw py::stop_iteration();
             }
         });
 }
-
-#pragma GCC diagnostic pop

python/src/module.cpp

@@ -22,6 +22,7 @@
 #include "bind_PixelMap.hpp"
 #include "bind_RawFile.hpp"
 #include "bind_calibration.hpp"
+#include "bind_testing.hpp"
 
 #include "StrixelRemap/bind_InclusiveROI.hpp"
 #include "StrixelRemap/bind_PredefinedVariables.hpp"
@@ -98,6 +99,7 @@ PYBIND11_MODULE(_aare, m) {
     define_jungfrau_data_file_io_bindings(m);
 
     bind_calibration(m);
+    define_testing_bindings(m);
 
     DEFINE_CLUSTER_BINDINGS(int, 3, 3, uint16_t, i);
     DEFINE_CLUSTER_BINDINGS(double, 3, 3, uint16_t, d);

python/src/module_config.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #pragma once
 #include <cstdint>
 

python/src/np_helper.hpp

@@ -41,31 +41,60 @@ template <typename T> py::array return_vector(std::vector<T> *vec) {
                           free_when_done); // numpy array references this parent
 }
 
-// todo rewrite generic
-template <class T, int Flags>
-auto get_shape_3d(const py::array_t<T, Flags> &arr) {
-    return aare::Shape<3>{arr.shape(0), arr.shape(1), arr.shape(2)};
+// Create a NDView of a numpy array. NDView assumes C-order strides so we
+// reject anything that is not Ndim dimensional and C-contiguous instead of
+// silently reading the data in the wrong order.
+template <ssize_t Ndim, class T, int Flags>
+auto get_checked_shape(const py::array_t<T, Flags> &arr) {
+    if (arr.ndim() != Ndim) {
+        throw py::value_error(
+            fmt::format("Expected {}D array, got {}D", Ndim, arr.ndim()));
+    }
+    if (!(arr.flags() & py::array::c_style)) {
+        throw py::value_error(
+            "Array is not C-contiguous. Use np.ascontiguousarray(arr)");
+    }
+    aare::Shape<Ndim> shape{};
+    for (ssize_t i = 0; i < Ndim; ++i) {
+        shape[i] = arr.shape(i);
+    }
+    return shape;
+}
+
+template <ssize_t Ndim, class T, int Flags>
+auto make_view(py::array_t<T, Flags> &arr) {
+    auto shape = get_checked_shape<Ndim>(arr);
+    return aare::NDView<T, Ndim>(arr.mutable_data(), shape);
+}
+
+// Read only view, also works for numpy arrays that are not writeable
+template <ssize_t Ndim, class T, int Flags>
+auto make_const_view(const py::array_t<T, Flags> &arr) {
+    auto shape = get_checked_shape<Ndim>(arr);
+    return aare::NDView<const T, Ndim>(arr.data(), shape);
 }
 
 template <class T, int Flags> auto make_view_3d(py::array_t<T, Flags> &arr) {
-    return aare::NDView<T, 3>(arr.mutable_data(), get_shape_3d<T, Flags>(arr));
+    return make_view<3>(arr);
 }
-
-template <class T, int Flags>
-auto get_shape_2d(const py::array_t<T, Flags> &arr) {
-    return aare::Shape<2>{arr.shape(0), arr.shape(1)};
-}
-
-template <class T, int Flags>
-auto get_shape_1d(const py::array_t<T, Flags> &arr) {
-    return aare::Shape<1>{arr.shape(0)};
-}
-
 template <class T, int Flags> auto make_view_2d(py::array_t<T, Flags> &arr) {
-    return aare::NDView<T, 2>(arr.mutable_data(), get_shape_2d<T, Flags>(arr));
+    return make_view<2>(arr);
 }
 template <class T, int Flags> auto make_view_1d(py::array_t<T, Flags> &arr) {
-    return aare::NDView<T, 1>(arr.mutable_data(), get_shape_1d<T, Flags>(arr));
+    return make_view<1>(arr);
+}
+
+template <class T, int Flags>
+auto make_const_view_3d(const py::array_t<T, Flags> &arr) {
+    return make_const_view<3>(arr);
+}
+template <class T, int Flags>
+auto make_const_view_2d(const py::array_t<T, Flags> &arr) {
+    return make_const_view<2>(arr);
+}
+template <class T, int Flags>
+auto make_const_view_1d(const py::array_t<T, Flags> &arr) {
+    return make_const_view<1>(arr);
 }
 
 template <typename ClusterType> struct fmt_format_trait; // forward declaration
@@ -83,8 +112,13 @@ struct fmt_format_trait<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
     }
 };
 
-template <typename ClusterType>
-auto fmt_format = fmt_format_trait<ClusterType>::value();
+// Function-local static rather than a variable template: MSVC does not run
+// the dynamic initializer of a variable template instantiated from inside
+// another template, which left the format string empty.
+template <typename ClusterType> const std::string &fmt_format() {
+    static const std::string value = fmt_format_trait<ClusterType>::value();
+    return value;
+}
 
 /**
  * Helper function to allocate image data given item size and shape

python/src/raw_master_file.hpp

@@ -67,9 +67,23 @@ void define_raw_master_file_bindings(py::module &m) {
         .def_property_readonly("total_frames_expected",
                                &RawMasterFile::total_frames_expected)
         .def_property_readonly("detector_layout",
-                               &RawMasterFile::detector_layout)
+                               &RawMasterFile::detector_layout, R"(
+            Number of modules along each axis of the detector
+
+            Returns
+            ----------
+                xy
+                    Module rows and columns
+            )")
         .def_property_readonly("udp_interfaces_per_module",
-                               &RawMasterFile::udp_interfaces_per_module)
+                               &RawMasterFile::udp_interfaces_per_module, R"(
+            Number of UDP interfaces along each axis of a module
+
+            Returns
+            ----------
+                xy
+                    UDP interface rows and columns
+            )")
         .def_property_readonly("analog_samples", &RawMasterFile::analog_samples,
                                R"(
             Number of analog samples

python/src/raw_sub_file.hpp

@@ -61,11 +61,6 @@ auto read_n_frames_from_RawSubFile(RawSubFile &self, size_t n_frames) {
     return py::make_tuple(header, image);
 }
 
-// Disable warnings for unused parameters, as we ignore some
-// in the __exit__ method
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wunused-parameter"
-
 void define_raw_sub_file_io_bindings(py::module &m) {
     py::class_<RawSubFile>(m, "RawSubFile")
         .def(py::init<const std::filesystem::path &, DetectorType, size_t,
@@ -89,17 +84,18 @@ void define_raw_sub_file_io_bindings(py::module &m) {
              })
         .def("__enter__", [](RawSubFile &self) { return &self; })
         .def("__exit__",
-             [](RawSubFile &self, const std::optional<pybind11::type> &exc_type,
-                const std::optional<pybind11::object> &exc_value,
-                const std::optional<pybind11::object> &traceback) {})
+             []([[maybe_unused]] RawSubFile &self,
+                [[maybe_unused]] const std::optional<pybind11::type> &exc_type,
+                [[maybe_unused]] const std::optional<pybind11::object>
+                    &exc_value,
+                [[maybe_unused]] const std::optional<pybind11::object>
+                    &traceback) {})
         .def("__iter__", [](RawSubFile &self) { return &self; })
         .def("__next__", [](RawSubFile &self) {
             try {
                 return read_frame_from_RawSubFile(self);
-            } catch (std::runtime_error &e) {
+            } catch (std::runtime_error &) {
                 throw py::stop_iteration();
             }
         });
 }
-
-#pragma GCC diagnostic pop

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

@@ -135,6 +135,43 @@ def test_cluster_finder():
     assert clusters.size == 0
 
 
+def test_cluster_finder_set_nSigma_before_pedestal_is_ready():
+    shape = [10, 10]
+    n_samples = 4
+    cf = _aare.ClusterFinder_Cluster3x3i(shape, min_pedestal_samples=n_samples)
+    assert not cf.pedestal_ready
+
+    cf.nSigma = 50
+    assert cf.nSigma == 50
+
+    # Pedestal of 100 with a noise of 2
+    for i in range(n_samples):
+        cf.push_pedestal_frame(np.full(shape, 98 + 4 * (i % 2), dtype=np.uint16))
+    assert cf.pedestal_ready
+
+    # 30 ADU is a cluster at the default nSigma of 5 but not at 50
+    frame = np.full(shape, 100, dtype=np.uint16)
+    frame[5, 5] = 130
+    cf.find_clusters(frame)
+    assert cf.steal_clusters().size == 0
+
+
+def test_cluster_finder_mt_set_nSigma_before_pedestal_is_ready():
+    shape = [10, 10]
+    n_samples = 4
+    cf = _aare.ClusterFinderMT_Cluster3x3i(
+        shape, n_threads=2, min_pedestal_samples=n_samples
+    )
+    assert not cf.pedestal_ready()
+
+    cf.set_nSigma(50)
+
+    for i in range(n_samples):
+        cf.push_pedestal_frame(np.full(shape, 98 + 4 * (i % 2), dtype=np.uint16))
+    assert cf.pedestal_ready()
+    cf.stop()
+
+
 def test_2x2_reduction(): 
     """Test 2x2 Reduction"""
     cluster = _aare.Cluster3x3i(5,5,np.array([1, 1, 1, 2, 3, 1, 2, 2, 1], dtype=np.int32))
@@ -162,3 +199,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)
 
 
@@ -157,5 +157,5 @@ def test_masking_requires_one_dimension():
 
     mask = np.array([[True, False]], dtype=bool)
 
-    with pytest.raises(ValueError, match="one-dimensional"):
+    with pytest.raises(ValueError):
         cv(mask)

python/tests/test_CtbRawFile.py

@@ -0,0 +1,59 @@
+# SPDX-License-Identifier: MPL-2.0
+
+
+import pytest
+import numpy as np
+from aare import CtbRawFile
+
+
+@pytest.mark.withdata
+def test_basic_properites(test_data_path):
+    f = CtbRawFile(test_data_path/'raw/ctb/run_master_0.json')
+    assert f.total_frames == 47
+    assert f.image_size_in_bytes == 640
+    assert f.master.analog_samples == 10
+    assert f.master.digital_samples is None
+    assert f.master.transceiver_samples is None
+
+
+@pytest.mark.withdata
+def test_seek_and_read(test_data_path):
+    f = CtbRawFile(test_data_path/'raw/ctb/run_master_0.json')
+    f.seek(0)
+    header, raw_data = f.read_frame()
+    assert header['frameNumber'] == 48 #we know that the first frame in this file is 48
+
+    #Seek a multiple of frames in file
+    f.seek(9)
+    header, raw_data = f.read_frame()
+    assert header['frameNumber'] == 57 
+
+@pytest.mark.withdata
+def test_read_all_frames_direct_and_iterate(test_data_path):
+    """
+    We want to make sure that iterating and reading give the same result.
+    """
+    f0 = CtbRawFile(test_data_path/'raw/ctb/run_master_0.json')
+    f1 = CtbRawFile(test_data_path/'raw/ctb/run_master_0.json')
+
+    for h0, d0 in f0:
+        h1, d1 = f1.read_frame()
+        assert h0 == h1
+        assert np.all(d0 == d1)
+
+@pytest.mark.withdata
+def test_read_all_frames_with_index_and_iterate(test_data_path):
+    """
+    We want to make sure that iterating and reading give the same result.
+    """
+    f0 = CtbRawFile(test_data_path/'raw/ctb/run_master_0.json')
+    f1 = CtbRawFile(test_data_path/'raw/ctb/run_master_0.json')
+
+    frames_read = 0
+    for i,(h0, d0) in enumerate(f0):
+        h1, d1 = f1.read_frame(i)
+        assert h0 == h1
+        assert np.all(d0 == d1)
+        frames_read += 1
+
+    assert frames_read == f0.total_frames

python/tests/test_FastPedestal.py

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 import numpy as np
 import pytest
 

python/tests/test_Fit.py

@@ -1,15 +1,401 @@
 # SPDX-License-Identifier: MPL-2.0
 import numpy as np
+import pytest
+
 import aare
 
+MINIMIZERS = [
+    aare.Minimizer.Migrad,
+    aare.Minimizer.Fumili,
+    aare.Minimizer.LevenbergMarquardt,
+    aare.Minimizer.VarPro,
+]
 
 
-def test_gaussian_model_evaluates_and_fits_data():
+def gaussian(x, A, mu, sigma):
+    return A * np.exp(-0.5 * ((x - mu) / sigma) ** 2)
+
+
+def uses_minuit2(minimizer):
+    return minimizer not in (
+        aare.Minimizer.LevenbergMarquardt,
+        aare.Minimizer.VarPro,
+    )
+
+
+@pytest.fixture(params=MINIMIZERS, ids=lambda m: m.name)
+def minimizer(request):
+    return request.param
+
+
+@pytest.fixture
+def gaussian_data():
+    x = np.linspace(-6.0, 6.0, 61)
+    y = gaussian(x, 120.0, 0.8, 1.3)
+    y_err = np.ones_like(y)
+    return x, y, y_err
+
+
+def test_gaussian_model_evaluates_and_fits_data(minimizer):
     x = np.linspace(-5.0, 5.0, 51)
     expected = np.array([20.0, 0.5, 1.2])
-    model = aare.Gaussian()
+    model = aare.Gaussian(minimizer=minimizer)
 
     y = model(x, expected)
     result = model.fit(x, y)
 
     np.testing.assert_allclose(result["par"], expected, atol=2e-3)
+
+
+def test_minimizer_can_be_selected():
+    assert aare.Gaussian().minimizer == aare.Minimizer.Migrad
+    assert set(aare.Minimizer.__members__) == {
+        "Migrad",
+        "Fumili",
+        "LevenbergMarquardt",
+        "VarPro",
+    }
+
+    for minimizer in MINIMIZERS:
+        assert aare.Gaussian(minimizer=minimizer).minimizer == minimizer
+
+    model = aare.Gaussian()
+    model.minimizer = aare.Minimizer.LevenbergMarquardt
+    assert model.minimizer == aare.Minimizer.LevenbergMarquardt
+
+
+def test_fit_returns_errors_and_chi2(gaussian_data, minimizer):
+    x, y, y_err = gaussian_data
+    model = aare.Gaussian(compute_errors=True, minimizer=minimizer)
+    result = model.fit(x, y, y_err)
+
+    assert set(result) == {"par", "par_err", "chi2"}
+    assert result["par"].shape == (3,)
+    assert result["par_err"].shape == (3,)
+    assert result["chi2"].shape == (1,)
+    np.testing.assert_allclose(result["par"], [120.0, 0.8, 1.3], rtol=1e-4)
+    np.testing.assert_allclose(
+        result["par_err"], [0.361, 0.00451, 0.00451], rtol=2e-2
+    )
+    assert result["chi2"][0] == pytest.approx(0.0, abs=1e-3)
+
+    without_errors = aare.fit(aare.Gaussian(minimizer=minimizer), x, y)
+    assert set(without_errors) == {"par", "chi2"}
+
+
+def test_pol1_weighted_errors_match_analytic_formula(minimizer):
+    x = np.arange(20, dtype=float)
+    y = 2.0 + 0.5 * x
+    y_err = 1.0 + 0.1 * x
+    w = 1.0 / y_err**2
+    S, Sx, Sxx = w.sum(), (w * x).sum(), (w * x * x).sum()
+    delta = S * Sxx - Sx * Sx
+
+    result = aare.Pol1(compute_errors=True, minimizer=minimizer).fit(x, y, y_err)
+
+    np.testing.assert_allclose(result["par"], [2.0, 0.5], rtol=1e-6)
+    np.testing.assert_allclose(
+        result["par_err"],
+        [np.sqrt(Sxx / delta), np.sqrt(S / delta)],
+        rtol=1e-6,
+    )
+
+
+def test_fixed_parameter_reports_zero_error(gaussian_data, minimizer):
+    x, y, y_err = gaussian_data
+    model = aare.Gaussian(compute_errors=True, minimizer=minimizer)
+    model.FixParameter("sigma", 1.5)
+    result = model.fit(x, y, y_err)
+
+    assert result["par"][2] == 1.5
+    assert result["par_err"][2] == 0.0
+    assert result["par"][0] == pytest.approx(111.15, rel=1e-3)
+    assert result["chi2"][0] == pytest.approx(1684.25, rel=1e-3)
+
+    model.ReleaseParameter("sigma")
+    result = model.fit(x, y, y_err)
+    assert result["par"][2] == pytest.approx(1.3, rel=1e-4)
+
+
+def test_active_limit_ends_on_the_limit(gaussian_data, minimizer):
+    x, y, y_err = gaussian_data
+    model = aare.Gaussian(compute_errors=True, minimizer=minimizer)
+    model.SetParLimits("mu", -5.0, 0.5)
+    result = model.fit(x, y, y_err)
+
+    margin = 1e-6 if uses_minuit2(minimizer) else 1e-9
+    assert result["par"][1] == pytest.approx(0.5, abs=margin)
+    assert result["par_err"][1] == 0.0
+    assert result["par"][0] == pytest.approx(116.88, rel=1e-3)
+    assert result["chi2"][0] == pytest.approx(4301.95, rel=1e-3)
+
+    inactive = aare.Gaussian(minimizer=minimizer)
+    inactive.SetParLimits(1, -5.0, 5.0)
+    np.testing.assert_allclose(
+        inactive.fit(x, y)["par"], [120.0, 0.8, 1.3], rtol=1e-4
+    )
+
+
+def test_user_start_value(gaussian_data, minimizer):
+    x, y, _ = gaussian_data
+    model = aare.Gaussian(minimizer=minimizer)
+    model.SetParameter("mu", 0.0)
+    np.testing.assert_allclose(model.fit(x, y)["par"], [120.0, 0.8, 1.3], rtol=1e-4)
+
+
+def test_max_calls_failure_returns_zeros(gaussian_data, minimizer):
+    x, y, y_err = gaussian_data
+    model = aare.Gaussian(max_calls=3, compute_errors=True, minimizer=minimizer)
+    # start far from the minimum so that no minimizer converges in 3 calls
+    model.SetParameter("mu", -3.0)
+    model.SetParameter("sigma", 3.0)
+    result = model.fit(x, y, y_err)
+    assert np.all(result["par"] == 0.0)
+    assert np.all(result["par_err"] == 0.0)
+    assert np.all(result["chi2"] == 0.0)
+
+
+def test_strategy_is_ignored_by_levenberg_marquardt(gaussian_data):
+    x, y, _ = gaussian_data
+    lm = aare.Minimizer.LevenbergMarquardt
+    a = aare.Gaussian(strategy=0, minimizer=lm).fit(x, y)
+    b = aare.Gaussian(strategy=1, minimizer=lm).fit(x, y)
+    np.testing.assert_array_equal(a["par"], b["par"])
+
+
+def test_invalid_parameter_access_raises():
+    model = aare.Gaussian()
+    assert model.par_names == ["A", "mu", "sigma"]
+    assert model.n_par == 3
+    with pytest.raises(RuntimeError):
+        model.SetParameter("nope", 1.0)
+    with pytest.raises(IndexError):
+        model.SetParameter(3, 1.0)
+    with pytest.raises(RuntimeError):
+        model.SetParLimits("mu", 1.0, 1.0)
+
+
+def test_mismatched_array_sizes_raise(gaussian_data, minimizer):
+    x, y, y_err = gaussian_data
+    model = aare.Gaussian(minimizer=minimizer)
+    with pytest.raises(RuntimeError):
+        model.fit(x, y[:-1])
+    with pytest.raises(RuntimeError):
+        model.fit(x, y, y_err[:-1])
+
+
+@pytest.mark.parametrize("cls", [aare.RisingScurve, aare.FallingScurve])
+def test_scurve_fit(cls, minimizer):
+    truth = np.array([10.0, 0.1, 50.0, 3.0, 1000.0, 2.0])
+    x = np.linspace(0.0, 100.0, 100)
+    model = cls(max_calls=500, compute_errors=True, minimizer=minimizer)
+    y = model(x, truth)
+    result = model.fit(x, y, np.ones_like(y))
+
+    np.testing.assert_allclose(result["par"], truth, rtol=1e-3)
+    assert np.all(result["par_err"] > 0.0)
+
+
+@pytest.mark.parametrize("with_errors", [False, True])
+@pytest.mark.parametrize("weighted", [False, True])
+def test_fit_3d_matches_fit_1d(with_errors, weighted, minimizer):
+    x = np.linspace(-6.0, 6.0, 61)
+    rows, cols = 4, 3
+    y = np.empty((rows, cols, x.size))
+    for r in range(rows):
+        for c in range(cols):
+            y[r, c] = gaussian(x, 80.0 + 10.0 * r, -0.6 + 0.3 * c, 0.9 + 0.1 * r)
+    y_err = np.ones_like(y) if weighted else None
+
+    model = aare.Gaussian(compute_errors=with_errors, minimizer=minimizer)
+    cube = model.fit(x, y, y_err, n_threads=3)
+
+    expected_keys = {"par", "par_err", "chi2"} if with_errors else {"par", "chi2"}
+    assert set(cube) == expected_keys
+    assert cube["par"].shape == (rows, cols, 3)
+    assert cube["chi2"].shape == (rows, cols)
+    if with_errors:
+        assert cube["par_err"].shape == (rows, cols, 3)
+    for r in range(rows):
+        for c in range(cols):
+            single = model.fit(x, y[r, c], None if y_err is None else y_err[r, c])
+            np.testing.assert_allclose(cube["par"][r, c], single["par"], rtol=1e-10)
+            assert cube["chi2"][r, c] == pytest.approx(single["chi2"][0], abs=1e-10)
+            if with_errors:
+                np.testing.assert_allclose(
+                    cube["par_err"][r, c], single["par_err"], rtol=1e-10
+                )
+    np.testing.assert_allclose(cube["par"][3, 2], [110.0, 0.0, 1.2], rtol=1e-4, atol=1e-6)
+
+
+def test_all_minimizers_agree_on_3d_data_with_errors():
+    x = np.linspace(-5.0, 5.0, 51)
+    expected = np.array([20.0, 0.5, 1.2])
+    single = aare.Gaussian()(x, expected)
+    y = np.tile(single, (2, 3, 1))
+    y_err = np.ones_like(y)
+
+    results = {
+        minimizer: aare.Gaussian(compute_errors=True, minimizer=minimizer).fit(
+            x, y, y_err, n_threads=2
+        )
+        for minimizer in MINIMIZERS
+    }
+    reference = results[aare.Minimizer.Migrad]
+    for minimizer, result in results.items():
+        assert result["par"].shape == (2, 3, 3)
+        assert result["par_err"].shape == (2, 3, 3)
+        np.testing.assert_allclose(result["par"], np.broadcast_to(expected, (2, 3, 3)), atol=2e-3)
+        np.testing.assert_allclose(result["chi2"], 0.0, atol=1e-3)
+        np.testing.assert_allclose(result["par_err"], reference["par_err"], rtol=2e-2)
+
+
+# Regressions of the built-in solvers: fits that used to be reported as valid
+# although they had stopped away from the minimum or returned garbage.
+
+
+def test_gaussian_fit_to_noise_fails_or_stays_bounded(minimizer):
+    # A width collapsing between two scan points let the amplitude diverge
+    # (up to 1e159) in a fit that was still reported as valid.
+    x = np.linspace(-6.0, 6.0, 61)
+    y = np.random.default_rng(0).normal(0.0, 1.0, (20, 20, x.size))
+    result = aare.Gaussian(minimizer=minimizer).fit(x, y, n_threads=2)
+
+    par = result["par"]
+    failed = np.all(par == 0.0, axis=-1)
+    bounded = np.all(np.isfinite(par) & (np.abs(par) < 1e6), axis=-1)
+    assert np.all(failed | bounded)
+    assert np.all(np.isfinite(result["chi2"]))
+
+
+def test_step_reflected_onto_start_is_not_convergence(minimizer):
+    # The step of +2 reflects at the limit 1 back onto the start 0; the
+    # truncated step must still be tried.
+    model = aare.Pol1(minimizer=minimizer)
+    model.SetParameter("p0", 0.0)
+    model.FixParameter("p1", 0.0)
+    model.SetParLimits("p0", -1.0, 1.0)
+    result = model.fit(np.arange(4.0), np.full(4, 2.2))
+
+    assert result["par"][0] == pytest.approx(1.0, abs=1e-6)
+    assert result["chi2"][0] == pytest.approx(5.76, rel=1e-6)
+
+
+@pytest.mark.parametrize(
+    "amplitude, offset", [(1e8, 0.0), (1e10, 0.0), (120.0, 1e9)]
+)
+def test_large_parameter_does_not_stop_the_fit(minimizer, amplitude, offset):
+    # One relative step norm over all parameters let a large amplitude or
+    # position stop the fit at the start estimate of the width.
+    x = np.linspace(-6.0, 6.0, 61)
+    y = gaussian(x, amplitude, 0.8, 1.3)
+    y_err = np.sqrt(y) + 1.0
+    result = aare.Gaussian(minimizer=minimizer, max_calls=500).fit(
+        x + offset, y, y_err
+    )
+
+    par = result["par"]
+    assert par[0] == pytest.approx(amplitude, rel=1e-3)
+    assert par[1] - offset == pytest.approx(0.8, abs=1e-3)
+    assert par[2] == pytest.approx(1.3, rel=1e-3)
+
+
+def test_rejected_trial_on_limit_does_not_freeze_interior_parameter(minimizer):
+    # A rejected trial truncated onto the lower limit froze mu at its start.
+    x = np.linspace(-0.1, 0.1, 11)
+    mu = np.sqrt(2.0 * np.log(2.0))
+    y = gaussian(x, 10.0, mu, 1.0)
+    model = aare.Gaussian(minimizer=minimizer)
+    model.FixParameter("A", 10.0)
+    model.FixParameter("sigma", 1.0)
+    model.SetParameter("mu", 3.0)
+    model.SetParLimits("mu", 0.0, 5.0)
+    result = model.fit(x, y)
+
+    assert result["par"][1] == pytest.approx(mu, abs=1e-3)
+    assert result["chi2"][0] < 1e-4
+
+
+# Samples with a zero error are ignored, also when they hold NaN (a masked
+# dead frame): a zero weight does not remove a NaN from a sum.
+
+
+@pytest.mark.parametrize("masked", [4, 0, 9])
+def test_masked_nan_sample_is_ignored(minimizer, masked):
+    # The start estimate of Pol1 uses the first and the last sample.
+    x = np.arange(10.0)
+    y = 2.0 + 0.5 * x
+    y_err = np.ones_like(y)
+    y[masked] = np.nan
+    y_err[masked] = 0.0
+    result = aare.Pol1(minimizer=minimizer, compute_errors=True).fit(x, y, y_err)
+
+    np.testing.assert_allclose(result["par"], [2.0, 0.5], atol=1e-5)
+    assert np.all(result["par_err"] > 0.0)
+    assert result["chi2"][0] == pytest.approx(0.0, abs=1e-6)
+
+
+def test_too_few_unmasked_samples_fail(minimizer):
+    x = np.arange(10.0)
+    y = np.full_like(x, np.nan)
+    y_err = np.zeros_like(x)
+    y[0], y_err[0] = 2.0, 1.0  # one sample left for two free parameters
+    result = aare.Pol1(minimizer=minimizer).fit(x, y, y_err)
+    assert np.all(result["par"] == 0.0)
+    assert np.all(result["chi2"] == 0.0)
+
+
+def test_masked_samples_in_cube_match_1d_fits(minimizer):
+    x = np.linspace(-6.0, 6.0, 61)
+    y = np.stack([gaussian(x, 80.0 + 10.0 * c, 0.8, 1.2) for c in range(4)])
+    y = y.reshape(1, 4, x.size)
+    y_err = np.ones_like(y)
+    for col, masked in [(1, 10), (2, 0), (3, 34)]:
+        y[0, col, masked] = np.nan
+        y_err[0, col, masked] = 0.0
+    model = aare.Gaussian(minimizer=minimizer, compute_errors=True)
+    cube = model.fit(x, y, y_err, n_threads=1)
+
+    for col in range(4):
+        single = model.fit(x, y[0, col], y_err[0, col])
+        np.testing.assert_allclose(cube["par"][0, col], single["par"], atol=1e-12)
+        np.testing.assert_allclose(cube["par_err"][0, col], single["par_err"], atol=1e-12)
+        np.testing.assert_allclose(
+            cube["par"][0, col], [80.0 + 10.0 * col, 0.8, 1.2], rtol=1e-3
+        )
+
+
+def test_parameter_whose_minimum_lies_on_its_limit_reports_no_error(minimizer):
+    # At a minimum that lies exactly on the limit the gradient is zero.
+    x = np.array([-1.0, -1.0, 1.0, 1.0])
+    model = aare.Pol1(minimizer=minimizer, compute_errors=True)
+    model.SetParLimits("p0", 0.0, 1.0)
+    result = model.fit(x, 1.0 + 2.0 * x)
+
+    assert result["par"][0] == pytest.approx(1.0, abs=1e-6)
+    assert result["par_err"][0] == 0.0
+    assert result["par_err"][1] == pytest.approx(0.5, rel=1e-6)
+
+
+@pytest.mark.parametrize("amplitude", [0.1, 1e-4])
+def test_fit_does_not_stop_at_start_without_a_step(minimizer, amplitude):
+    # The EDM tolerance is absolute in chi2 units: an unweighted fit of small
+    # data passes it at the start estimate (sigma 1.3617) already.
+    x = np.linspace(-6.0, 6.0, 61)
+    y = gaussian(x, amplitude, 0.8, 1.3)
+    result = aare.Gaussian(minimizer=minimizer).fit(x, y)
+    assert abs(result["par"][2] - 1.3) < 0.03  # Migrad ends at 1.3205
+
+    one_call = aare.Gaussian(minimizer=minimizer, max_calls=1).fit(x, y)
+    assert np.all(one_call["par"] == 0.0)
+
+
+@pytest.mark.parametrize("lo, span", [(8000, 1000), (20000, 500), (60000, 200)])
+def test_polynomial_far_from_zero_keeps_its_precision(minimizer, lo, span):
+    # Nearly collinear columns 1, x, x^2: VarPro's normal equations lost six
+    # digits at x = 20000.
+    x = np.linspace(lo, lo + span, 40)
+    c = np.array([1.5 - 2e-3 * lo + 3e-7 * lo**2, 2e-3 - 6e-7 * lo, 3e-7])
+    y = c[0] + c[1] * x + c[2] * x**2
+    result = aare.Pol2(minimizer=minimizer).fit(x, y)
+    np.testing.assert_allclose(result["par"], c, rtol=1e-9)

python/tests/test_InterpolationAPI.py

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 import pytest 
 
 from aare import Interpolator, ClusterVector, Etai, Cluster

python/tests/test_Pedestal.py

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 import numpy as np
 import pytest
 
@@ -171,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":
@@ -191,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":

python/tests/test_ROI.py

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 import pytest 
 
 import numpy as np

python/tests/test_RawFile.py

@@ -1,7 +1,7 @@
 # SPDX-License-Identifier: MPL-2.0
 import pytest
 import json
-from aare import File, RawFile, RawSubFile, DetectorType, ROI, UDPPortPosition
+from aare import File, RawFile, RawSubFile, DetectorType, ROI, UDPPortPosition, xy
 import numpy as np
 
 from aare import strixelremap
@@ -31,6 +31,61 @@ def small_raw_file(tmp_path):
     return master_path
 
 
+@pytest.mark.parametrize("reader_type", [RawFile, File])
+@pytest.mark.parametrize("legacy_master", [False, True])
+@pytest.mark.parametrize("padding, policy, supported", [
+    (0, "nodiscard", False),
+    (0, "discard", False),
+    (0, "discardpartial", True),
+    (1, "nodiscard", True),
+    (1, "discard", True),
+    (1, "discardpartial", True),
+])
+def test_raw_frame_policy(small_raw_file, reader_type, legacy_master,
+                          padding, policy, supported):
+    master_path = small_raw_file
+    if legacy_master:
+        master_path = master_path.with_suffix(".raw")
+        master_path.write_text(
+            "Version : 6.4\n"
+            "Detector Type : Jungfrau\n"
+            "Timing Mode : auto\n"
+            "Geometry : [1, 1]\n"
+            "Image Size : 12\n"
+            "Pixels : [3, 2]\n"
+            "Dynamic Range : 16\n"
+            "Max Frames Per File : 1\n"
+            "Total Frames : 2\n"
+            "Frames in File : 2\n"
+            f"Frame Padding : {padding}\n"
+            f"Frame Discard Policy : {policy}\n"
+        )
+    else:
+        metadata = json.loads(master_path.read_text())
+        metadata["Frame Padding"] = padding
+        metadata["Frame Discard Policy"] = policy
+        master_path.write_text(json.dumps(metadata))
+
+    if supported:
+        reader = reader_type(master_path)
+        assert reader.total_frames == 2
+        frame = reader.read_frame()
+        if reader_type is RawFile:
+            _, frame = frame
+        np.testing.assert_array_equal(frame, np.arange(6).reshape(2, 3))
+    else:
+        message = "requires frame padding or discardpartial"
+        with pytest.raises(RuntimeError, match=message) as error:
+            reader_type(master_path)
+        assert str(master_path) in str(error.value)
+
+        for index in range(2):
+            (master_path.parent / f"run_d0_f{index}_0.raw").unlink()
+        with pytest.raises(RuntimeError, match=message) as error:
+            reader_type(master_path)
+        assert str(master_path) in str(error.value)
+
+
 @pytest.mark.parametrize("method, args, kwargs", [
     ("read_frame", (), {}),
     ("read_n", (1,), {}),
@@ -185,6 +240,20 @@ def test_raw_frame_count_across_rois(small_raw_file, short_module):
             reader.read_roi(roi_index)
 
 
+def test_raw_file_geometry(small_raw_file):
+    metadata = json.loads(small_raw_file.read_text())
+    metadata["Geometry"] = {"x": 1, "y": 2}
+    small_raw_file.write_text(json.dumps(metadata))
+    for index in range(2):
+        data = (small_raw_file.parent / f"run_d0_f{index}_0.raw").read_bytes()
+        (small_raw_file.parent / f"run_d1_f{index}_0.raw").write_bytes(data)
+
+    reader = RawFile(small_raw_file)
+    assert reader.geometry == xy(row=2, col=1)
+    assert reader.master.detector_layout == xy(row=2, col=1)
+    assert reader.master.udp_interfaces_per_module == xy(row=1, col=1)
+
+
 @pytest.mark.parametrize("disabled_port", [0, 1])
 def test_raw_frame_count_with_disabled_udp_port(small_raw_file, disabled_port):
     metadata = json.loads(small_raw_file.read_text())

python/tests/test_RawMasterFile.py

@@ -1,9 +1,10 @@
+# SPDX-License-Identifier: MPL-2.0
 
 
 import json
 
 import pytest
-from aare import RawMasterFile, ReadoutMode, DetectorType
+from aare import RawMasterFile, ReadoutMode, DetectorType, FrameDiscardPolicy
 
 
 @pytest.mark.withdata
@@ -16,7 +17,15 @@ def test_read_rawfile_quad_eiger_and_compare_to_numpy(test_data_path):
     assert(f.detector_type == DetectorType.Jungfrau)
 
 
-def test_raw_master_file_context_manager(tmp_path):
+@pytest.mark.parametrize(
+    "policy, expected_policy",
+    [
+        ("nodiscard", FrameDiscardPolicy.NoDiscard),
+        ("discard", FrameDiscardPolicy.Discard),
+        ("discardpartial", FrameDiscardPolicy.DiscardPartial),
+    ],
+)
+def test_raw_master_file_context_manager(tmp_path, policy, expected_policy):
     file_name = tmp_path / "run_master_0.json"
     file_name.write_text(json.dumps({
         "Version": 7.2,
@@ -29,10 +38,13 @@ def test_raw_master_file_context_manager(tmp_path):
         "Total Frames": 2,
         "Frames in File": 2,
         "Frame Padding": 1,
-        "Frame Discard Policy": "nodiscard",
+        "Frame Discard Policy": policy,
     }))
 
     
     with RawMasterFile(file_name) as context_file:
         assert context_file.reading_mode == ReadoutMode.UNKNOWN
         assert context_file.detector_type == DetectorType.Jungfrau
+        frame_discard_policy = context_file.frame_discard_policy
+        assert isinstance(frame_discard_policy, FrameDiscardPolicy)
+        assert frame_discard_policy == expected_policy

python/tests/test_Transformation.py

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 import pytest 
 import numpy as np
 from aare import CtbRawFile, transform

python/tests/test_interpolation_simulated.py

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 
 #test script to test interpolation on simulated data 
 

python/tests/test_make_view.py

@@ -0,0 +1,74 @@
+# SPDX-License-Identifier: MPL-2.0
+"""Probe make_view in np_helper.hpp through the bindings in _aare.testing"""
+import numpy as np
+import pytest
+
+from aare._aare import testing
+
+
+
+
+def test_make_view_roundtrip():
+    """Send an array to C++ and check that we get the same array back"""
+    arr = np.arange(30, dtype=np.float64).reshape(5, 6)
+    out = testing.view_roundtrip_2d(arr)
+    np.testing.assert_array_equal(out, arr)
+
+    #Sanity check to make sure we got a copy an not just the same thing back
+    arr[3,2] = 31655
+    assert arr[3,2] != out[3,2]
+
+def test_make_view_const_roundtrip():
+    """Const works if the array is not writeable"""
+    arr = np.arange(30, dtype=np.float64).reshape(5, 6)
+    arr.flags.writeable = False
+    out = testing.const_view_roundtrip_2d(arr)
+    np.testing.assert_array_equal(out, arr)
+
+def test_make_view_const_roundtrip_writable():
+    """Const works also if the array is  writeable"""
+    arr = np.arange(30, dtype=np.float64).reshape(5, 6)
+    out = testing.const_view_roundtrip_2d(arr)
+    np.testing.assert_array_equal(out, arr)
+
+def test_make_view_fill():
+    arr = np.arange(30, dtype=np.float64).reshape(5, 6)
+    testing.view_fill_2d(arr, 7)
+    assert (arr == 7).all()
+
+def test_make_view_fill_rejects_nonwriteable():
+    arr = np.arange(30, dtype=np.float64).reshape(5, 6)
+    arr.flags.writeable = False
+    with pytest.raises(ValueError, match="not writeable"):
+        testing.view_fill_2d(arr, 7)
+
+def test_make_view_2d_rejects_1d():
+    arr = np.arange(30, dtype=np.uint16)
+    with pytest.raises(ValueError, match="Expected 2D array, got 1D"):
+        testing.view_roundtrip_2d(arr)
+
+def test_make_view_2d_rejects_3d():
+    arr = np.arange(30, dtype=np.uint16).reshape(3, 5, 2)
+    with pytest.raises(ValueError, match="Expected 2D array, got 3D"):
+        testing.view_roundtrip_2d(arr)
+
+def test_make_view_rejects_noncontiguous():
+    arr = np.arange(30, dtype=np.float64).reshape(5, 6)
+    arr = arr[::2, ::3]  # non-contiguous
+    assert not arr.flags.c_contiguous
+    with pytest.raises(ValueError, match="not C-contiguous"):
+        testing.view_roundtrip_2d(arr)
+
+def test_make_view_rejects_noncontiguous_const():
+    arr = np.arange(30, dtype=np.float64).reshape(5, 6)
+    arr = arr[::2, ::3]  # non-contiguous
+    assert not arr.flags.c_contiguous
+    with pytest.raises(ValueError, match="not C-contiguous"):
+        testing.const_view_roundtrip_2d(arr)
+
+def test_make_view_rejects_fortran_order():
+    arr = np.asfortranarray(np.arange(30, dtype=np.float64).reshape(5, 6))
+    assert not arr.flags.c_contiguous
+    with pytest.raises(ValueError, match="not C-contiguous"):
+        testing.view_roundtrip_2d(arr)
+

python/tests/test_xy.py

@@ -0,0 +1,26 @@
+# SPDX-License-Identifier: MPL-2.0
+from aare import xy
+
+
+def test_xy_fields():
+    layout = xy(row=2, col=3)
+    assert layout.row == 2
+    assert layout.col == 3
+
+    layout.row = 4
+    assert layout.row == 4
+
+
+def test_xy_unpacks_as_row_col():
+    row, col = xy(2, 3)
+    assert (row, col) == (2, 3)
+
+
+def test_xy_comparison():
+    assert xy(2, 3) == xy(2, 3)
+    assert xy(2, 3) != xy(3, 2)
+    assert xy(2, 3) != (2, 3)
+
+
+def test_xy_repr():
+    assert repr(xy(2, 3)) == "xy(row=2, col=3)"