diff --git a/binds/python/bind_immutable.cpp b/binds/python/bind_immutable.cpp index 307367435..bfd9f440d 100644 --- a/binds/python/bind_immutable.cpp +++ b/binds/python/bind_immutable.cpp @@ -80,9 +80,8 @@ void bind_immutable_module(py::module& m) { // Property accessors .def_property_readonly( "points", - [](morphio::Morphology* morpho) { - const auto& data = morpho->points(); - return py::array(static_cast(data.size()), data.data()); + [](const morphio::Morphology& morpho) { + return internal_vector_as_readonly_array(morpho.points(), morpho); }, "Returns a list with all points from all sections (soma points are not included)\n" "Note: points belonging to the n'th section are located at indices:\n" @@ -96,17 +95,15 @@ void bind_immutable_module(py::module& m) { .def_property_readonly( "diameters", [](const morphio::Morphology& morpho) { - const auto& data = morpho.diameters(); - return py::array(static_cast(data.size()), data.data()); + return internal_vector_as_readonly_array(morpho.diameters(), morpho); }, "Returns a list with all diameters from all sections (soma points are not included)\n" "Note: diameters belonging to the n'th section are located at indices:\n" "[Morphology.sectionOffsets(n), Morphology.sectionOffsets(n+1)[") .def_property_readonly( "perimeters", - [](const morphio::Morphology& obj) { - const auto& data = obj.perimeters(); - return py::array(static_cast(data.size()), data.data()); + [](const morphio::Morphology& morpho) { + return internal_vector_as_readonly_array(morpho.perimeters(), morpho); }, "Returns a list with all perimeters from all sections (soma points are not included)\n" "Note: perimeters belonging to the n'th section are located at indices:\n" @@ -125,8 +122,7 @@ void bind_immutable_module(py::module& m) { .def_property_readonly( "section_types", [](const morphio::Morphology& morph) { - const auto& data = morph.sectionTypes(); - return py::array(static_cast(data.size()), data.data()); + return internal_vector_as_readonly_array(morph.sectionTypes(), morph); }, "Returns a vector with the section type of every section") .def_property_readonly("connectivity", @@ -268,7 +264,9 @@ void bind_immutable_module(py::module& m) { "(dendrite, axon, ...)") .def_property_readonly( "points", - [](morphio::Section* section) { return span_array_to_ndarray(section->points()); }, + [](const morphio::Section& self) { + return internal_vector_as_readonly_array(self.points(), self); + }, "Returns list of section's point coordinates") .def_property_readonly( @@ -278,11 +276,15 @@ void bind_immutable_module(py::module& m) { .def_property_readonly( "diameters", - [](morphio::Section* section) { return span_to_ndarray(section->diameters()); }, + [](const morphio::Section& self) { + return internal_vector_as_readonly_array(self.diameters(), self); + }, "Returns list of section's point diameters") .def_property_readonly( "perimeters", - [](morphio::Section* section) { return span_to_ndarray(section->perimeters()); }, + [](const morphio::Section& self) { + return internal_vector_as_readonly_array(self.perimeters(), self); + }, "Returns list of section's point perimeters") .def("is_heterogeneous", @@ -411,18 +413,16 @@ void bind_immutable_module(py::module& m) { // Property accessors .def_property_readonly( "points", - [](morphio::DendriticSpine* morpho) { - const auto& data = morpho->points(); - return py::array(static_cast(data.size()), data.data()); + [](const morphio::DendriticSpine& morph) { + return internal_vector_as_readonly_array(morph.points(), morph); }, "Returns a list with all points from all sections\n" "Note: points belonging to the n'th section are located at indices:\n" "[DendriticSpine.sectionOffsets(n), DendriticSpine.sectionOffsets(n+1)[") .def_property_readonly( "diameters", - [](const morphio::DendriticSpine& morpho) { - const auto& data = morpho.diameters(); - return py::array(static_cast(data.size()), data.data()); + [](const morphio::DendriticSpine& morph) { + return internal_vector_as_readonly_array(morph.diameters(), morph); }, "Returns a list with all diameters from all sections\n" "Note: diameters belonging to the n'th section are located at indices:\n" @@ -443,8 +443,7 @@ void bind_immutable_module(py::module& m) { .def_property_readonly( "section_types", [](const morphio::DendriticSpine& morph) { - const auto& data = morph.sectionTypes(); - return py::array(static_cast(data.size()), data.data()); + return internal_vector_as_readonly_array(morph.diameters(), morph); }, "Returns a vector with the section type of every section") .def_property_readonly("connectivity", diff --git a/binds/python/bind_vasculature.cpp b/binds/python/bind_vasculature.cpp index 82e394b29..7f2c45358 100644 --- a/binds/python/bind_vasculature.cpp +++ b/binds/python/bind_vasculature.cpp @@ -54,9 +54,8 @@ void bind_vasculature(py::module& m) { // Property accessors .def_property_readonly( "points", - [](morphio::vasculature::Vasculature* morpho) { - return py::array(static_cast(morpho->points().size()), - morpho->points().data()); + [](const morphio::vasculature::Vasculature& self) { + return internal_vector_as_readonly_array(self.points(), self); }, "Returns a list with all points from all sections") @@ -67,24 +66,21 @@ void bind_vasculature(py::module& m) { .def_property_readonly( "diameters", - [](morphio::vasculature::Vasculature* morpho) { - auto diameters = morpho->diameters(); - return py::array(static_cast(diameters.size()), diameters.data()); + [](const morphio::vasculature::Vasculature& self) { + return internal_vector_as_readonly_array(self.diameters(), self); }, "Returns a list with all diameters from all sections") .def_property_readonly( "section_types", - [](morphio::vasculature::Vasculature* obj) { - auto data = obj->sectionTypes(); - return py::array(static_cast(data.size()), data.data()); + [](const morphio::vasculature::Vasculature& self) { + return internal_vector_as_readonly_array(self.sectionTypes(), self); }, "Returns a vector with the section type of every section") .def_property_readonly( "section_connectivity", - [](morphio::vasculature::Vasculature* morpho) { - return py::array(static_cast(morpho->sectionConnectivity().size()), - morpho->sectionConnectivity().data()); + [](const morphio::vasculature::Vasculature& self) { + return internal_vector_as_readonly_array(self.sectionConnectivity(), self); }, "Returns a 2D array of the section connectivity") @@ -128,8 +124,8 @@ void bind_vasculature(py::module& m) { "Returns the morphological type of this section") .def_property_readonly( "points", - [](morphio::vasculature::Section* section) { - return span_array_to_ndarray(section->points()); + [](const morphio::vasculature::Section& self) { + return internal_vector_as_readonly_array(self.points(), self); }, "Returns list of section's point coordinates") @@ -140,8 +136,8 @@ void bind_vasculature(py::module& m) { .def_property_readonly( "diameters", - [](morphio::vasculature::Section* section) { - return span_to_ndarray(section->diameters()); + [](const morphio::vasculature::Section& self) { + return internal_vector_as_readonly_array(self.points(), self); }, "Returns list of section's point diameters") diff --git a/binds/python/bindings_utils.h b/binds/python/bindings_utils.h index 4d6ff2ebb..9816cbf68 100644 --- a/binds/python/bindings_utils.h +++ b/binds/python/bindings_utils.h @@ -48,3 +48,21 @@ inline py::array_t as_pyarray(Sequence&& seq) { capsule // numpy array references this parent ); } + +/** + * @brief Exposes an internal vector (not temporary) as a non-writeable numpy array. + */ +template +inline py::array internal_vector_as_readonly_array(const Sequence& seq, const Parent& parent) { + // The correct base must be used here so that refcounting is done correctly on the python + // side. The parent that holds that returned data should be passed as a base. See: + // https://github.com/pybind/pybind11/issues/2271#issuecomment-650740842 + auto res = py::array(static_cast(seq.size()), seq.data(), py::cast(parent)); + + // pybind11 casts away const-ness in return values. Thus, we need to explicitly set the + // numpy array flag to not writeable to prevent mutating the underlying data. See limitations: + // https://pybind11.readthedocs.io/en/stable/limitations.html + py::detail::array_proxy(res.ptr())->flags &= ~py::detail::npy_api::NPY_ARRAY_WRITEABLE_; + + return res; +} diff --git a/tests/test_4_immut.py b/tests/test_4_immut.py index 8b1a31e3f..f5fe3d5fc 100644 --- a/tests/test_4_immut.py +++ b/tests/test_4_immut.py @@ -46,6 +46,7 @@ def test_components(): for cell in CELLS.values(): assert cell.n_points == len(cell.points) assert cell.section(0).n_points == len(cell.section(0).points) + assert not cell.points.flags.writeable def test_is_root(): diff --git a/tests/test_readonly_views.py b/tests/test_readonly_views.py new file mode 100644 index 000000000..c3705cf15 --- /dev/null +++ b/tests/test_readonly_views.py @@ -0,0 +1,223 @@ +""" +Readonly Morphology returns readonly numpy views (not copies) to its internal +points, diameters and perimeters datasets. Here it is ensured that refcounting +of the Morphology is increased when these views are available on the python side. +Morphology can be released only when all references to itself and its views are +freed. The standard datasets of all the readonly morphologies are: + - points + - diameters + - section_types +""" +import sys +import pytest +import morphio +from pathlib import Path + +import numpy as np +from numpy import testing as npt + +DATA_DIR = Path(__file__).parent.resolve() / "data" + +READONLY_CLASSES_PATHS = [ + (morphio.Morphology, DATA_DIR / "simple-heterogeneous-neurite.swc"), + (morphio.GlialCell, DATA_DIR / "astrocyte.h5"), + (morphio.vasculature.Vasculature, DATA_DIR / "h5/vasculature1.h5"), + (morphio.DendriticSpine, DATA_DIR / "h5/v1/simple-dendritric-spine.h5") +] + + +def _assert_refcount(obj, expected_refcount): + """Check that refcount of the obj has the expected value. An object + passed in _assert_refcount has 3 extra refcounts: + - From the obj variable + - From passing in _assert_refcount + - From passing in sys.getrefcount + """ + actual_refcount = sys.getrefcount(obj) - 3 + assert actual_refcount == expected_refcount, ( + f"\nActual: {actual_refcount}\n" + f"Expected: {expected_refcount}" + ) + + +@pytest.mark.parametrize("MorphologyClass, path", READONLY_CLASSES_PATHS) +def test_morphology_level_immutability(MorphologyClass, path): + + morph = MorphologyClass(path) + + # ensure we cannot overwrite the attributes + with pytest.raises(AttributeError): + morph.points = np.ones(morph.points.shape) + + with pytest.raises(AttributeError): + morph.diameters = np.ones(morph.diameters.shape) + + with pytest.raises(AttributeError): + morph.section_types = np.ones(morph.section_types.shape) + + # ensure we cannot change the underlying data + with pytest.raises(ValueError): + morph.points[:] = np.ones(morph.points.shape) + + with pytest.raises(ValueError): + morph.diameters[:] = np.ones(morph.diameters.shape) + + with pytest.raises(ValueError): + morph.section_types[:] = np.ones(morph.section_types.shape) + + # check if flag set correctly + assert not morph.points.flags["WRITEABLE"] + assert not morph.diameters.flags["WRITEABLE"] + assert not morph.section_types.flags["WRITEABLE"] + + # check that the parent of the array is the morphology + assert morph.points.base is morph + assert morph.diameters.base is morph + assert morph.section_types.base is morph + + # make the same tests when the datasets are assigned to a variable + points = morph.points + _assert_refcount(points, 1) + assert points.base is morph + assert points.flags["WRITEABLE"] == False + + diameters = morph.diameters + _assert_refcount(diameters, 1) + assert diameters.base is morph + assert diameters.flags["WRITEABLE"] == False + + section_types = morph.section_types + _assert_refcount(section_types, 1) + assert section_types.base is morph + assert section_types.flags["WRITEABLE"] == False + + +@pytest.mark.parametrize("MorphologyClass, path", READONLY_CLASSES_PATHS) +def test_section_level_immutability__iter(MorphologyClass, path): + + morph = MorphologyClass(path) + _assert_refcount(morph, 1) + + for section in morph.iter(): + + _assert_refcount(section, 1) + # 1 (morph) + 1 (iter()) + _assert_refcount(morph, 2) + + points = section.points + assert points.base is section + _assert_refcount(section, 2) + _assert_refcount(morph, 2) + + diameters = section.diameters + assert diameters.base is section + _assert_refcount(section, 3) + _assert_refcount(morph, 2) + + #del points, diameters + #_assert_refcount(section, 1) + + +@pytest.mark.parametrize("MorphologyClass, path", READONLY_CLASSES_PATHS) +def test_section_level_immutability__sections(MorphologyClass, path): + + morph = MorphologyClass(path) + _assert_refcount(morph, 1) + + for section in morph.sections: + + # 1 (section) and 1 (morph.sections container) + _assert_refcount(section, 2) + _assert_refcount(morph, 1) + + points = section.points + assert points.base is section + _assert_refcount(section, 3) + _assert_refcount(morph, 1) + + diameters = section.diameters + assert diameters.base is section + _assert_refcount(section, 4) + _assert_refcount(morph, 1) + + del points, diameters + _assert_refcount(section, 2) + +@pytest.mark.parametrize("MorphologyClass, path", READONLY_CLASSES_PATHS) +def test_copies_mutability(MorphologyClass, path): + + morph = MorphologyClass(path) + + _assert_refcount(morph, 1) + + points_copy = morph.points.copy() + diameters_copy = morph.diameters.copy() + section_types_copy = morph.section_types.copy() + + # we copy the points, morph refcount should not change + _assert_refcount(morph, 1) + + # sanity check + npt.assert_array_almost_equal(points_copy, morph.points) + npt.assert_array_almost_equal(diameters_copy, morph.diameters) + npt.assert_array_almost_equal(section_types_copy, morph.section_types) + + # check that new datasets are writeable + new_points = morph.points + np.zeros(morph.points.shape) + assert new_points.flags["WRITEABLE"] + + new_diameters = morph.diameters + np.zeros(morph.diameters.shape) + assert new_diameters.flags["WRITEABLE"] + + new_section_types = morph.section_types + np.zeros(morph.section_types.shape) + assert new_section_types.flags["WRITEABLE"] + + +@pytest.mark.parametrize("MorphologyClass, path", READONLY_CLASSES_PATHS) +def test_morphology_refcount(MorphologyClass, path): + morph = MorphologyClass(path) + + points = morph.points + diameters = morph.diameters + section_types = morph.section_types + + # morph refcount should be 1 (morph) + 3 (points, diameters, section_types) + _assert_refcount(morph, 4) + + # removing the points, should decrease morph refcount by 1 + del points + _assert_refcount(morph, 3) + + # removing the diameters, should decrease morph refcount by 1 + del diameters + _assert_refcount(morph, 2) + + del section_types + _assert_refcount(morph, 1) + + +@pytest.mark.parametrize("MorphologyClass, path", READONLY_CLASSES_PATHS) +def test_morphology_not_released_while_datasets_in_scope(MorphologyClass, path): + + morph = MorphologyClass(path) + + # morphology should not be released if it's del/outscoped in python + points = morph.points + diameters = morph.diameters + section_types = morph.section_types + + points_copy = points.copy() + diameters_copy = diameters.copy() + section_types_copy = section_types.copy() + + del morph + _assert_refcount(points, 1) + npt.assert_array_almost_equal(points_copy, points) + del points + + _assert_refcount(diameters, 1) + npt.assert_array_almost_equal(diameters_copy, diameters) + del diameters + + _assert_refcount(section_types, 1) + npt.assert_array_almost_equal(section_types_copy, section_types)