Merge remote-tracking branch 'origin/main' into dev/multifit

# Conflicts:
#	RELEASE.md
#	docs/src/pyFit.rst

RELEASE.md

@@ -56,6 +56,14 @@
 
 ### API Changes:
 
+- The C++ ``NDArray`` view constructor and ``NDArray::copy_from()`` take
+  ``NDView<const T, Ndim>``. Mutable views still convert implicitly, so
+  existing callers are unaffected, and read-only views can now be copied.
+- Added C++ overloads of ``adc_sar_05_06_07_08decode64to16``,
+  ``adc_sar_05_decode64to16`` and ``adc_sar_04_decode64to16`` that take the
+  packed samples as ``NDView<const uint8_t, 2>``. Words are assembled with
+  ``memcpy``, so the buffer does not need to be 8-byte aligned. The existing
+  ``NDView<uint64_t, 2>`` overloads are unchanged.
 - ``NDArray`` ``+ - * /`` with a scalar now returns a lazy expression instead
   of an ``NDArray`` and no longer converts the scalar to the element type.
 - ``FitModel`` is now plain data without a pimpl (C++ users need to rebuild);
@@ -114,9 +122,23 @@
 - ``RawMasterFile::rois()`` always returns a list of rois (no optional). Per default it returns a list of one ROI element spawing the entire detector 
 - ``TimingMode::Auto`` changed to ``TimingMode::AUTO_TIMING``, ``TimingMode::Trigger`` changed to ``TimingMode::TRIGGER_EXPOSURE``
 
+
 ### Bugfixes:
 - Fitting a 3D data cube in Python without ``y_err`` now returns ``par_err``
   when ``compute_errors`` is set, as fitting a single pixel already did.
+- The Python ``Cluster`` constructors validate that the data array holds
+  exactly one value per pixel and raise ``ValueError`` otherwise. Previously a
+  longer array wrote past the end of the cluster data.
+- The Python CTB decoding helpers ``adc_sar_*decode64to16``,
+  ``apply_custom_weights``, ``expand24to32bit``, ``expand4to8bit`` and
+  ``decode_my302`` validate their input through ``make_view``. They require
+  C-contiguous arrays of the documented rank and reject other dtypes with
+  ``TypeError`` instead of converting a copy. The ADC SAR decoders no longer
+  reinterpret the byte buffer as aligned 64-bit words; a row length that is
+  not a multiple of 8 bytes raises ``ValueError`` instead of dropping the
+  trailing bytes. Read-only arrays are accepted where the input is only read.
+  ``decode_my302`` raises ``ValueError`` instead of ``RuntimeError`` for a
+  wrong input size.
 - Mismatched operators inhibited vectorization in gcc of NDArray math operators
 - ``NDArray``/``NDView`` math operators and expressions were not vectorized
   for ``uint8_t`` and 64 bit integers, up to 30x slower than a plain loop.
@@ -162,6 +184,10 @@
   without implicit conversion. Both ``push()`` and ``push_with_threshold()``
   validate that frames and thresholds are two-dimensional before constructing
   views, preventing incorrect results from unsupported array layouts or ranks.
+- Python ``ClusterFinder.push_pedestal_frame()`` and ``find_clusters()`` now
+  require C-contiguous ``uint16`` frames without implicit conversion, raising
+  ``TypeError`` instead of silently copying and casting the input. The frame
+  argument can also be passed by keyword as ``frame``.
 
 - Fixed a leaked empty ``ClusterVector`` at the end of Python ``ClusterFile``
   iteration. Chunk iteration now rejects a zero chunk size.

docs/src/NumpyArrays.rst

