added template for making views

python/src/bind_ClusterVector.hpp

@@ -48,10 +48,7 @@ void define_ClusterVector(py::module &m, const std::string &typestr) {
             "__call__",
             [](ClusterVector<ClusterType> &self,
                py::array_t<bool, py::array::c_style> mask) {
-                if (mask.ndim() != 1) {
-                    throw py::value_error("Mask must be one-dimensional");
-                }
-                return self(make_view_1d(mask));
+                return self(make_view<1>(mask));
             },
             py::arg("mask").noconvert(), R"doc(
             Return a filtered copy of this ClusterVector.

python/src/bind_testing.hpp

@@ -0,0 +1,145 @@
+// SPDX-License-Identifier: MPL-2.0
+// Probes used by the python test-suite to exercise the helpers in
+// np_helper.hpp. Internal, no API guarantees.
+#pragma once
+
+#include "aare/NDArray.hpp"
+#include "aare/NDView.hpp"
+#include "np_helper.hpp"
+
+#include <algorithm>
+#include <cstdint>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <type_traits>
+
+namespace py = pybind11;
+using namespace ::aare;
+
+// Arrays are taken with the default flags (no c_style) so that nothing is
+// made contiguous by pybind11 before it reaches make_view
+template <typename T, ssize_t Ndim>
+void define_make_view_probes(py::module &m, const std::string &suffix) {
+
+    // Copy what C++ sees through the view back to python
+    m.def(("view_roundtrip_" + suffix).c_str(),
+          [](py::array_t<T> arr) {
+              auto view = make_view<Ndim>(arr);
+              return return_image_data(new NDArray<T, Ndim>(view));
+          },
+          py::arg("arr").noconvert());
+
+    // shape, strides (in elements) and address of the data seen by the view
+    m.def(("view_info_" + suffix).c_str(),
+          [](py::array_t<T> arr) {
+              auto view = make_view<Ndim>(arr);
+              return py::make_tuple(view.shape(), view.strides(),
+                                    reinterpret_cast<uintptr_t>(view.data()));
+          },
+          py::arg("arr").noconvert());
+
+    // Same as above but through a read only view
+    m.def(("const_view_roundtrip_" + suffix).c_str(),
+          [](py::array_t<T> arr) {
+              auto view = make_const_view<Ndim>(arr);
+              static_assert(
+                  std::is_same_v<decltype(view), NDView<const T, Ndim>>);
+              auto out = new NDArray<T, Ndim>(view.shape());
+              std::copy(view.begin(), view.end(), out->begin());
+              return return_image_data(out);
+          },
+          py::arg("arr").noconvert());
+
+    m.def(("const_view_info_" + suffix).c_str(),
+          [](py::array_t<T> arr) {
+              auto view = make_const_view<Ndim>(arr);
+              return py::make_tuple(view.shape(), view.strides(),
+                                    reinterpret_cast<uintptr_t>(view.data()));
+          },
+          py::arg("arr").noconvert());
+
+    // Write through the view
+    m.def(("view_fill_" + suffix).c_str(),
+          [](py::array_t<T> arr, T value) {
+              auto view = make_view<Ndim>(arr);
+              view = value;
+          },
+          py::arg("arr").noconvert(), py::arg("value"));
+}
+
+void define_testing_bindings(py::module &m) {
+    auto testing = m.def_submodule(
+        "testing", "Internal helpers for the test-suite, no API guarantees");
+
+    // noconvert + overloads gives dispatch on dtype
+    define_make_view_probes<double, 1>(testing, "1d");
+    define_make_view_probes<double, 2>(testing, "2d");
+    define_make_view_probes<double, 3>(testing, "3d");
+    define_make_view_probes<uint16_t, 1>(testing, "1d");
+    define_make_view_probes<uint16_t, 2>(testing, "2d");
+    define_make_view_probes<uint16_t, 3>(testing, "3d");
+
+    // The wrappers used by the bindings
+    testing.def(
+        "make_view_1d_info",
+        [](py::array_t<double> arr) {
+            auto view = make_view_1d(arr);
+            return py::make_tuple(view.shape(), view.strides());
+        },
+        py::arg("arr").noconvert());
+    testing.def(
+        "make_view_2d_info",
+        [](py::array_t<double> arr) {
+            auto view = make_view_2d(arr);
+            return py::make_tuple(view.shape(), view.strides());
+        },
+        py::arg("arr").noconvert());
+    testing.def(
+        "make_view_3d_info",
+        [](py::array_t<double> arr) {
+            auto view = make_view_3d(arr);
+            return py::make_tuple(view.shape(), view.strides());
+        },
+        py::arg("arr").noconvert());
+
+    testing.def(
+        "make_const_view_1d_info",
+        [](py::array_t<double> arr) {
+            auto view = make_const_view_1d(arr);
+            return py::make_tuple(view.shape(), view.strides());
+        },
+        py::arg("arr").noconvert());
+    testing.def(
+        "make_const_view_2d_info",
+        [](py::array_t<double> arr) {
+            auto view = make_const_view_2d(arr);
+            return py::make_tuple(view.shape(), view.strides());
+        },
+        py::arg("arr").noconvert());
+    testing.def(
+        "make_const_view_3d_info",
+        [](py::array_t<double> arr) {
+            auto view = make_const_view_3d(arr);
+            return py::make_tuple(view.shape(), view.strides());
+        },
+        py::arg("arr").noconvert());
+
+    // With c_style in the signature pybind11 hands us a contiguous temporary
+    // so make_view accepts any layout, but writes don't reach the caller
+    testing.def(
+        "view_roundtrip_2d_c_style",
+        [](py::array_t<double, py::array::c_style | py::array::forcecast> arr) {
+            auto view = make_view_2d(arr);
+            return return_image_data(new NDArray<double, 2>(view));
+        },
+        py::arg("arr"));
+    testing.def(
+        "view_fill_2d_c_style",
+        [](py::array_t<double, py::array::c_style | py::array::forcecast> arr,
+           double value) {
+            auto view = make_view_2d(arr);
+            view = value;
+        },
+        py::arg("arr"), py::arg("value"));
+}

python/src/module.cpp

@@ -22,6 +22,7 @@
 #include "bind_PixelMap.hpp"
 #include "bind_RawFile.hpp"
 #include "bind_calibration.hpp"
+#include "bind_testing.hpp"
 
 // TODO! migrate the other names
 #include "ctb_raw_file.hpp"
@@ -88,6 +89,7 @@ PYBIND11_MODULE(_aare, m) {
     define_jungfrau_data_file_io_bindings(m);
 
     bind_calibration(m);
+    define_testing_bindings(m);
 
     DEFINE_CLUSTER_BINDINGS(int, 3, 3, uint16_t, i);
     DEFINE_CLUSTER_BINDINGS(double, 3, 3, uint16_t, d);

python/src/np_helper.hpp

@@ -41,31 +41,60 @@ template <typename T> py::array return_vector(std::vector<T> *vec) {
                           free_when_done); // numpy array references this parent
 }
 
-// todo rewrite generic
-template <class T, int Flags>
-auto get_shape_3d(const py::array_t<T, Flags> &arr) {
-    return aare::Shape<3>{arr.shape(0), arr.shape(1), arr.shape(2)};
+// Create a NDView of a numpy array. NDView assumes C-order strides so we
+// reject anything that is not Ndim dimensional and C-contiguous instead of
+// silently reading the data in the wrong order.
+template <ssize_t Ndim, class T, int Flags>
+auto get_checked_shape(const py::array_t<T, Flags> &arr) {
+    if (arr.ndim() != Ndim) {
+        throw py::value_error(
+            fmt::format("Expected {}D array, got {}D", Ndim, arr.ndim()));
+    }
+    if (!(arr.flags() & py::array::c_style)) {
+        throw py::value_error(
+            "Array is not C-contiguous. Use np.ascontiguousarray(arr)");
+    }
+    aare::Shape<Ndim> shape{};
+    for (ssize_t i = 0; i < Ndim; ++i) {
+        shape[i] = arr.shape(i);
+    }
+    return shape;
+}
+
+template <ssize_t Ndim, class T, int Flags>
+auto make_view(py::array_t<T, Flags> &arr) {
+    auto shape = get_checked_shape<Ndim>(arr);
+    return aare::NDView<T, Ndim>(arr.mutable_data(), shape);
+}
+
+// Read only view, also works for numpy arrays that are not writeable
+template <ssize_t Ndim, class T, int Flags>
+auto make_const_view(const py::array_t<T, Flags> &arr) {
+    auto shape = get_checked_shape<Ndim>(arr);
+    return aare::NDView<const T, Ndim>(arr.data(), shape);
 }
 
 template <class T, int Flags> auto make_view_3d(py::array_t<T, Flags> &arr) {
-    return aare::NDView<T, 3>(arr.mutable_data(), get_shape_3d<T, Flags>(arr));
+    return make_view<3>(arr);
 }
-
-template <class T, int Flags>
-auto get_shape_2d(const py::array_t<T, Flags> &arr) {
-    return aare::Shape<2>{arr.shape(0), arr.shape(1)};
-}
-
-template <class T, int Flags>
-auto get_shape_1d(const py::array_t<T, Flags> &arr) {
-    return aare::Shape<1>{arr.shape(0)};
-}
-
 template <class T, int Flags> auto make_view_2d(py::array_t<T, Flags> &arr) {
-    return aare::NDView<T, 2>(arr.mutable_data(), get_shape_2d<T, Flags>(arr));
+    return make_view<2>(arr);
 }
 template <class T, int Flags> auto make_view_1d(py::array_t<T, Flags> &arr) {
-    return aare::NDView<T, 1>(arr.mutable_data(), get_shape_1d<T, Flags>(arr));
+    return make_view<1>(arr);
+}
+
+template <class T, int Flags>
+auto make_const_view_3d(const py::array_t<T, Flags> &arr) {
+    return make_const_view<3>(arr);
+}
+template <class T, int Flags>
+auto make_const_view_2d(const py::array_t<T, Flags> &arr) {
+    return make_const_view<2>(arr);
+}
+template <class T, int Flags>
+auto make_const_view_1d(const py::array_t<T, Flags> &arr) {
+    return make_const_view<1>(arr);
 }
 
 template <typename ClusterType> struct fmt_format_trait; // forward declaration

python/tests/test_ClusterVector.py

@@ -157,5 +157,5 @@ def test_masking_requires_one_dimension():
 
     mask = np.array([[True, False]], dtype=bool)
 
-    with pytest.raises(ValueError, match="one-dimensional"):
+    with pytest.raises(ValueError):
         cv(mask)

python/tests/test_make_view.py

@@ -0,0 +1,226 @@
+# SPDX-License-Identifier: MPL-2.0
+"""Probe make_view in np_helper.hpp through the bindings in _aare.testing"""
+import numpy as np
+import pytest
+
+from aare._aare import testing
+
+SHAPES = {1: (7,), 2: (5, 6), 3: (3, 5, 6)}
+DTYPES = [np.float64, np.uint16]
+
+
+def _probe(name, ndim):
+    return getattr(testing, f"{name}_{ndim}d")
+
+
+def _array(ndim, dtype):
+    shape = SHAPES[ndim]
+    return np.arange(np.prod(shape), dtype=dtype).reshape(shape)
+
+
+def _c_strides(shape):
+    # strides in elements, as used by NDView
+    return tuple(int(np.prod(shape[i + 1 :])) for i in range(len(shape)))
+
+
+def _noncontiguous(ndim, dtype, layout):
+    arr = _array(ndim, dtype)
+    if layout == "sliced":
+        out = np.repeat(arr, 2, axis=-1)[..., ::2]
+    elif layout == "reversed":
+        out = np.ascontiguousarray(arr[..., ::-1])[..., ::-1]
+    elif layout == "fortran":
+        out = np.asfortranarray(arr)
+    elif layout == "transposed":
+        out = np.ascontiguousarray(arr.T).T
+    np.testing.assert_array_equal(out, arr)
+    assert not out.flags.c_contiguous
+    return out
+
+
+VIEWS = ["view", "const_view"]
+
+
+@pytest.mark.parametrize("view", VIEWS)
+@pytest.mark.parametrize("dtype", DTYPES)
+@pytest.mark.parametrize("ndim", [1, 2, 3])
+def test_view_of_contiguous_array(ndim, dtype, view):
+    arr = _array(ndim, dtype)
+
+    out = _probe(f"{view}_roundtrip", ndim)(arr)
+    assert out.dtype == dtype
+    np.testing.assert_array_equal(out, arr)
+
+    shape, strides, ptr = _probe(f"{view}_info", ndim)(arr)
+    assert tuple(shape) == arr.shape
+    assert tuple(strides) == _c_strides(arr.shape)
+    assert ptr == arr.ctypes.data  # no copy
+
+
+@pytest.mark.parametrize("dtype", DTYPES)
+@pytest.mark.parametrize("ndim", [1, 2, 3])
+def test_write_through_view_reaches_caller(ndim, dtype):
+    arr = _array(ndim, dtype)
+    _probe("view_fill", ndim)(arr, 7)
+    assert (arr == 7).all()
+
+
+@pytest.mark.parametrize("dtype", DTYPES)
+@pytest.mark.parametrize(
+    "ndim, layout",
+    [
+        (1, "sliced"),
+        (1, "reversed"),
+        (2, "sliced"),
+        (2, "reversed"),
+        (2, "fortran"),
+        (2, "transposed"),
+        (3, "sliced"),
+        (3, "fortran"),
+        (3, "transposed"),
+    ],
+)
+def test_noncontiguous_array_is_rejected(ndim, layout, dtype):
+    arr = _noncontiguous(ndim, dtype, layout)
+    for name in [
+        "view_roundtrip",
+        "view_info",
+        "const_view_roundtrip",
+        "const_view_info",
+    ]:
+        with pytest.raises(ValueError, match="not C-contiguous"):
+            _probe(name, ndim)(arr)
+
+    # and nothing is written to it
+    before = arr.copy()
+    with pytest.raises(ValueError, match="not C-contiguous"):
+        _probe("view_fill", ndim)(arr, 7)
+    np.testing.assert_array_equal(arr, before)
+
+    # making it contiguous is the way out
+    out = _probe("view_roundtrip", ndim)(np.ascontiguousarray(arr))
+    np.testing.assert_array_equal(out, arr)
+
+
+@pytest.mark.parametrize("view", VIEWS)
+@pytest.mark.parametrize("dtype", DTYPES)
+@pytest.mark.parametrize(
+    "expected, got", [(1, 2), (1, 3), (2, 1), (2, 3), (3, 1), (3, 2)]
+)
+def test_wrong_number_of_dimensions_is_rejected(expected, got, dtype, view):
+    arr = _array(got, dtype)
+    with pytest.raises(ValueError, match=f"Expected {expected}D array, got {got}D"):
+        _probe(f"{view}_roundtrip", expected)(arr)
+
+
+def test_zero_dimensional_array_is_rejected():
+    with pytest.raises(ValueError, match="Expected 1D array, got 0D"):
+        testing.view_roundtrip_1d(np.array(1.0))
+
+
+def test_dimensions_are_checked_before_layout():
+    arr = _noncontiguous(3, np.float64, "fortran")
+    with pytest.raises(ValueError, match="Expected 2D array, got 3D"):
+        testing.view_roundtrip_2d(arr)
+
+
+@pytest.mark.parametrize("view", VIEWS)
+@pytest.mark.parametrize(
+    "shape", [(0,), (1,), (0, 4), (4, 0), (1, 4), (4, 1), (1, 1), (1, 4, 1), (2, 0, 3)]
+)
+def test_empty_and_unit_length_dimensions(shape, view):
+    arr = np.arange(np.prod(shape), dtype=np.float64).reshape(shape)
+    out = _probe(f"{view}_roundtrip", len(shape))(arr)
+    assert out.shape == shape
+    np.testing.assert_array_equal(out, arr)
+
+
+def test_unit_length_dimensions_that_numpy_flags_as_contiguous():
+    # Fortran order and slicing a length one axis still give a C-contiguous array
+    for arr in [
+        np.asfortranarray(np.arange(4.0).reshape(1, 4)),
+        np.asfortranarray(np.arange(4.0).reshape(4, 1)),
+        np.arange(8.0).reshape(2, 4)[:1],
+    ]:
+        assert arr.flags.c_contiguous
+        np.testing.assert_array_equal(testing.view_roundtrip_2d(arr), arr)
+
+
+def test_contiguous_slices_are_viewed_in_place():
+    base = np.arange(60, dtype=np.float64).reshape(10, 6)
+    rows = base[2:5]
+    assert rows.flags.c_contiguous and not rows.flags.owndata
+
+    shape, _, ptr = testing.view_info_2d(rows)
+    assert tuple(shape) == (3, 6)
+    assert ptr == rows.ctypes.data
+    np.testing.assert_array_equal(testing.view_roundtrip_2d(rows), rows)
+
+    testing.view_fill_2d(rows, -1)
+    assert (base[2:5] == -1).all()
+    assert (base[:2] != -1).all() and (base[5:] != -1).all()
+
+
+def test_read_only_array_is_rejected_by_mutable_view():
+    arr = _array(2, np.float64)
+    arr.flags.writeable = False
+    with pytest.raises(ValueError, match="not writeable"):
+        testing.view_roundtrip_2d(arr)
+
+
+@pytest.mark.parametrize("dtype", DTYPES)
+@pytest.mark.parametrize("ndim", [1, 2, 3])
+def test_read_only_array_is_accepted_by_const_view(ndim, dtype):
+    arr = _array(ndim, dtype)
+    arr.flags.writeable = False
+
+    out = _probe("const_view_roundtrip", ndim)(arr)
+    np.testing.assert_array_equal(out, arr)
+
+    _, _, ptr = _probe("const_view_info", ndim)(arr)
+    assert ptr == arr.ctypes.data  # no copy
+
+
+@pytest.mark.parametrize("view", VIEWS)
+@pytest.mark.parametrize("ndim", [1, 2, 3])
+def test_make_view_nd_wrappers(ndim, view):
+    wrapper = getattr(testing, f"make_{view}_{ndim}d_info")
+    arr = _array(ndim, np.float64)
+    shape, strides = wrapper(arr)
+    assert tuple(shape) == arr.shape
+    assert tuple(strides) == _c_strides(arr.shape)
+
+    wrong = _array(ndim % 3 + 1, np.float64)
+    with pytest.raises(ValueError, match=f"Expected {ndim}D array, got {wrong.ndim}D"):
+        wrapper(wrong)
+
+    with pytest.raises(ValueError, match="not C-contiguous"):
+        wrapper(_noncontiguous(ndim, np.float64, "sliced"))
+
+
+def test_probes_dispatch_on_dtype_without_converting():
+    with pytest.raises(TypeError):
+        testing.view_roundtrip_2d(_array(2, np.float32))
+    with pytest.raises(TypeError):
+        testing.view_roundtrip_2d(_array(2, np.float64).tolist())
+
+
+@pytest.mark.parametrize("layout", ["sliced", "reversed", "fortran", "transposed"])
+def test_c_style_signature_converts_before_make_view(layout):
+    # When the binding asks pybind11 for c_style the check in make_view never
+    # fires, the view is of a temporary and writes are lost
+    arr = _noncontiguous(2, np.float64, layout)
+    np.testing.assert_array_equal(testing.view_roundtrip_2d_c_style(arr), arr)
+
+    before = arr.copy()
+    testing.view_fill_2d_c_style(arr, 7)
+    np.testing.assert_array_equal(arr, before)
+
+    contiguous = _array(2, np.float64)
+    testing.view_fill_2d_c_style(contiguous, 7)
+    assert (contiguous == 7).all()
+
+
+def test_c_style_signature_still_checks_dimensions():
+    with pytest.raises(ValueError, match="Expected 2D array, got 3D"):
+        testing.view_roundtrip_2d_c_style(_array(3, np.float64))