Skip to content

Commit d65ac17

Browse files
authored
Add batched simplices_containing query and Rust default_loss for LearnerND hot loops (#13)
* Add batched simplices_containing query and Rust default_loss Profiling LearnerND with the Rust backend (adaptive PR #493) showed only 14-25% of runtime left inside this extension; two of the remaining Python hotspots are triangulation-shaped and move here: - Triangulation.simplices_containing(point, simplex=None, candidates=None, eps=1e-8): all simplices containing a point in one call. Instead of one barycentric solve per neighbouring simplex (the loop LearnerND.tell_pending runs in Python today, ~14 checks per point in 2D, ~70 in 3D), it solves once for a simplex known to contain the point (the hint, or locate_point), reduces it to the face the point lies on, and looks up the simplices containing that face in the facet index. - default_loss(simplex, values, value_scale=None): the LearnerND default loss (embedded simplex volume) taking the scaled vertex/value arrays directly, signature-compatible with loss_per_simplex. End-to-end with both wired into LearnerND (examples/learnernd_batched_apis.py): 1.18x in 2D and 1.39x in 3D on top of the Rust backend, with identical sampled points. * docs: refresh README benchmarks against adaptive 1.5.0 adaptive 1.5.0 ships the automatic Rust-backend selection, so the usage section now leads with `pip install "adaptive[rust]"` (monkey-patching is only needed for < 1.5.0) and the example points at the release instead of the PR branch. Re-measured all tables against 1.5.0: standalone numbers are unchanged within noise; LearnerND's pure-Python path got faster in 1.5.0, so the honest end-to-end ratio is now 3.3x (was 3.7x). Adds the batched-API numbers (1.17x 2D / 1.40x 3D on top of the backend) to the performance section.
1 parent c6d6de4 commit d65ac17

10 files changed

Lines changed: 506 additions & 21 deletions

File tree

‎README.md‎

Lines changed: 35 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -6,31 +6,39 @@
66
[![License](https://img.shields.io/badge/license-BSD--3--Clause-blue.svg)](LICENSE)
77

88
Fast N-dimensional Delaunay triangulation in Rust with Python bindings (PyO3).
9-
Drop-in replacement for [adaptive](https://github.com/python-adaptive/adaptive)'s `Triangulation` class — **30-300× faster** standalone, **3.7×** end-to-end in `LearnerND` (where adaptive's own Python code dominates).
9+
Drop-in replacement for [adaptive](https://github.com/python-adaptive/adaptive)'s `Triangulation` class — **30-300× faster** standalone, **~3.3×** end-to-end in `LearnerND` (where adaptive's own Python code dominates). Used automatically by adaptive ≥ 1.5 when installed.
1010

1111
## Performance
1212

13-
Measured with the scripts in [`examples/`](examples/), best of 3 for the standalone runs.
13+
Measured with the scripts in [`examples/`](examples/) against adaptive 1.5.0, best of 3.
1414
Absolute times are machine-dependent; the ratios are representative.
1515

1616
### Standalone triangulation (incremental insertion)
1717
| Case | Rust | Python | Speedup |
1818
|---|---:|---:|---:|
19-
| 2D, 1K pts | 18 ms | 731 ms | **40×** |
20-
| 2D, 5K pts | 134 ms | 14,611 ms | **109×** |
21-
| 3D, 500 pts | 32 ms | 3,001 ms | **94×** |
22-
| 3D, 2K pts | 152 ms | 44,262 ms | **291×** |
19+
| 2D, 1K pts | 18 ms | 730 ms | **41×** |
20+
| 2D, 5K pts | 138 ms | 14,442 ms | **105×** |
21+
| 3D, 500 pts | 32 ms | 3,037 ms | **94×** |
22+
| 3D, 2K pts | 150 ms | 44,531 ms | **297×** |
2323

2424
### LearnerND integration (end-to-end, `ring_of_fire` 2D)
2525
| N pts | Learner2D (scipy) | LearnerND (Python) | LearnerND (Rust) |
2626
|---|---:|---:|---:|
27-
| 1,000 | 0.23 s | 0.59 s | **0.16 s** |
28-
| 2,000 | 0.90 s | 1.16 s | **0.32 s** |
29-
| 5,000 | 5.64 s | 2.95 s | **0.81 s** |
27+
| 1,000 | 0.23 s | 0.50 s | **0.16 s** |
28+
| 2,000 | 0.90 s | 1.01 s | **0.32 s** |
29+
| 5,000 | 5.69 s | 2.57 s | **0.79 s** |
3030

31-
LearnerND + Rust is **3.7× faster** than LearnerND + Python, and **7× faster** than Learner2D at 5K points.
31+
LearnerND + Rust is **3.3× faster** than LearnerND + Python, and **7× faster** than Learner2D at 5K points.
3232
The end-to-end ratio is smaller than the standalone one because adaptive's own Python-side loss machinery dominates once the triangulation is fast.
3333

34+
### Batched LearnerND APIs (not yet wired into adaptive)
35+
36+
`simplices_containing` and `default_loss` move two of the remaining `LearnerND`
37+
Python hot loops into Rust. Wired in the way a future adaptive release would use
38+
them ([`examples/learnernd_batched_apis.py`](examples/learnernd_batched_apis.py)),
39+
they add **1.17×** (2D, 3000 pts) to **1.40×** (3D, 1500 pts) on top of the
40+
table above, while sampling identical points.
41+
3442
## Installation
3543

3644
```bash
@@ -59,8 +67,17 @@ print(tri.reference_invariant()) # True
5967

6068
## Usage with adaptive's LearnerND
6169

62-
This is a drop-in replacement for `adaptive`'s built-in triangulation.
63-
Monkey-patch the module to use Rust triangulation everywhere:
70+
Since adaptive 1.5.0 this package is detected and used automatically — no code changes needed:
71+
72+
```bash
73+
pip install "adaptive[rust]"
74+
```
75+
76+
Per learner, the backend can be selected explicitly with
77+
`LearnerND(..., triangulation_backend="auto" | "python" | "rust")`, or globally
78+
with the `ADAPTIVE_TRIANGULATION_BACKEND` environment variable.
79+
80+
For adaptive < 1.5.0, monkey-patch the module instead:
6481

6582
```python
6683
import adaptive_triangulation as at
@@ -92,6 +109,10 @@ tri.circumscribed_circle(simplex) # → (center, radius)
92109
tri.volume(simplex) # Simplex volume
93110
tri.volumes() # All simplex volumes
94111
tri.point_in_simplex(point, simplex) # Containment test
112+
tri.simplices_containing(point) # All simplices containing a point, in one call
113+
# instead of a point_in_simplex loop; pass a known
114+
# containing simplex via simplex=... to skip the
115+
# locate step, or restrict with candidates=...
95116
tri.point_in_circumcircle(pt, simplex) # Circumcircle test
96117
tri.bowyer_watson(pt_index) # Direct Bowyer-Watson
97118
tri.get_opposing_vertices(simplex) # Facet neighbours' opposite vertices
@@ -111,6 +132,8 @@ from adaptive_triangulation import (
111132
point_in_simplex, # Containment test
112133
volume, # Simplex volume
113134
simplex_volume_in_embedding, # Volume in embedding space
135+
default_loss, # LearnerND's default loss (embedded simplex volume),
136+
# signature-compatible with loss_per_simplex
114137
orientation, # Face orientation
115138
)
116139
```

‎examples/learnernd_batched_apis.py‎

Lines changed: 116 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,116 @@
1+
"""Measure the LearnerND speedup from `simplices_containing` and `default_loss`.
2+
3+
This monkeypatches adaptive (>= 1.5.0, with the Rust backend active) to use
4+
the two batched APIs, the way the adaptive-side integration would:
5+
6+
- ``tell_pending`` distributes a pending point over the simplices containing
7+
it with one ``tri.simplices_containing(point, simplex=hint)`` call instead
8+
of a Python loop of per-neighbour ``tri.point_in_simplex`` checks.
9+
- The default loss is the Rust ``default_loss`` instead of the Python
10+
wrapper around ``simplex_volume_in_embedding``.
11+
12+
Requirements:
13+
pip install "adaptive[rust]>=1.5"
14+
"""
15+
16+
from __future__ import annotations
17+
18+
import time
19+
from typing import TYPE_CHECKING
20+
21+
import adaptive_triangulation as at
22+
import numpy as np
23+
from adaptive.learner import learnerND as lnd_mod
24+
from adaptive.learner.learnerND import LearnerND
25+
26+
if TYPE_CHECKING:
27+
from collections.abc import Callable, Sequence
28+
29+
reference_default_loss = lnd_mod.default_loss
30+
31+
32+
def ring_of_fire(xy: Sequence[float]) -> float:
33+
x, y = xy
34+
a, d = 0.2, 0.5
35+
return x + np.exp(-((x**2 + y**2 - d**2) ** 2) / a**4)
36+
37+
38+
def sphere_of_fire(xyz: Sequence[float]) -> float:
39+
x, y, z = xyz
40+
a, d = 0.2, 0.5
41+
return x + np.exp(-((x**2 + y**2 + z**2 - d**2) ** 2) / a**4)
42+
43+
44+
def _add_pending_to_simplex(self: LearnerND, point: tuple, simplex: tuple) -> tuple:
45+
"""`_try_adding_pending_point_to_simplex` without the containment check
46+
(the caller already knows `simplex` contains `point`)."""
47+
if simplex not in self._subtriangulations:
48+
vertices = self.tri.get_vertices(simplex)
49+
self._subtriangulations[simplex] = self._triangulation_class(vertices)
50+
self._pending_to_simplex[point] = simplex
51+
return self._subtriangulations[simplex].add_point(point)
52+
53+
54+
def tell_pending_batched(self: LearnerND, point: tuple, *, simplex: tuple | None = None) -> None:
55+
"""tell_pending with the point_in_simplex loop replaced by one
56+
`simplices_containing` call (the `simplex` argument becomes the hint)."""
57+
point = tuple(point)
58+
if not self.inside_bounds(point):
59+
return
60+
self.pending_points.add(point)
61+
if self.tri is None:
62+
return
63+
containing = self.tri.simplices_containing(point, simplex=simplex)
64+
for simpl in containing:
65+
_, to_add = _add_pending_to_simplex(self, point, simpl)
66+
if to_add is None:
67+
continue
68+
self._update_subsimplex_losses(simpl, to_add)
69+
70+
71+
def run(
72+
func: Callable[[Sequence[float]], float],
73+
bounds: list[tuple[float, float]],
74+
n_points: int,
75+
*,
76+
rust_loss: bool,
77+
batched_tell: bool,
78+
) -> tuple[float, LearnerND]:
79+
lnd_mod.default_loss = at.default_loss if rust_loss else reference_default_loss
80+
learner = LearnerND(func, bounds=bounds)
81+
if batched_tell:
82+
learner.tell_pending = tell_pending_batched.__get__(learner)
83+
t0 = time.perf_counter()
84+
for _ in range(n_points):
85+
points, _ = learner.ask(1)
86+
for p in points:
87+
learner.tell(p, func(p))
88+
elapsed = time.perf_counter() - t0
89+
return elapsed, learner
90+
91+
92+
CONFIGS = [
93+
("baseline (Rust backend, adaptive >= 1.5)", False, False),
94+
("+ rust default_loss", True, False),
95+
("+ simplices_containing tell_pending", False, True),
96+
("+ both", True, True),
97+
]
98+
99+
for name, func, bounds, n in [
100+
("2D ring_of_fire", ring_of_fire, [(-1, 1), (-1, 1)], 3000),
101+
("3D sphere_of_fire", sphere_of_fire, [(-1, 1), (-1, 1), (-1, 1)], 1500),
102+
]:
103+
print(f"\n{name}, {n} points:")
104+
baseline = None
105+
base_points = None
106+
for label, rust_loss, batched_tell in CONFIGS:
107+
best = min(
108+
run(func, bounds, n, rust_loss=rust_loss, batched_tell=batched_tell)[0]
109+
for _ in range(3)
110+
)
111+
_, learner = run(func, bounds, n, rust_loss=rust_loss, batched_tell=batched_tell)
112+
pts = sorted(learner.data)
113+
if baseline is None:
114+
baseline, base_points = best, pts
115+
same = "identical points" if pts == base_points else "DIFFERENT POINTS"
116+
print(f" {label:40s} {best:6.2f}s ({baseline / best:4.2f}x, {same})")

‎python/adaptive_triangulation/__init__.py‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,7 @@
2929
VerticesProxy,
3030
__version__,
3131
circumsphere,
32+
default_loss,
3233
fast_2d_circumcircle,
3334
fast_2d_point_in_simplex,
3435
fast_3d_circumcircle,
@@ -47,6 +48,7 @@
4748
"VerticesProxy",
4849
"__version__",
4950
"circumsphere",
51+
"default_loss",
5052
"fast_2d_circumcircle",
5153
"fast_2d_point_in_simplex",
5254
"fast_3d_circumcircle",

‎python/adaptive_triangulation/__init__.pyi‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,9 @@ from ._rust import (
1616
from ._rust import (
1717
circumsphere as circumsphere,
1818
)
19+
from ._rust import (
20+
default_loss as default_loss,
21+
)
1922
from ._rust import (
2023
fast_2d_circumcircle as fast_2d_circumcircle,
2124
)
@@ -51,6 +54,7 @@ __all__: list[str] = [
5154
"VerticesProxy",
5255
"__version__",
5356
"circumsphere",
57+
"default_loss",
5458
"fast_2d_circumcircle",
5559
"fast_2d_point_in_simplex",
5660
"fast_3d_circumcircle",

‎python/adaptive_triangulation/_rust.pyi‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -94,6 +94,13 @@ class Triangulation:
9494
simplex: SimplexLike,
9595
eps: float = 1e-8,
9696
) -> bool: ...
97+
def simplices_containing(
98+
self,
99+
point: PointLike,
100+
simplex: SimplexLike | None = None,
101+
candidates: Iterable[SimplexLike] | None = None,
102+
eps: float = 1e-8,
103+
) -> list[Simplex]: ...
97104
def bowyer_watson(
98105
self,
99106
pt_index: int,
@@ -104,6 +111,11 @@ class Triangulation:
104111
def convex_invariant(self, vertex: int) -> bool: ...
105112

106113
def circumsphere(points: PointsLike) -> tuple[tuple[float, ...], float]: ...
114+
def default_loss(
115+
simplex: PointsLike,
116+
values: Sequence[float] | PointsLike,
117+
value_scale: float | None = None,
118+
) -> float: ...
107119
def fast_2d_circumcircle(points: PointsLike) -> tuple[tuple[float, ...], float]: ...
108120
def fast_2d_point_in_simplex(point: PointLike, simplex: PointsLike, eps: float = 1e-8) -> bool: ...
109121
def fast_3d_circumcircle(points: PointsLike) -> tuple[tuple[float, ...], float]: ...

‎src/geometry.rs‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -666,6 +666,19 @@ pub fn simplex_volume_in_embedding(vertices: &[Vec<f64>]) -> Result<f64, Geometr
666666
Ok(vol_square.sqrt())
667667
}
668668

669+
/// The loss adaptive's `LearnerND` assigns to a simplex by default: the
670+
/// volume of the simplex embedded in (input + output)-dimensional space,
671+
/// where each vertex is extended with its function value(s). `simplex` and
672+
/// `values` are zipped pairwise, like the reference implementation.
673+
pub fn default_loss(simplex: &[Vec<f64>], values: &[Vec<f64>]) -> Result<f64, GeometryError> {
674+
let embedded: Vec<Vec<f64>> = simplex
675+
.iter()
676+
.zip(values)
677+
.map(|(vertex, value)| vertex.iter().chain(value.iter()).copied().collect())
678+
.collect();
679+
simplex_volume_in_embedding(&embedded)
680+
}
681+
669682
#[cfg(test)]
670683
mod tests {
671684
use super::*;

‎src/lib.rs‎

Lines changed: 50 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@ pub mod py;
88
pub mod tolerances;
99
pub mod triangulation;
1010

11+
use numpy::{PyReadonlyArray1, PyReadonlyArray2};
1112
use pyo3::exceptions::{PyValueError, PyZeroDivisionError};
1213
use pyo3::prelude::*;
1314
use pyo3::types::{PyAny, PyModule, PyTuple};
@@ -47,6 +48,7 @@ fn _rust(m: &Bound<'_, PyModule>) -> PyResult<()> {
4748
m.add_function(wrap_pyfunction!(py_fast_norm, m)?)?;
4849
m.add_function(wrap_pyfunction!(py_volume, m)?)?;
4950
m.add_function(wrap_pyfunction!(py_simplex_volume_in_embedding, m)?)?;
51+
m.add_function(wrap_pyfunction!(py_default_loss, m)?)?;
5052
m.add_function(wrap_pyfunction!(py_orientation, m)?)?;
5153

5254
Ok(())
@@ -176,6 +178,54 @@ fn py_simplex_volume_in_embedding(vertices: &Bound<'_, PyAny>) -> PyResult<f64>
176178
.map_err(|err| PyValueError::new_err(err.to_string()))
177179
}
178180

181+
/// Parse per-vertex function values: one scalar per vertex or one vector per
182+
/// vertex, with bulk-copy fast paths for 1-D and 2-D f64 numpy arrays (what
183+
/// `LearnerND._compute_loss` passes after scaling).
184+
fn parse_values(obj: &Bound<'_, PyAny>) -> PyResult<Vec<Vec<f64>>> {
185+
if let Ok(array) = obj.extract::<PyReadonlyArray1<'_, f64>>() {
186+
return Ok(array.as_array().iter().map(|&value| vec![value]).collect());
187+
}
188+
if let Ok(array) = obj.extract::<PyReadonlyArray2<'_, f64>>() {
189+
return Ok(array
190+
.as_array()
191+
.outer_iter()
192+
.map(|row| row.to_vec())
193+
.collect());
194+
}
195+
let Ok(iter) = obj.try_iter() else {
196+
return Err(PyValueError::new_err(
197+
"Expected one scalar or one vector of floats per vertex",
198+
));
199+
};
200+
let mut values = Vec::new();
201+
for item in iter {
202+
let item = item?;
203+
match item.extract::<f64>() {
204+
Ok(value) => values.push(vec![value]),
205+
Err(_) => values.push(parse_point(&item)?),
206+
}
207+
}
208+
Ok(values)
209+
}
210+
211+
/// The default `LearnerND` loss: the volume of the simplex embedded in
212+
/// (input + output)-dimensional space, where each vertex is extended with
213+
/// its function value(s). Signature-compatible with adaptive's
214+
/// `loss_per_simplex` functions; `value_scale` is accepted and unused, like
215+
/// the reference (`values` arrive already scaled).
216+
#[pyfunction]
217+
#[pyo3(name = "default_loss", signature = (simplex, values, value_scale=None))]
218+
fn py_default_loss(
219+
simplex: &Bound<'_, PyAny>,
220+
values: &Bound<'_, PyAny>,
221+
value_scale: Option<&Bound<'_, PyAny>>,
222+
) -> PyResult<f64> {
223+
let _ = value_scale;
224+
let simplex = parse_points_sized(simplex, "Please provide a 2-dimensional list of points")?;
225+
let values = parse_values(values)?;
226+
geom::default_loss(&simplex, &values).map_err(|err| PyValueError::new_err(err.to_string()))
227+
}
228+
179229
#[pyfunction]
180230
#[pyo3(name = "orientation")]
181231
fn py_orientation(face: &Bound<'_, PyAny>, origin: &Bound<'_, PyAny>) -> PyResult<i32> {

0 commit comments

Comments
 (0)