@@ -0,0 +1,147 @@
+****************
+Numpy arrays
+****************
+
+Most Python functions in aare take NumPy arrays as input. The bindings in
+``python/src/`` receive them as ``py::array_t<T>`` and wrap the buffer in an
+``NDView`` before calling the C++ implementation. No data is copied, so the
+C++ code reads and writes the memory that the caller's array owns. This page
+explains the choices we make when declaring such an argument and what they
+mean for the user.
+
+The preferred pattern
+~~~~~~~~~~~~~~~~~~~~~~~~
+
+Declare the argument with only the element type, take it by value, mark it
+``noconvert`` in the ``py::arg`` list and let ``make_view`` check the rest:
+
+.. code-block:: cpp
+
+    m.def(
+        "push",
+        [](Pedestal &self, py::array_t<uint16_t> frame) {
+            auto view = make_view_2d(frame);
+            self.push(view);
+        },
+        py::arg("frame").noconvert());
+
+The responsibilities are split between pybind11 and ``make_view``:
+
+- pybind11 checks that the argument is a NumPy array with exactly the
+  requested dtype. Anything else fails with a ``TypeError``, and since the
+  argument is ``noconvert`` pybind11 never tries to build a matching array
+  from it.
+- ``make_view`` (``python/src/np_helper.hpp``) checks that the array has the
+  expected number of dimensions and is C-contiguous. It raises a
+  ``ValueError`` with a message that tells the user what to fix, for example
+  ``Array is not C-contiguous. Use np.ascontiguousarray(arr)``. Writing
+  functions use ``make_view`` which also rejects read-only arrays. Read-only
+  input goes through ``make_const_view``.
+
+The result is a binding that either operates on the caller's memory or
+refuses the call. It never silently works on a temporary copy.
+
+Take the array by value
+~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+``py::array_t`` is a handle to the Python object, not the data. Passing it by
+value copies the handle, which increments the reference count of the NumPy
+array, and nothing else. The pixel data stays where it is. We therefore
+declare the parameter as ``py::array_t<T> frame`` rather than
+``py::array_t<T> &frame``:
+
+- A reference saves one reference-count increment per call, which is not
+  measurable next to the work done on a frame.
+- A by-value parameter is a plain lvalue inside the lambda, so
+  ``make_view``, which calls ``mutable_data()``, works without further
+  qualification. A ``const py::array_t<T> &`` does not compile with it.
+- It matches the usual pybind11 style for ``py::object`` and its subclasses,
+  so the binding reads the same as examples in the pybind11 documentation.
+
+Why not let pybind11 convert
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Without ``noconvert``, pybind11 tries every overload once with conversions
+disabled and then once more with conversions enabled. In the second pass a
+``py::array_t<T>`` argument accepts anything NumPy can turn into an array of
+type ``T``: an array with another dtype, a non-contiguous array, or even a
+Python list. pybind11 creates a new array, hands the converted copy to the
+binding and drops it when the call returns. For detector data this has
+several consequences:
+
+- The copy costs time and memory, and the user gets no indication that it
+  happened. For a stack of frames this can be gigabytes.
+- Writes are lost. A function that fills or updates its input, for example
+  a pedestal subtraction in place, writes into the temporary copy.
+- Values are cast without warning. ``forcecast`` allows unsafe casts, so a
+  ``float64`` frame is truncated to ``uint16`` instead of being rejected.
+- Overloads stop working as intended. ``apply_calibration`` is registered for
+  both ``double`` and ``float``. With conversions allowed the first overload
+  would accept everything by casting and the ``float`` overload would never be
+  chosen.
+
+Argument conversion is only worth it for small inputs where convenience beats
+control, for example the parameter vector passed to a fit. Even then the
+function must not write into the argument.
+
+.. note::
+
+    When a binding intentionally lets pybind11 convert, say so with a comment
+    next to the argument. Otherwise a missing ``noconvert`` looks like an
+    oversight and the next person to touch the binding will add it.
+
+    .. code-block:: cpp
+
+        // Conversion allowed: small read-only vector, accept lists and any dtype
+        py::array_t<double, py::array::c_style | py::array::forcecast> par
+
+What the template flags mean
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+``py::array_t<T, Flags>`` takes an optional second template argument. The
+flags do two things: they add requirements that the array must fulfil for the
+*non-converting* check to pass, and they tell NumPy what to produce when a
+conversion *is* allowed.
+
+``py::array::c_style``
+    The array must be C-contiguous. With ``noconvert`` a non-contiguous array
+    is rejected by pybind11 with a ``TypeError`` listing the accepted
+    signatures. Without ``noconvert`` pybind11 makes a contiguous copy.
+    Since ``make_view`` already enforces contiguity and gives a more helpful
+    ``ValueError``, the flag is redundant in the preferred pattern. It is
+    still present in some older bindings and is harmless there.
+
+``py::array::f_style``
+    Same as above for Fortran order. ``NDView`` assumes C order, so we do not
+    use it.
+
+``py::array::forcecast``
+    Allows NumPy to perform unsafe casts, such as ``float64`` to ``uint16``,
+    when converting. It is the **default** value of the template argument, so
+    ``py::array_t<uint16_t>`` and ``py::array_t<uint16_t, py::array::forcecast>``
+    are the same type. The flag has no effect together with ``noconvert``
+    because no conversion takes place. It only matters when conversions are
+    allowed, where it decides whether a lossy cast is accepted or the call is
+    rejected.
+
+The common spelling ``py::array_t<T, py::array::c_style | py::array::forcecast>``
+therefore means: when converting, produce a contiguous array of type ``T`` and
+cast whatever you were given. Combined with ``noconvert`` it collapses to
+"C-contiguous array of exactly type ``T``", which is what the plain
+``py::array_t<T>`` plus ``make_view`` already gives us.
+
+Summary
+~~~~~~~~~~~~~~~~~~
+
+- Declare the argument as ``py::array_t<T>`` by value and add
+  ``.noconvert()`` to its ``py::arg``. The handle is reference counted, so
+  passing it by value does not copy the data.
+- Create the view with ``make_view_1d/2d/3d`` or ``make_const_view_*`` and
+  let it validate dimensions and contiguity. Add an explicit check only when
+  the binding accepts several shapes, as ``apply_calibration`` does for 2D
+  and 3D pedestals.
+- Return data to Python with ``return_image_data`` or ``return_vector`` so that
+  the NumPy array owns the C++ allocation.
+- Only allow conversion for small, read-only inputs where accepting lists or
+  other dtypes is a real convenience for the user, and mark it with a comment
+  in the code.

