From 941325ea151e627f61c9e732e34a751ecfffa510 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Erik=20Fr=C3=B6jdh?= Date: Fri, 25 Sep 2026 14:40:33 +0200 Subject: [PATCH] More uniform binding of Numpy arrays (#381) - Introduce a standard way to bind numpy arrays: py::array_t and .noconvert() on the argument - Take py::array_t by value since it's just a small wrapper - Update bindings --- RELEASE.md | 26 ++ docs/src/NumpyArrays.rst | 147 ++++++++++ docs/src/index.rst | 1 + docs/src/pyFit.rst | 10 + include/aare/NDArray.hpp | 15 +- include/aare/decode.hpp | 17 ++ python/src/bind_Cluster.hpp | 11 +- python/src/bind_ClusterFile.hpp | 10 +- python/src/bind_ClusterFinder.hpp | 16 +- python/src/bind_ClusterFinderMT.hpp | 16 +- python/src/bind_ClusterVector.hpp | 3 +- python/src/bind_FastPedestal.hpp | 6 +- python/src/bind_Interpolator.hpp | 19 +- python/src/bind_Pedestal.hpp | 17 +- .../bind_PedestalTrackingPixelHistogram.hpp | 6 +- python/src/bind_PixelHistogram.hpp | 2 +- python/src/bind_calibration.hpp | 32 +-- python/src/ctb_raw_file.hpp | 254 +++++++----------- python/src/fit.hpp | 18 +- python/src/var_cluster.hpp | 58 ++-- python/tests/test_Cluster.py | 18 ++ python/tests/test_ClusterVector.py | 2 +- python/tests/test_Pedestal.py | 5 +- src/NDArray.test.cpp | 34 +++ src/decode.cpp | 52 ++++ src/decode.test.cpp | 82 ++++++ 26 files changed, 597 insertions(+), 280 deletions(-) create mode 100644 docs/src/NumpyArrays.rst diff --git a/RELEASE.md b/RELEASE.md index 256cbb63..6cf4c316 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -35,6 +35,14 @@ ### API Changes: +- The C++ ``NDArray`` view constructor and ``NDArray::copy_from()`` take + ``NDView``. 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``. Words are assembled with + ``memcpy``, so the buffer does not need to be 8-byte aligned. The existing + ``NDView`` 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. - ``FastPedestal`` variance is now a private ``double`` intermediate. Removed @@ -80,7 +88,21 @@ - ``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: +- 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. @@ -126,6 +148,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. diff --git a/docs/src/NumpyArrays.rst b/docs/src/NumpyArrays.rst new file mode 100644 index 00000000..b837771c --- /dev/null +++ b/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`` 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 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 frame`` rather than +``py::array_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 &`` 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`` 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 par + +What the template flags mean +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +``py::array_t`` 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`` and ``py::array_t`` + 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`` +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`` plus ``make_view`` already gives us. + +Summary +~~~~~~~~~~~~~~~~~~ + +- Declare the argument as ``py::array_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. diff --git a/docs/src/index.rst b/docs/src/index.rst index ab22e57d..4f0151b7 100644 --- a/docs/src/index.rst +++ b/docs/src/index.rst @@ -68,3 +68,4 @@ AARE Philosophy Workflow Tests + NumpyArrays diff --git a/docs/src/pyFit.rst b/docs/src/pyFit.rst index 257f7414..a8375772 100644 --- a/docs/src/pyFit.rst +++ b/docs/src/pyFit.rst @@ -27,4 +27,14 @@ 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. + .. autofunction:: fit diff --git a/include/aare/NDArray.hpp b/include/aare/NDArray.hpp index 68f88ad6..3a2fbda3 100644 --- a/include/aare/NDArray.hpp +++ b/include/aare/NDArray.hpp @@ -66,14 +66,20 @@ class NDArray : public ArrayExpr, 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 and NDView are + * accepted. * * @param v view of data to initialize the NDArray with */ - explicit NDArray(const NDView v) : NDArray(v.shape()) { + explicit NDArray(NDView 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 v) : NDArray(NDView(v)) {} + /** * @brief Construct a new NDArray object from an std::array. */ @@ -149,13 +155,14 @@ class NDArray : public ArrayExpr, 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 and NDView. * @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 v) { + void copy_from(NDView v) { if (v.shape() != shape_) { throw std::runtime_error(LOCATION + "Shape mismatch in NDArray::copy_from"); diff --git a/include/aare/decode.hpp b/include/aare/decode.hpp index a6efd187..00433109 100644 --- a/include/aare/decode.hpp +++ b/include/aare/decode.hpp @@ -16,6 +16,23 @@ void adc_sar_05_decode64to16(NDView input, void adc_sar_04_decode64to16(NDView input, NDView 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 input, + NDView output); +void adc_sar_05_decode64to16(NDView input, + NDView output); +void adc_sar_04_decode64to16(NDView input, + NDView output); + /** * @brief Called with a 32 bit unsigned integer, shift by offset * and then return the lower 24 bits as an 32 bit integer diff --git a/python/src/bind_Cluster.hpp b/python/src/bind_Cluster.hpp index 5a749bfc..afa9bfd1 100644 --- a/python/src/bind_Cluster.hpp +++ b/python/src/bind_Cluster.hpp @@ -24,14 +24,21 @@ void define_Cluster(py::module &m, const std::string &typestr) { py::class_>( 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 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 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; diff --git a/python/src/bind_ClusterFile.hpp b/python/src/bind_ClusterFile.hpp index 09fc649b..c8544e9d 100644 --- a/python/src/bind_ClusterFile.hpp +++ b/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 &self, - py::array_t noise_map) { + py::array_t + 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 &self, - py::array_t gain_map) { + py::array_t + gain_map) { auto view = make_view_2d(gain_map); self.set_gain_map(view); }, diff --git a/python/src/bind_ClusterFinder.hpp b/python/src/bind_ClusterFinder.hpp index 2e0bf8b3..661b9d75 100644 --- a/python/src/bind_ClusterFinder.hpp +++ b/python/src/bind_ClusterFinder.hpp @@ -43,12 +43,14 @@ void define_ClusterFinder(py::module &m, const std::string &typestr) { &ClusterFinder::set_nSigma, R"(number of sigma above the pedestal to consider a photon during cluster finding.)") - .def("push_pedestal_frame", - [](ClusterFinder &self, - py::array_t frame) { - auto view = make_view_2d(frame); - self.push_pedestal_frame(view); - }) + .def( + "push_pedestal_frame", + [](ClusterFinder &self, + py::array_t frame) { + auto view = make_view_2d(frame); + self.push_pedestal_frame(view); + }, + py::arg("frame").noconvert()) .def("clear_pedestal", &ClusterFinder::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 \ No newline at end of file diff --git a/python/src/bind_ClusterFinderMT.hpp b/python/src/bind_ClusterFinderMT.hpp index 37ab8d47..2c166a1a 100644 --- a/python/src/bind_ClusterFinderMT.hpp +++ b/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 &self, - py::array_t frame) { - auto view = make_view_2d(frame); - self.push_pedestal_frame(view); - }) + .def( + "push_pedestal_frame", + [](ClusterFinderMT &self, + py::array_t frame) { + auto view = make_view_2d(frame); + self.push_pedestal_frame(view); + }, + py::arg("frame").noconvert()) .def( "find_clusters", [](ClusterFinderMT &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 &self) { diff --git a/python/src/bind_ClusterVector.hpp b/python/src/bind_ClusterVector.hpp index c6888a4a..9b031e5d 100644 --- a/python/src/bind_ClusterVector.hpp +++ b/python/src/bind_ClusterVector.hpp @@ -46,8 +46,7 @@ void define_ClusterVector(py::module &m, const std::string &typestr) { .def( "__call__", - [](ClusterVector &self, - py::array_t mask) { + [](ClusterVector &self, py::array_t mask) { return self(make_view_1d(mask)); }, py::arg("mask").noconvert(), R"doc( diff --git a/python/src/bind_FastPedestal.hpp b/python/src/bind_FastPedestal.hpp index b2b4fb8f..5c278490 100644 --- a/python/src/bind_FastPedestal.hpp +++ b/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 &pedestal, - py::array_t &frame) { + [](FastPedestal &pedestal, py::array_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 &pedestal, - py::array_t &frame) { + [](FastPedestal &pedestal, py::array_t frame) { pedestal.add_init_frame(make_view_2d(frame)); }, py::arg("frame").noconvert(), diff --git a/python/src/bind_Interpolator.hpp b/python/src/bind_Interpolator.hpp index a98e539d..3b7a015b 100644 --- a/python/src/bind_Interpolator.hpp +++ b/python/src/bind_Interpolator.hpp @@ -101,9 +101,7 @@ void define_interpolation_bindings(py::module &m) { auto interpolator = py::class_(m, "Interpolator") .def(py::init( - [](py::array_t - etacube, + [](py::array_t etacube, py::array_t xbins, py::array_t ybins, py::array_t 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 xbins, py::array_t 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 + py::array_t 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{}; diff --git a/python/src/bind_Pedestal.hpp b/python/src/bind_Pedestal.hpp index 602ee13e..9e8af455 100644 --- a/python/src/bind_Pedestal.hpp +++ b/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 &pedestal, - py::array_t &f) { - if (f.ndim() != 2) { - throw py::value_error("Frame must be 2-dimensional"); - } + [](Pedestal &pedestal, py::array_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 &pedestal, - py::array_t &f, - py::array_t &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 &pedestal, py::array_t f, + py::array_t threshold) { auto frame_view = make_view_2d(f); auto threshold_view = make_view_2d(threshold); pedestal.push_with_threshold(frame_view, threshold_view); diff --git a/python/src/bind_PedestalTrackingPixelHistogram.hpp b/python/src/bind_PedestalTrackingPixelHistogram.hpp index 35b8b24a..139b58a1 100644 --- a/python/src/bind_PedestalTrackingPixelHistogram.hpp +++ b/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 - frame) { + py::array_t 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 - image) { + py::array_t 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. diff --git a/python/src/bind_PixelHistogram.hpp b/python/src/bind_PixelHistogram.hpp index d205bce2..d3a1b462 100644 --- a/python/src/bind_PixelHistogram.hpp +++ b/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 image) { + [](Hist &self, py::array_t 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. diff --git a/python/src/bind_calibration.hpp b/python/src/bind_calibration.hpp index c9479e8a..6efa87bc 100644 --- a/python/src/bind_calibration.hpp +++ b/python/src/bind_calibration.hpp @@ -11,38 +11,34 @@ namespace py = pybind11; template -py::array_t pybind_apply_calibration( - py::array_t data, - py::array_t pedestal, - py::array_t - calibration, - int n_threads = 4) { +py::array_t +pybind_apply_calibration(py::array_t data, + py::array_t pedestal, + py::array_t 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(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(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(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 pybind_count_switching_pixels( - py::array_t data, - ssize_t n_threads = 4) { +py::array_t pybind_count_switching_pixels(py::array_t data, + ssize_t n_threads = 4) { auto data_span = make_view_3d(data); auto arr = new NDArray{}; @@ -51,9 +47,8 @@ py::array_t pybind_count_switching_pixels( } template -py::array_t pybind_calculate_pedestal( - py::array_t data, - ssize_t n_threads) { +py::array_t pybind_calculate_pedestal(py::array_t data, + ssize_t n_threads) { auto data_span = make_view_3d(data); auto arr = new NDArray{}; @@ -62,9 +57,8 @@ py::array_t pybind_calculate_pedestal( } template -py::array_t pybind_calculate_pedestal_g0( - py::array_t data, - ssize_t n_threads) { +py::array_t pybind_calculate_pedestal_g0(py::array_t data, + ssize_t n_threads) { auto data_span = make_view_3d(data); auto arr = new NDArray{}; diff --git a/python/src/ctb_raw_file.hpp b/python/src/ctb_raw_file.hpp index 5a6953b7..b6338bdd 100644 --- a/python/src/ctb_raw_file.hpp +++ b/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 , NDView)> +py::array_t decode64to16(py::array_t input) { + auto bytes = make_const_view_2d(input); + constexpr ssize_t bytes_per_word = sizeof(uint64_t); + auto output = py::array_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 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, + py::arg("input").noconvert(), decode64to16_doc); - // Create a 2D output array with the same shape as the input - std::vector shape{input.shape(0), - input.shape(1) / - static_cast(bits_per_byte)}; - py::array_t output(shape); + m.def("adc_sar_05_decode64to16", &decode64to16, + py::arg("input").noconvert(), decode64to16_doc); - // Create a view of the input and output arrays - NDView input_view( - reinterpret_cast(input.mutable_data()), - {output.shape(0), output.shape(1)}); - NDView output_view(output.mutable_data(), - {output.shape(0), output.shape(1)}); + m.def("adc_sar_04_decode64to16", &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 input, py::array_t weights) { + auto input_view = make_view_1d(input); + auto weights_view = make_view_1d(weights); + auto output = py::array_t(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 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(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 input) { - if (input.ndim() != 2) { - throw std::runtime_error( - "Only 2D arrays are supported at this moment"); - } + m.def( + "expand4to8bit", + [](py::array_t input) { + auto input_view = make_const_view_1d(input); + auto output = py::array_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 shape{input.shape(0), - input.shape(1) / - static_cast(bits_per_byte)}; - py::array_t output(shape); + m.def( + "decode_my302", + [](py::array_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 input_view( - reinterpret_cast(input.mutable_data()), - {output.shape(0), output.shape(1)}); - NDView 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 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 shape{input.shape(0), - input.shape(1) / - static_cast(bits_per_byte)}; - py::array_t output(shape); + auto output = py::array_t(channels * counters); + auto output_view = make_view_1d(output); - // Create a view of the input and output arrays - NDView input_view( - reinterpret_cast(input.mutable_data()), - {output.shape(0), output.shape(1)}); - NDView 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 in_part( + input_view.data() + in_step * i, {in_step}); + NDView 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 - &input, - py::array_t - &weights) { - // Create new array with same shape as the input array - // (uninitialized values) - py::buffer_info buf = input.request(); - py::array_t output(buf.shape); - - // Use NDViews to call into the C++ library - auto weights_view = make_view_1d(weights); - NDView input_view(input.mutable_data(), - {input.size()}); - NDView output_view(output.mutable_data(), - {output.size()}); - - apply_custom_weights(input_view, output_view, weights_view); - return output; - }); - - m.def("expand24to32bit", - [](py::array_t - &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 output(buf.size / bytes_per_channel); - - NDView input_view(input.mutable_data(), - {input.size()}); - NDView output_view(output.mutable_data(), - {output.size()}); - - aare::expand24to32bit(input_view, output_view, bitoff); - return output; - }); - - m.def("expand4to8bit", - [](py::array_t - &input) { - py::buffer_info buf = input.request(); - - py::array_t output(buf.size * 2); - - NDView input_view(input.mutable_data(), - {input.size()}); - NDView output_view(output.mutable_data(), - {output.size()}); - - aare::expand4to8bit(input_view, output_view); - return output; - }); - - m.def("decode_my302", - [](py::array_t - &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 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 input_view(input.mutable_data() + step * i, - {input.size() / n_outputs}); - NDView 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_(m, "CtbRawFile") .def(py::init()) diff --git a/python/src/fit.hpp b/python/src/fit.hpp index 82bff325..1e8d9c39 100644 --- a/python/src/fit.hpp +++ b/python/src/fit.hpp @@ -60,7 +60,8 @@ template void bind_fit_model(py::module &m, const char *name) { .def_property("compute_errors", &FM::compute_errors, &FM::SetComputeErrors) .def( - "__call__", + "__call__", // conversion ok, we want to be able to call with any + // dtype [](const FM & /*self*/, py::array_t x, py::array_t @@ -78,7 +79,7 @@ template 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 x, py::array_t y, @@ -163,11 +164,6 @@ fit_dispatch(const aare::FitModel &model, py::array_t>( 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({y.shape(0), y.shape(1), npar}, 0.0); auto y_view_err = make_view_3d(y_err); @@ -208,11 +204,6 @@ fit_dispatch(const aare::FitModel &model, py::array_t>( 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, x_view, y_view, y_view_err); } else { @@ -241,7 +232,8 @@ void define_fit_bindings(py::module &m) { m.def( "fit", - [](py::object model_obj, + [](py::object model_obj, // conversion ok, we want to be able to call + // with any dtype py::array_t x, py::array_t y, py::object y_err_obj, int n_threads) -> py::object { diff --git a/python/src/var_cluster.hpp b/python/src/var_cluster.hpp index 44d70db1..1a3b674c 100644 --- a/python/src/var_cluster.hpp +++ b/python/src/var_cluster.hpp @@ -21,40 +21,42 @@ void define_var_cluster_finder_bindings(py::module &m) { auto *ptr = new NDArray(self.labeled()); return return_image_data(ptr); }) - .def("set_noiseMap", - [](VarClusterFinder &self, - py::array_t - noise_map) { - auto noise_map_span = make_view_2d(noise_map); - self.set_noiseMap(noise_map_span); - }) + .def( + "set_noiseMap", // conversions ok + [](VarClusterFinder &self, + py::array_t + 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::set_numberOfNeighbours) .def("set_empty_surroundingPixels", &VarClusterFinder::set_empty_surroundingPixels) .def("set_peripheralThresholdFactor", &VarClusterFinder::set_peripheralThresholdFactor) - .def("find_clusters", - [](VarClusterFinder &self, - py::array_t - img) { - auto view = make_view_2d(img); - self.find_clusters(view); - }) - .def("find_clusters_X", - [](VarClusterFinder &self, - py::array_t - img) { - auto img_span = make_view_2d(img); - self.find_clusters_X(img_span); - }) - .def("single_pass", - [](VarClusterFinder &self, - py::array_t - img) { - auto img_span = make_view_2d(img); - self.single_pass(img_span); - }) + .def( + "find_clusters", + [](VarClusterFinder &self, py::array_t img) { + auto view = make_view_2d(img); + self.find_clusters(view); + }, + py::arg("img").noconvert()) + .def( + "find_clusters_X", + [](VarClusterFinder &self, py::array_t img) { + auto img_span = make_view_2d(img); + self.find_clusters_X(img_span); + }, + py::arg("img").noconvert()) + .def( + "single_pass", + [](VarClusterFinder &self, py::array_t img) { + auto img_span = make_view_2d(img); + self.single_pass(img_span); + }, + py::arg("img").noconvert()) .def("hits", [](VarClusterFinder &self) { auto ptr = new std::vector::Hit>( diff --git a/python/tests/test_Cluster.py b/python/tests/test_Cluster.py index 027cc203..3978ab0c 100644 --- a/python/tests/test_Cluster.py +++ b/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() diff --git a/python/tests/test_ClusterVector.py b/python/tests/test_ClusterVector.py index dd75366e..b0e9cffa 100644 --- a/python/tests/test_ClusterVector.py +++ b/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) diff --git a/python/tests/test_Pedestal.py b/python/tests/test_Pedestal.py index 4248d157..d2a3ddb3 100644 --- a/python/tests/test_Pedestal.py +++ b/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": diff --git a/src/NDArray.test.cpp b/src/NDArray.test.cpp index 3657c7b0..872f22cd 100644 --- a/src/NDArray.test.cpp +++ b/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 source({2, 3}, 0); + for (ssize_t i = 0; i < source.size(); ++i) + source(i) = static_cast(i); + const auto &const_source = source; + + NDArray from_const(const_source.view()); + NDArray 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 source({4}, 1.5); + const auto &const_source = source; + + NDArray 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 wrong_shape({3}, 0.0); + REQUIRE_THROWS(wrong_shape.copy_from(const_source.view())); +} + TEST_CASE("Indexing of a 2D image") { std::array shape{{3, 7}}; NDArray img(shape, 5); diff --git a/src/decode.cpp b/src/decode.cpp index afe261f9..f522f5c1 100644 --- a/src/decode.cpp +++ b/src/decode.cpp @@ -1,6 +1,7 @@ // SPDX-License-Identifier: MPL-2.0 #include "aare/decode.hpp" #include +#include #include namespace aare { @@ -102,6 +103,57 @@ void adc_sar_04_decode64to16(NDView 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 +void decode64to16_from_bytes(NDView input, + NDView 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 input, + NDView output) { + decode64to16_from_bytes(input, output); +} + +void adc_sar_05_decode64to16(NDView input, + NDView output) { + decode64to16_from_bytes(input, output); +} + +void adc_sar_04_decode64to16(NDView input, + NDView output) { + decode64to16_from_bytes(input, output); +} + double apply_custom_weights(uint16_t input, const NDView weights) { if (weights.size() > 16) { throw std::invalid_argument( diff --git a/src/decode.test.cpp b/src/decode.test.cpp index bab7c38b..09d6c61c 100644 --- a/src/decode.test.cpp +++ b/src/decode.test.cpp @@ -5,6 +5,8 @@ #include #include using Catch::Matchers::WithinAbs; +#include +#include #include 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 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 bytes( + reinterpret_cast(words.data()), {2, 16}); + aare::NDArray 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 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, + aare::NDView); + 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 buffer(8); + std::memcpy(buffer.data(), &word, sizeof(word)); + aare::NDView bytes(buffer.data(), {1, 8}); + aare::NDArray 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 buffer(3 + 8, 0); + uint64_t word = (1UL << 29) | (1UL << 19); + std::memcpy(buffer.data() + 3, &word, sizeof(word)); + + aare::NDView bytes(buffer.data() + 3, {1, 8}); + aare::NDArray 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 buffer(24, 0); + aare::NDArray out({1, 1}); + + // Row length that is not a multiple of 8 + aare::NDView 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 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;