Numpy triangles array

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/)
 

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")

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))