docs/src/index.rst

@@ -68,3 +68,4 @@ AARE
     Philosophy
     Workflow
     Tests
+    NumpyArrays

docs/src/pyFit.rst

@@ -28,6 +28,16 @@ For three-dimensional data, pass an array with shape
 The result dictionary contains ``par`` and ``chi2``. It also contains
 ``par_err`` when ``compute_errors`` is enabled.
 
+.. note::
+
+    The fit works on ``float64`` data. If ``x``, ``y`` or ``y_err`` is not
+    already a C-contiguous ``float64`` array, the Python bindings convert it
+    before fitting, which allocates a full copy. Integer or ``float32``
+    input therefore works, but for a large three-dimensional ``y`` the
+    copy can be several times the size of the original array. Pass
+    ``float64`` arrays to avoid it. A ``float64`` C-contiguous array is used
+    in place without copying.
+
 Choosing the minimizer
 ~~~~~~~~~~~~~~~~~~~~~~
 

include/aare/NDArray.hpp

@@ -66,14 +66,20 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
 
     /**
      * @brief Construct a new NDArray object from a NDView.
-     * @note The data is copied from the view to the NDArray.
+     * @note The data is copied from the view to the NDArray. The view is
+     * only read, so both NDView<T, Ndim> and NDView<const T, Ndim> are
+     * accepted.
      *
      * @param v view of data to initialize the NDArray with
      */
