diff --git a/README.md b/README.md index e874810..b741c5e 100644 --- a/README.md +++ b/README.md @@ -16,6 +16,31 @@ pip3 install . pytest ./cdt_bindings_test.py ``` +## Usage + +```python +import numpy as np +import PythonCDT as cdt + +t = cdt.Triangulation(cdt.VertexInsertionOrder.AUTO, cdt.IntersectingConstraintEdges.TRY_RESOLVE, 0.0) +t.insert_vertices(np.array([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]])) +t.insert_edges(np.array([[0, 3]], dtype=np.uintc)) +t.erase_super_triangle() + +vertices = t.vertices_array() # numpy array with fields 'x' and 'y' +triangles = t.triangles_array() # numpy array with fields 'vertices' and 'neighbors' +triangles["vertices"] # (T, 3) vertex indices into vertices +``` + +`vertices_array()` and `triangles_array()` return copies. With `copy=False` they return read-only views of the +triangulation's memory without copying; a view is invalidated by any call that modifies the triangulation. + +### Threads + +Triangulating releases the GIL, so separate triangulations can be built in parallel on Python threads. +One triangulation can be shared between threads: calls on it wait for each other. +Iterators (`*_iter()`) and `copy=False` views are not protected: don't use them while another thread modifies the triangulation. + ## License [Mozilla Public License, v. 2.0](https://www.mozilla.org/en-US/MPL/2.0/FAQ/) diff --git a/benchmark_gil_release.py b/benchmark_gil_release.py new file mode 100644 index 0000000..8aa9398 --- /dev/null +++ b/benchmark_gil_release.py @@ -0,0 +1,47 @@ +"""Shows the effect of releasing the GIL during triangulation""" + +import os +import threading +import time +from concurrent.futures import ThreadPoolExecutor + +import numpy as np + +import PythonCDT as cdt + +N_THREADS = min(8, os.cpu_count()) +N_VERTICES = 200_000 + + +def triangulate(vertices) -> None: + t = cdt.Triangulation(cdt.VertexInsertionOrder.AUTO, cdt.IntersectingConstraintEdges.NOT_ALLOWED, 0.0) + t.insert_vertices(vertices) + + +def seconds(f) -> float: + start = time.perf_counter() + f() + return time.perf_counter() - start + + +rng = np.random.default_rng(0) +batches = [rng.random((N_VERTICES, 2)) for _ in range(N_THREADS)] + +# separate triangulations built in parallel threads +sequential = seconds(lambda: [triangulate(b) for b in batches]) +with ThreadPoolExecutor(N_THREADS) as pool: + threaded = seconds(lambda: list(pool.map(triangulate, batches))) +print(f"{N_THREADS} triangulations of {N_VERTICES} vertices each") +print(f" sequential {sequential:8.3f} s") +print(f" {N_THREADS} threads {threaded:8.3f} s ({sequential / threaded:.1f}x speedup)") + +# a Python thread keeps running while another thread triangulates +worker = threading.Thread(target=triangulate, args=(batches[0],)) +ticks = [time.perf_counter()] +worker.start() +while worker.is_alive(): + time.sleep(0.001) + ticks.append(time.perf_counter()) +longest_stall = max(np.diff(ticks)) +print(f"main thread ticking every 1 ms during {ticks[-1] - ticks[0]:.3f} s triangulation") +print(f" ticks {len(ticks) - 1:>6}, longest stall {1e3 * longest_stall:8.1f} ms") diff --git a/benchmark_list_vs_array_access.py b/benchmark_list_vs_array_access.py new file mode 100644 index 0000000..aaad84a --- /dev/null +++ b/benchmark_list_vs_array_access.py @@ -0,0 +1,17 @@ +"""Compares reading a triangulation as Python lists vs numpy arrays""" + +import timeit + +import numpy as np + +import PythonCDT as cdt + +t = cdt.Triangulation(cdt.VertexInsertionOrder.AUTO, cdt.IntersectingConstraintEdges.NOT_ALLOWED, 0.0) +t.insert_vertices(np.random.default_rng(0).random((200_000, 2))) + +print(f"{'':<10}{'count':>9}{'list':>12}{'copy':>12}{'view':>12}") +for name in ["vertices", "triangles"]: + array = getattr(t, f"{name}_array") + ways = [lambda: getattr(t, name), array, lambda: array(copy=False)] + ms = [1e3 * min(timeit.repeat(way, number=1, repeat=5)) for way in ways] + print(f"{name:<10}{len(array()):>9}" + "".join(f"{m:>9.3f} ms" for m in ms)) diff --git a/cdt_bindings.cpp b/cdt_bindings.cpp index 692dc4b..02e598c 100644 --- a/cdt_bindings.cpp +++ b/cdt_bindings.cpp @@ -11,6 +11,7 @@ #include #include +#include #include #include #include @@ -23,11 +24,145 @@ namespace py = pybind11; using coord_t = double; using NearPointLocator_t = CDT::LocatorKDTree; using V2d = CDT::V2d; -using Triangulation = CDT::Triangulation; +using CdtTriangulation = CDT::Triangulation; namespace { +struct Triangulation : CdtTriangulation +{ + using CdtTriangulation::CdtTriangulation; + std::mutex mutex; +}; + +struct LockWithoutGil +{ + explicit LockWithoutGil(Triangulation& t) + : lock(t.mutex) + {} + py::gil_scoped_release release; // constructed first, destroyed last + std::lock_guard lock; +}; + +struct LockWithGil +{ + explicit LockWithGil(Triangulation& t) + : lock(t.mutex, std::try_to_lock) + { + if(!lock.owns_lock()) + { + py::gil_scoped_release release; + lock.lock(); + } + } + std::unique_lock lock; +}; + +// Binds a container member as property `name`, `name_count()` and `name_iter()` +template +void def_container( + py::class_& cls, + const std::string& name, + Container CdtTriangulation::*member) +{ + cls.def_property_readonly( + name.c_str(), + [member](Triangulation& t) { + LockWithGil lock(t); + return py::cast( + t.*member, + py::return_value_policy::reference_internal, + py::cast(&t)); + }, + "Items reference the triangulation's memory; they are invalidated " + "by any call that modifies the triangulation.") + .def( + (name + "_count").c_str(), + [member](Triangulation& t) { + LockWithGil lock(t); + return (t.*member).size(); + }) + .def( + (name + "_iter").c_str(), + [member](Triangulation& t) { + LockWithGil lock(t); + const Container& container = t.*member; + return py::make_iterator(container.begin(), container.end()); + }, + py::keep_alive<0, 1>(), + "Iterates the triangulation's memory; the iterator is invalidated " + "by any call that modifies the triangulation."); +} + +// Binds a vector member as `name_array(*, copy=True)` returning a numpy array +template +void def_array( + py::class_& cls, + const std::string& name, + std::vector CdtTriangulation::*member) +{ + cls.def( + (name + "_array").c_str(), + [member](Triangulation& t, bool copy) { + LockWithGil lock(t); + const std::vector& items = t.*member; + const auto size = static_cast(items.size()); + if (copy) + return py::array_t(size, items.data()); + // the view's base is the Python wrapper of t, keeping t alive + py::array_t view(size, items.data(), py::cast(&t)); + view.attr("setflags")(py::arg("write") = false); + return view; + }, + py::kw_only(), + py::arg("copy") = true, + "Copy as a numpy structured array. copy=False returns a read-only view " + "of the triangulation's memory instead; it is invalidated by any call " + "that modifies the triangulation."); +} + +template +struct BufferPair +{ + T v[2]; +}; + +// Reads a C-contiguous buffer of shape (2N,) or (N, 2) as N pairs +template +std::pair*, std::size_t> buffer_pairs( + const py::buffer_info& info, + const std::string& type_name, + const std::string& item_name) +{ + if (info.format != py::format_descriptor::format()) + { + throw std::runtime_error( + "Incompatible format: expected a " + type_name + " array!"); + } + if (info.ndim != 1 && info.ndim != 2) + { + throw std::runtime_error("Incompatible buffer dimension!"); + } + if (info.ndim == 2 ? info.shape[1] != 2 : info.shape[0] % 2 != 0) + { + throw std::runtime_error( + "Buffer must hold " + type_name + " pairs (2 per " + item_name + + "): shape (N, 2) or (2N,)!"); + } + // strides of dimensions with a single element are irrelevant + py::ssize_t expected_stride = info.itemsize; + for (py::ssize_t d = info.ndim - 1; d >= 0; --d) + { + if (info.shape[d] > 1 && info.strides[d] != expected_stride) + { + throw std::runtime_error( + "Buffer must be C-contiguous: use numpy.ascontiguousarray!"); + } + expected_stride *= info.shape[d]; + } + return {static_cast*>(info.ptr), info.size / 2}; +} + std::string TriInd2str(CDT::TriInd it) { return it != CDT::noNeighbor ? std::to_string(it) : "-"; @@ -172,208 +307,134 @@ PYBIND11_MODULE(PythonCDT, m) return oss.str(); }); - py::class_(m, "Triangulation") - .def( - py::init< - CDT::VertexInsertionOrder::Enum, - CDT::IntersectingConstraintEdges::Enum, - coord_t>(), - py::arg("vertex_insertion_order"), - py::arg("intersecting_edges_strategy"), - py::arg("min_dist_to_constraint_edge")) - // vertices - .def_readonly("vertices", &Triangulation::vertices) - .def( - "vertices_count", - [](const Triangulation& t) { return t.vertices.size(); }) - .def( - "vertices_iter", - [](const Triangulation& t) -> py::iterator { - return py::make_iterator(t.vertices.begin(), t.vertices.end()); - }, - py::keep_alive<0, 1>()) - // triangles - .def_readonly("triangles", &Triangulation::triangles) - .def( - "triangles_count", - [](const Triangulation& t) { return t.triangles.size(); }) - .def( - "triangles_iter", - [](const Triangulation& t) -> py::iterator { - return py::make_iterator( - t.triangles.begin(), t.triangles.end()); - }, - py::keep_alive<0, 1>()) - // fixed edges - .def_readonly("fixed_edges", &Triangulation::fixedEdges) - .def( - "fixed_edges_count", - [](const Triangulation& t) { return t.fixedEdges.size(); }) - .def( - "fixed_edges_iter", - [](const Triangulation& t) -> py::iterator { - return py::make_iterator( - t.fixedEdges.begin(), t.fixedEdges.end()); - }, - py::keep_alive<0, 1>()) - // overlaps - .def_readonly("overlap_count", &Triangulation::overlapCount) - .def( - "overlap_count_count", - [](const Triangulation& t) { return t.overlapCount.size(); }) - .def( - "overlap_count_iter", - [](const Triangulation& t) -> py::iterator { - return py::make_iterator( - t.overlapCount.begin(), t.overlapCount.end()); - }, - py::keep_alive<0, 1>()) - // piece to originals mapping for edges - .def_readonly("piece_to_originals", &Triangulation::pieceToOriginals) - .def( - "piece_to_originals_count", - [](const Triangulation& t) { return t.pieceToOriginals.size(); }) - .def( - "piece_to_originals_iter", - [](const Triangulation& t) -> py::iterator { - return py::make_iterator( - t.pieceToOriginals.begin(), t.pieceToOriginals.end()); - }, - py::keep_alive<0, 1>()) + py::class_ triangulation(m, "Triangulation"); + triangulation.def( + py::init< + CDT::VertexInsertionOrder::Enum, + CDT::IntersectingConstraintEdges::Enum, + coord_t>(), + py::arg("vertex_insertion_order"), + py::arg("intersecting_edges_strategy"), + py::arg("min_dist_to_constraint_edge")); + def_container(triangulation, "vertices", &CdtTriangulation::vertices); + def_container(triangulation, "triangles", &CdtTriangulation::triangles); + def_container(triangulation, "fixed_edges", &CdtTriangulation::fixedEdges); + def_container( + triangulation, "overlap_count", &CdtTriangulation::overlapCount); + def_container( + triangulation, + "piece_to_originals", + &CdtTriangulation::pieceToOriginals); + def_array(triangulation, "vertices", &CdtTriangulation::vertices); + def_array(triangulation, "triangles", &CdtTriangulation::triangles); + triangulation // methods .def( "insert_vertices", - static_cast&)>( - &Triangulation::insertVertices), + [](Triangulation& t, const std::vector& vertices) { + LockWithoutGil lock(t); + t.insertVertices(vertices); + }, py::arg("vertices")) .def( "insert_vertices", [](Triangulation& t, py::buffer b) { const py::buffer_info info = b.request(); - // sanity checks - if (info.format != py::format_descriptor::format()) - { - throw std::runtime_error( - "Incompatible format: expected a double array!"); - } - if (info.ndim != 1 && info.ndim != 2) - { - throw std::runtime_error("Incompatible buffer dimension!"); - } - if (info.size % 2 != 0) - { - throw std::runtime_error("Buffer must hold even number of " - "coordinates (2 per vertex)!"); - } - // create from buffer - struct XY - { - coord_t xy[2]; - }; - const std::size_t n_vert = info.size / 2; - const XY* const ptr = static_cast(info.ptr); + const auto [ptr, n_vert] = + buffer_pairs(info, "double", "vertex"); + LockWithoutGil lock(t); t.insertVertices( ptr, ptr + n_vert, - [](const XY& v) { return v.xy[0]; }, - [](const XY& v) { return v.xy[1]; }); + [](const BufferPair& v) { return v.v[0]; }, + [](const BufferPair& v) { return v.v[1]; }); }, py::arg("vertex_buffer")) .def( "insert_edges", - static_cast&)>( - &Triangulation::insertEdges), + [](Triangulation& t, const std::vector& edges) { + LockWithoutGil lock(t); + t.insertEdges(edges); + }, py::arg("edges")) .def( "insert_edges", [](Triangulation& t, py::buffer b) { const py::buffer_info info = b.request(); - // sanity checks - if (info.format != - py::format_descriptor::format()) - { - throw std::runtime_error( - "Incompatible format: expected a CDT::VertInd array!"); - } - if (info.ndim != 1 && info.ndim != 2) - { - throw std::runtime_error("Incompatible buffer dimension!"); - } - if (info.size % 2 != 0) - { - throw std::runtime_error("Buffer must hold even number of " - "CDT::VertInd (2 per edge)!"); - } - // create from buffer - struct EdgeData - { - CDT::VertInd vv[2]; - }; - const std::size_t n_vert = info.size / 2; - const EdgeData* const ptr = static_cast(info.ptr); + const auto [ptr, n_edges] = + buffer_pairs(info, "CDT::VertInd", "edge"); + LockWithoutGil lock(t); t.insertEdges( ptr, - ptr + n_vert, - [](const EdgeData& e) { return e.vv[0]; }, - [](const EdgeData& e) { return e.vv[1]; }); + ptr + n_edges, + [](const BufferPair& e) { return e.v[0]; }, + [](const BufferPair& e) { return e.v[1]; }); }, py::arg("edge_buffer")) .def( "conform_to_edges", - static_cast&)>( - &Triangulation::conformToEdges), + [](Triangulation& t, const std::vector& edges) { + LockWithoutGil lock(t); + t.conformToEdges(edges); + }, py::arg("edges")) .def( "conform_to_edges", [](Triangulation& t, py::buffer b) { const py::buffer_info info = b.request(); - // sanity checks - if (info.format != - py::format_descriptor::format()) - { - throw std::runtime_error( - "Incompatible format: expected a CDT::VertInd array!"); - } - if (info.ndim != 1 && info.ndim != 2) - { - throw std::runtime_error("Incompatible buffer dimension!"); - } - if (info.size % 2 != 0) - { - throw std::runtime_error("Buffer must hold even number of " - "CDT::VertInd (2 per edge)!"); - } - // create from buffer - struct EdgeData - { - CDT::VertInd vv[2]; - }; - const std::size_t n_vert = info.size / 2; - const EdgeData* const ptr = static_cast(info.ptr); + const auto [ptr, n_edges] = + buffer_pairs(info, "CDT::VertInd", "edge"); + LockWithoutGil lock(t); t.conformToEdges( ptr, - ptr + n_vert, - [](const EdgeData& e) { return e.vv[0]; }, - [](const EdgeData& e) { return e.vv[1]; }); + ptr + n_edges, + [](const BufferPair& e) { return e.v[0]; }, + [](const BufferPair& e) { return e.v[1]; }); }, py::arg("edge_buffer")) - .def("erase_super_triangle", &Triangulation::eraseSuperTriangle) - .def("erase_outer_triangles", &Triangulation::eraseOuterTriangles) + .def( + "erase_super_triangle", + [](Triangulation& t) { + LockWithoutGil lock(t); + t.eraseSuperTriangle(); + }) + .def( + "erase_outer_triangles", + [](Triangulation& t) { + LockWithoutGil lock(t); + t.eraseOuterTriangles(); + }) .def( "erase_outer_triangles_and_holes", - &Triangulation::eraseOuterTrianglesAndHoles) - .def("is_finalized", &Triangulation::isFinalized) + [](Triangulation& t) { + LockWithoutGil lock(t); + t.eraseOuterTrianglesAndHoles(); + }) + .def( + "is_finalized", + [](Triangulation& t) { + LockWithGil lock(t); + return t.isFinalized(); + }) .def( "calculate_triangle_depths", - &Triangulation::calculateTriangleDepths) + [](Triangulation& t) { + LockWithoutGil lock(t); + return t.calculateTriangleDepths(); + }) .def( "remove_triangles", - static_cast( - &Triangulation::removeTriangles), + [](Triangulation& t, const CDT::TriIndUSet& triangle_indices) { + LockWithoutGil lock(t); + t.removeTriangles(triangle_indices); + }, py::arg("triangle_indices")); m.def( "verify_topology", - &CDT::verifyTopology, + [](Triangulation& t) { + LockWithoutGil lock(t); + return CDT::verifyTopology(t); + }, py::arg("triangulation")); } diff --git a/cdt_bindings_test.py b/cdt_bindings_test.py index c98d4e9..0a0954a 100644 --- a/cdt_bindings_test.py +++ b/cdt_bindings_test.py @@ -10,6 +10,7 @@ import pytest import tempfile import hashlib +from concurrent.futures import ThreadPoolExecutor import PythonCDT as cdt @@ -73,6 +74,12 @@ def test_Triangulation() -> None: assert len(t.triangles) == 2, "Wrong triangle count in triangulation" assert len(t.fixed_edges) == 1, "Wrong fixed edge count in triangulation" + # properties give references into the triangulation + vertices = t.vertices + vertices[3].x = 42 + assert t.vertices[3] == cdt.V2d(42, -0.5), "Vertex property must reference triangulation's vertices" + vertices[3].x = 0 + # test retrieving triangulation data using iterators assert t.vertices_count() == len(t.vertices), "Wrong vertex count" assert t.triangles_count() == len(t.triangles), "Wrong triangle count" @@ -190,3 +197,63 @@ def test_insert_conform_edges(ee) -> None: assert len(t.triangles) == 19, "Wrong triangle count in triangulation" assert len(t.fixed_edges) == 6, "Wrong fixed edge count in triangulation" assert triangulation_md5_checksum(t) == '9c87b435e247c1658ec0f04af3340dc7', "Wrong OFF file contents" + + +@pytest.mark.parametrize("copy", [True, False]) +def test_arrays(copy) -> None: + t = cdt.Triangulation(cdt.VertexInsertionOrder.AS_PROVIDED, cdt.IntersectingConstraintEdges.NOT_ALLOWED, 0.0) + assert t.vertices_array(copy=copy).shape == t.triangles_array(copy=copy).shape == (0,), \ + "Empty triangulation must give empty arrays" + + t.insert_vertices([cdt.V2d(-1, 0), cdt.V2d(0, 0.5), cdt.V2d(1, 0), cdt.V2d(0, -0.5)]) + vertices = t.vertices_array(copy=copy) + triangles = t.triangles_array(copy=copy) + assert vertices.tolist() == [(v.x, v.y) for v in t.vertices] + assert triangles["vertices"].tolist() == [list(tri.vertices) for tri in t.triangles] + assert triangles["neighbors"].tolist() == [list(tri.neighbors) for tri in t.triangles] + + expected = [vertices.copy(), triangles.copy()] + del t + for arr, before in zip([vertices, triangles], expected): + assert arr.flags.owndata == arr.flags.writeable == copy, "Copy must be owned and writeable, view neither" + assert np.array_equal(arr, before), "Array must stay valid after the triangulation is deleted" + + +def test_shared_triangulation() -> None: + """Threads can insert into and read from one triangulation concurrently""" + vertices = np.random.default_rng(0).random((160_000, 2)) + t = cdt.Triangulation(cdt.VertexInsertionOrder.AUTO, cdt.IntersectingConstraintEdges.NOT_ALLOWED, 0.0) + + def insert_and_read(batch) -> None: + t.insert_vertices(batch) + assert t.triangles_array()["vertices"].max() < t.vertices_count() + + with ThreadPoolExecutor(max_workers=8) as pool: + list(pool.map(insert_and_read, np.array_split(vertices, 8))) + assert t.vertices_count() == len(vertices) + 3, "All vertices plus the super-triangle must be inserted" + assert cdt.verify_topology(t) + + +def test_insert_buffers_must_be_contiguous_pairs() -> None: + """Buffer overloads reject buffers that aren't a C-contiguous run of pairs""" + vertices = np.array([[0, 0], [9, 9], [1, 0], [9, 9], [0, 1], [9, 9], [1, 1], [9, 9]], dtype=np.float64) + edges = np.array([[0, 1], [9, 9], [1, 3], [9, 9]], dtype=np.uintc) + + def triangulation(): + return cdt.Triangulation(cdt.VertexInsertionOrder.AS_PROVIDED, cdt.IntersectingConstraintEdges.NOT_ALLOWED, 0.0) + + for bad in [vertices[::2], vertices[:, ::-1], np.asfortranarray(vertices[::2]), np.zeros((2, 3))]: + with pytest.raises(RuntimeError): + triangulation().insert_vertices(bad) + for insert in [cdt.Triangulation.insert_edges, cdt.Triangulation.conform_to_edges]: + t = triangulation() + t.insert_vertices(vertices[::2].copy()) + with pytest.raises(RuntimeError): + insert(t, edges[::2]) + + for good in [vertices[::2].copy(), vertices[::2].ravel()]: + t = triangulation() + t.insert_vertices(good) + assert [(v.x, v.y) for v in t.vertices][3:] == [(0, 0), (1, 0), (0, 1), (1, 1)], "Wrong vertices inserted" + t.insert_edges(edges[::2].copy()) + assert t.fixed_edges == {cdt.Edge(3, 4), cdt.Edge(4, 6)}, "Wrong edges inserted"