diff --git a/python/src/bind_ClusterVector.hpp b/python/src/bind_ClusterVector.hpp index 095c3f8a..32a1fc6c 100644 --- a/python/src/bind_ClusterVector.hpp +++ b/python/src/bind_ClusterVector.hpp @@ -48,10 +48,7 @@ void define_ClusterVector(py::module &m, const std::string &typestr) { "__call__", [](ClusterVector &self, py::array_t 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. diff --git a/python/src/bind_testing.hpp b/python/src/bind_testing.hpp new file mode 100644 index 00000000..a7904c65 --- /dev/null +++ b/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 +#include +#include +#include +#include +#include + +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 +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 arr) { + auto view = make_view(arr); + return return_image_data(new NDArray(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 arr) { + auto view = make_view(arr); + return py::make_tuple(view.shape(), view.strides(), + reinterpret_cast(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 arr) { + auto view = make_const_view(arr); + static_assert( + std::is_same_v>); + auto out = new NDArray(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 arr) { + auto view = make_const_view(arr); + return py::make_tuple(view.shape(), view.strides(), + reinterpret_cast(view.data())); + }, + py::arg("arr").noconvert()); + + // Write through the view + m.def(("view_fill_" + suffix).c_str(), + [](py::array_t arr, T value) { + auto view = make_view(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(testing, "1d"); + define_make_view_probes(testing, "2d"); + define_make_view_probes(testing, "3d"); + define_make_view_probes(testing, "1d"); + define_make_view_probes(testing, "2d"); + define_make_view_probes(testing, "3d"); + + // The wrappers used by the bindings + testing.def( + "make_view_1d_info", + [](py::array_t 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 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 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 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 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 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 arr) { + auto view = make_view_2d(arr); + return return_image_data(new NDArray(view)); + }, + py::arg("arr")); + testing.def( + "view_fill_2d_c_style", + [](py::array_t arr, + double value) { + auto view = make_view_2d(arr); + view = value; + }, + py::arg("arr"), py::arg("value")); +} diff --git a/python/src/module.cpp b/python/src/module.cpp index 8e39d064..b83deaca 100644 --- a/python/src/module.cpp +++ b/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); diff --git a/python/src/np_helper.hpp b/python/src/np_helper.hpp index c0f5c3fd..d21aceff 100644 --- a/python/src/np_helper.hpp +++ b/python/src/np_helper.hpp @@ -41,31 +41,60 @@ template py::array return_vector(std::vector *vec) { free_when_done); // numpy array references this parent } -// todo rewrite generic -template -auto get_shape_3d(const py::array_t &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 +auto get_checked_shape(const py::array_t &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 shape{}; + for (ssize_t i = 0; i < Ndim; ++i) { + shape[i] = arr.shape(i); + } + return shape; +} + +template +auto make_view(py::array_t &arr) { + auto shape = get_checked_shape(arr); + return aare::NDView(arr.mutable_data(), shape); +} + +// Read only view, also works for numpy arrays that are not writeable +template +auto make_const_view(const py::array_t &arr) { + auto shape = get_checked_shape(arr); + return aare::NDView(arr.data(), shape); } template auto make_view_3d(py::array_t &arr) { - return aare::NDView(arr.mutable_data(), get_shape_3d(arr)); + return make_view<3>(arr); } - -template -auto get_shape_2d(const py::array_t &arr) { - return aare::Shape<2>{arr.shape(0), arr.shape(1)}; -} - -template -auto get_shape_1d(const py::array_t &arr) { - return aare::Shape<1>{arr.shape(0)}; -} - template auto make_view_2d(py::array_t &arr) { - return aare::NDView(arr.mutable_data(), get_shape_2d(arr)); + return make_view<2>(arr); } template auto make_view_1d(py::array_t &arr) { - return aare::NDView(arr.mutable_data(), get_shape_1d(arr)); + return make_view<1>(arr); +} + +template +auto make_const_view_3d(const py::array_t &arr) { + return make_const_view<3>(arr); +} +template +auto make_const_view_2d(const py::array_t &arr) { + return make_const_view<2>(arr); +} +template +auto make_const_view_1d(const py::array_t &arr) { + return make_const_view<1>(arr); } template struct fmt_format_trait; // forward declaration diff --git a/python/tests/test_ClusterVector.py b/python/tests/test_ClusterVector.py index 8b077dc0..dd75366e 100644 --- a/python/tests/test_ClusterVector.py +++ b/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) diff --git a/python/tests/test_make_view.py b/python/tests/test_make_view.py new file mode 100644 index 00000000..3c70ddf3 --- /dev/null +++ b/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))