-    explicit NDArray(const NDView<T, Ndim> v) : NDArray(v.shape()) {
+    explicit NDArray(NDView<const T, Ndim> v) : NDArray(v.shape()) {
         std::copy(v.begin(), v.end(), begin());
     }
 
+    // Keeps an exact match for mutable views so that they do not pick the
+    // element-wise ArrayExpr constructor.
+    explicit NDArray(NDView<T, Ndim> v) : NDArray(NDView<const T, Ndim>(v)) {}
+
     /**
      * @brief Construct a new NDArray object from an std::array.
      */
@@ -149,13 +155,14 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
     /**
      * @brief Copy data from a view of matching shape into this array without
      * reallocating.
-     * @param v view to copy from, must have the same shape as this array
+     * @param v view to copy from, must have the same shape as this array.
+     * Accepts NDView<T, Ndim> and NDView<const T, Ndim>.
      * @throws std::runtime_error if the shapes differ
      *
      * Use this instead of assigning a new NDArray when the buffer needs to be
      * kept, for example when the array is part of a preallocated pool.
      */
-    void copy_from(const NDView<T, Ndim> v) {
+    void copy_from(NDView<const T, Ndim> v) {
         if (v.shape() != shape_) {
             throw std::runtime_error(LOCATION +
                                      "Shape mismatch in NDArray::copy_from");

include/aare/decode.hpp

@@ -16,6 +16,23 @@ void adc_sar_05_decode64to16(NDView<uint64_t, 2> input,
 void adc_sar_04_decode64to16(NDView<uint64_t, 2> input,
                              NDView<uint16_t, 2> output);
 
+/**
+ * @brief Decode packed 64-bit ADC SAR samples stored as raw bytes.
+ *
+ * Every 8 bytes of a row form one 64-bit word in native byte order. Words
+ * are assembled with memcpy, so the buffer does not need to be 8-byte
+ * aligned. The output must have the same number of rows as the input and
+ * input.shape(1) / 8 columns.
+ * @throws std::invalid_argument if the row length is not a multiple of 8
+ * bytes or the output shape does not match.
+ */
+void adc_sar_05_06_07_08decode64to16(NDView<const uint8_t, 2> input,
+                                     NDView<uint16_t, 2> output);
+void adc_sar_05_decode64to16(NDView<const uint8_t, 2> input,
+                             NDView<uint16_t, 2> output);
+void adc_sar_04_decode64to16(NDView<const uint8_t, 2> input,
+                             NDView<uint16_t, 2> output);
+
 /**
  * @brief Called with a 32 bit unsigned integer, shift by offset
  * and then return the lower 24 bits as an 32 bit integer

python/src/bind_Cluster.hpp

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

python/src/bind_ClusterFile.hpp

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

python/src/bind_ClusterFinder.hpp

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

python/src/bind_ClusterFinderMT.hpp

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

python/src/bind_ClusterVector.hpp

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

python/src/bind_FastPedestal.hpp

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

python/src/bind_Interpolator.hpp

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

python/src/bind_Pedestal.hpp

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

python/src/bind_PedestalTrackingPixelHistogram.hpp

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

python/src/bind_PixelHistogram.hpp

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

python/src/bind_calibration.hpp

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

python/src/ctb_raw_file.hpp

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

python/src/fit.hpp

@@ -88,7 +88,8 @@ template <typename Model> void bind_fit_model(py::module &m, const char *name) {
                       &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>
@@ -106,7 +107,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,
@@ -195,11 +196,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);
@@ -255,11 +251,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 {
@@ -349,7 +340,8 @@ void define_fit_bindings(py::module &m) {
 
     m.def(
         "fit",
-        [](py::object model_obj,
+        [](py::object model_obj, // conversion ok, we want to be able to call
+                                 // with any dtype
            py::array_t<double, py::array::c_style | py::array::forcecast> x,
            py::array_t<double, py::array::c_style | py::array::forcecast> y,
            py::object y_err_obj, int n_threads) -> py::object {

python/src/var_cluster.hpp

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

python/tests/test_Cluster.py

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

python/tests/test_ClusterVector.py

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

python/tests/test_Pedestal.py

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

src/NDArray.test.cpp

@@ -86,6 +86,40 @@ TEST_CASE("A const NDArray returns a read-only view") {
     REQUIRE(view(1, 2) == 1);
 }
 
+TEST_CASE("NDArray can be constructed from a read-only view") {
+    NDArray<int, 2> source({2, 3}, 0);
+    for (ssize_t i = 0; i < source.size(); ++i)
+        source(i) = static_cast<int>(i);
+    const auto &const_source = source;
+
+    NDArray<int, 2> from_const(const_source.view());
+    NDArray<int, 2> from_mutable(source.view());
+
+    REQUIRE(from_const.shape() == source.shape());
+    REQUIRE(from_mutable.shape() == source.shape());
+    for (ssize_t i = 0; i < source.size(); ++i) {
+        REQUIRE(from_const(i) == source(i));
+        REQUIRE(from_mutable(i) == source(i));
+    }
+
+    // Independent copy, not a view of the source
+    source(0) = 42;
+    REQUIRE(from_const(0) == 0);
+}
+
+TEST_CASE("NDArray::copy_from accepts a read-only view") {
+    NDArray<double, 1> source({4}, 1.5);
+    const auto &const_source = source;
+
+    NDArray<double, 1> target({4}, 0.0);
+    target.copy_from(const_source.view());
+    for (ssize_t i = 0; i < target.size(); ++i)
+        REQUIRE(target(i) == 1.5);
+
+    NDArray<double, 1> wrong_shape({3}, 0.0);
+    REQUIRE_THROWS(wrong_shape.copy_from(const_source.view()));
+}
+
 TEST_CASE("Indexing of a 2D image") {
     std::array<ssize_t, 2> shape{{3, 7}};
     NDArray<long> img(shape, 5);

src/decode.cpp

@@ -1,6 +1,7 @@
 // SPDX-License-Identifier: MPL-2.0
 #include "aare/decode.hpp"
 #include <cmath>
+#include <cstring>
 #include <fmt/format.h>
 namespace aare {
 
@@ -102,6 +103,57 @@ void adc_sar_04_decode64to16(NDView<uint64_t, 2> input,
     }
 }
 
+namespace {
+
+constexpr ssize_t bytes_per_word = sizeof(uint64_t);
+
+// Shared loop for the byte-buffer overloads. Words are copied out with
+// memcpy so the input does not need to be aligned to 8 bytes.
+template <uint16_t (*decode_word)(uint64_t)>
+void decode64to16_from_bytes(NDView<const uint8_t, 2> input,
+                             NDView<uint16_t, 2> output) {
+    if (input.shape(1) % bytes_per_word != 0) {
+        throw std::invalid_argument(
+            fmt::format("{} input row length of {} bytes is not a multiple "
+                        "of {}",
+                        LOCATION, input.shape(1), bytes_per_word));
+    }
+    if (input.shape(0) != output.shape(0) ||
+        input.shape(1) / bytes_per_word != output.shape(1)) {
+        throw std::invalid_argument(fmt::format(
+            "{} output shape must be ({}, {}) for an input of "
+            "shape ({}, {})",
+            LOCATION, input.shape(0), input.shape(1) / bytes_per_word,
+            output.shape(0), output.shape(1)));
+    }
+
+    for (ssize_t i = 0; i < output.shape(0); ++i) {
+        const uint8_t *row = input.data() + i * input.shape(1);
+        for (ssize_t j = 0; j < output.shape(1); ++j) {
+            uint64_t word{};
+            std::memcpy(&word, row + j * bytes_per_word, sizeof(word));
+            output(i, j) = decode_word(word);
+        }
+    }
+}
+
+} // namespace
+
+void adc_sar_05_06_07_08decode64to16(NDView<const uint8_t, 2> input,
+                                     NDView<uint16_t, 2> output) {
+    decode64to16_from_bytes<adc_sar_05_06_07_08decode64to16>(input, output);
+}
+
+void adc_sar_05_decode64to16(NDView<const uint8_t, 2> input,
+                             NDView<uint16_t, 2> output) {
+    decode64to16_from_bytes<adc_sar_05_decode64to16>(input, output);
+}
+
+void adc_sar_04_decode64to16(NDView<const uint8_t, 2> input,
+                             NDView<uint16_t, 2> output) {
+    decode64to16_from_bytes<adc_sar_04_decode64to16>(input, output);
+}
+
 double apply_custom_weights(uint16_t input, const NDView<double, 1> weights) {
     if (weights.size() > 16) {
         throw std::invalid_argument(

src/decode.test.cpp

@@ -5,6 +5,8 @@
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/matchers/catch_matchers_floating_point.hpp>
 using Catch::Matchers::WithinAbs;
+#include <cstring>
+#include <stdexcept>
 #include <vector>
 
 using aare::BitOffset;
@@ -49,6 +51,86 @@ TEST_CASE("test_adc_sar_05_decode64to16") {
     CHECK(output == 3072UL);
 }
 
+TEST_CASE("Decode 64 bit ADC SAR samples from a byte buffer") {
+    // Two rows of two words each, viewed as raw bytes
+    aare::NDArray<uint64_t, 2> words({2, 2});
+    words(0, 0) = 1UL << 29;
+    words(0, 1) = (1UL << 25) | (1UL << 22);
+    words(1, 0) = 0;
+    words(1, 1) = (1UL << 29) | (1UL << 19) | (1UL << 28);
+
+    aare::NDView<const uint8_t, 2> bytes(
+        reinterpret_cast<const uint8_t *>(words.data()), {2, 16});
+    aare::NDArray<uint16_t, 2> out({2, 2});
+    aare::adc_sar_05_decode64to16(bytes, out.view());
+
+    CHECK(out(0, 0) == 1);
+    CHECK(out(0, 1) == 3072);
+    CHECK(out(1, 0) == 0);
+    CHECK(out(1, 1) == 7);
+
+    // Matches the overload that takes the words directly
+    aare::NDArray<uint16_t, 2> reference({2, 2});
+    aare::adc_sar_05_decode64to16(words.view(), reference.view());
+    for (ssize_t i = 0; i < 2; ++i) {
+        for (ssize_t j = 0; j < 2; ++j) {
+            CHECK(out(i, j) == reference(i, j));
+        }
+    }
+}
+
+TEST_CASE("Each byte buffer decoder maps to its own bit order") {
+    // First two bits of each decoder's bit list
+    struct Case {
+        void (*decode)(aare::NDView<const uint8_t, 2>,
+                       aare::NDView<uint16_t, 2>);
+        int bit0;
+        int bit1;
+    };
+    const Case cases[] = {
+        {aare::adc_sar_05_06_07_08decode64to16, 29, 17},
+        {aare::adc_sar_05_decode64to16, 29, 19},
+        {aare::adc_sar_04_decode64to16, 15, 17},
+    };
+
+    for (const auto &c : cases) {
+        uint64_t word = (1UL << c.bit0) | (1UL << c.bit1);
+        std::vector<uint8_t> buffer(8);
+        std::memcpy(buffer.data(), &word, sizeof(word));
+        aare::NDView<const uint8_t, 2> bytes(buffer.data(), {1, 8});
+        aare::NDArray<uint16_t, 2> out({1, 1});
+        c.decode(bytes, out.view());
+        CHECK(out(0, 0) == 3);
+    }
+}
+
+TEST_CASE("Byte buffer decoding does not require an aligned buffer") {
+    // Place the word three bytes into the buffer
+    std::vector<uint8_t> buffer(3 + 8, 0);
+    uint64_t word = (1UL << 29) | (1UL << 19);
+    std::memcpy(buffer.data() + 3, &word, sizeof(word));
+
+    aare::NDView<const uint8_t, 2> bytes(buffer.data() + 3, {1, 8});
+    aare::NDArray<uint16_t, 2> out({1, 1});
+    aare::adc_sar_05_decode64to16(bytes, out.view());
+    CHECK(out(0, 0) == 3);
+}
+
+TEST_CASE("Byte buffer decoding rejects invalid shapes") {
+    std::vector<uint8_t> buffer(24, 0);
+    aare::NDArray<uint16_t, 2> out({1, 1});
+
+    // Row length that is not a multiple of 8
+    aare::NDView<const uint8_t, 2> odd_row(buffer.data(), {2, 12});
+    REQUIRE_THROWS_AS(aare::adc_sar_05_decode64to16(odd_row, out.view()),
+                      std::invalid_argument);
+
+    // Output shape that does not match rows x (bytes / 8)
+    aare::NDView<const uint8_t, 2> two_words(buffer.data(), {1, 16});
+    REQUIRE_THROWS_AS(aare::adc_sar_05_decode64to16(two_words, out.view()),
+                      std::invalid_argument);
+}
+
 TEST_CASE("test_apply_custom_weights") {
 
     uint16_t input = 1;