Skip to content

Commit 2a6c247

Browse files
basnijholtbasnijholt
committed
Fix index out of bounds error in high-dimensional circumsphere calculation
- Fixed incorrect indexing in Cramer's rule implementation for N-dimensional circumsphere - Added direct triangulation performance benchmarks showing 11-16x speedup - Verified Rust implementation works correctly in release mode Performance results: - Triangulation creation/addition: 11.64x speedup - Point location queries: 16.13x speedup
1 parent ece68e6 commit 2a6c247

File tree

2 files changed

+223
-1
lines changed

2 files changed

+223
-1
lines changed

‎src/geometry.rs

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -224,15 +224,18 @@ pub fn circumsphere(points: &[Vec<f64>]) -> Result<(Vec<f64>, f64), GeometryErro
224224
// Build matrix with appropriate column replaced
225225
let mut num_mat = DMatrix::<f64>::zeros(dim + 1, dim + 1);
226226
for i in 0..=dim {
227-
for j in 0..=dim {
227+
for j in 0..dim {
228228
if j == coord_idx {
229229
num_mat[(i, j)] = mat[i][0]; // Replace with sum of squares
230230
} else if j < coord_idx {
231231
num_mat[(i, j)] = mat[i][j + 1];
232232
} else {
233+
// j > coord_idx
233234
num_mat[(i, j)] = mat[i][j + 2];
234235
}
235236
}
237+
// Last column is always 1
238+
num_mat[(i, dim)] = 1.0;
236239
}
237240

238241
center[coord_idx] = factor * num_mat.determinant();

‎test_triangulation_direct.py

Lines changed: 219 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,219 @@
1+
#!/usr/bin/env python3
2+
"""Direct comparison of Python vs Rust triangulation performance."""
3+
4+
import time
5+
6+
import numpy as np
7+
from adaptive_rust import RustTriangulation
8+
9+
from adaptive.learner.triangulation import Triangulation as PyTriangulation
10+
11+
12+
def benchmark_triangulation_creation(n_points=100, dim=2):
13+
"""Benchmark triangulation creation and point addition."""
14+
print(f"\nBenchmarking with {n_points} points in {dim}D:")
15+
print("-" * 50)
16+
17+
# Generate random points
18+
np.random.seed(42)
19+
initial_points = np.random.rand(dim + 1, dim) # Initial simplex
20+
additional_points = np.random.rand(n_points, dim)
21+
22+
# Python implementation
23+
print("Python Triangulation:")
24+
start = time.perf_counter()
25+
py_tri = PyTriangulation(initial_points)
26+
py_create_time = time.perf_counter() - start
27+
print(f" Creation time: {py_create_time * 1000:.3f}ms")
28+
29+
start = time.perf_counter()
30+
for _i, pt in enumerate(additional_points):
31+
try:
32+
py_tri.add_point(pt)
33+
except ValueError:
34+
pass # Point already exists
35+
py_add_time = time.perf_counter() - start
36+
print(f" Adding {n_points} points: {py_add_time * 1000:.3f}ms")
37+
print(f" Final simplices: {len(py_tri.simplices)}")
38+
py_total = py_create_time + py_add_time
39+
40+
# Rust implementation
41+
print("\nRust Triangulation:")
42+
start = time.perf_counter()
43+
rust_tri = RustTriangulation(initial_points)
44+
rust_create_time = time.perf_counter() - start
45+
print(f" Creation time: {rust_create_time * 1000:.3f}ms")
46+
47+
start = time.perf_counter()
48+
for _i, pt in enumerate(additional_points):
49+
try:
50+
rust_tri.add_point(pt.tolist())
51+
except ValueError:
52+
pass # Point already exists
53+
rust_add_time = time.perf_counter() - start
54+
print(f" Adding {n_points} points: {rust_add_time * 1000:.3f}ms")
55+
print(f" Final simplices: {len(rust_tri.get_simplices())}")
56+
rust_total = rust_create_time + rust_add_time
57+
58+
# Speedup
59+
speedup = py_total / rust_total
60+
print(f"\nSpeedup: {speedup:.2f}x")
61+
62+
return py_total, rust_total, speedup
63+
64+
65+
def benchmark_point_location(n_points=100, n_queries=1000):
66+
"""Benchmark point location performance."""
67+
print(f"\nBenchmarking point location with {n_points} points, {n_queries} queries:")
68+
print("-" * 50)
69+
70+
# Create triangulations with points
71+
np.random.seed(42)
72+
dim = 2
73+
initial_points = np.random.rand(dim + 1, dim)
74+
additional_points = np.random.rand(n_points, dim)
75+
76+
# Build Python triangulation
77+
py_tri = PyTriangulation(initial_points)
78+
for pt in additional_points:
79+
try:
80+
py_tri.add_point(pt)
81+
except ValueError:
82+
pass
83+
84+
# Build Rust triangulation
85+
rust_tri = RustTriangulation(initial_points)
86+
for pt in additional_points:
87+
try:
88+
rust_tri.add_point(pt.tolist())
89+
except ValueError:
90+
pass
91+
92+
# Generate query points
93+
query_points = np.random.rand(n_queries, dim)
94+
95+
# Python point location
96+
print("Python point location:")
97+
start = time.perf_counter()
98+
py_found = 0
99+
for pt in query_points:
100+
simplex = py_tri.locate_point(pt)
101+
if simplex:
102+
py_found += 1
103+
py_time = time.perf_counter() - start
104+
print(f" Time: {py_time * 1000:.3f}ms")
105+
print(f" Points found: {py_found}/{n_queries}")
106+
107+
# Rust point location
108+
print("\nRust point location:")
109+
start = time.perf_counter()
110+
rust_found = 0
111+
for pt in query_points:
112+
simplex = rust_tri.locate_point(pt.tolist())
113+
if simplex:
114+
rust_found += 1
115+
rust_time = time.perf_counter() - start
116+
print(f" Time: {rust_time * 1000:.3f}ms")
117+
print(f" Points found: {rust_found}/{n_queries}")
118+
119+
speedup = py_time / rust_time
120+
print(f"\nSpeedup: {speedup:.2f}x")
121+
122+
return py_time, rust_time, speedup
123+
124+
125+
def benchmark_scaling():
126+
"""Benchmark how performance scales with number of points."""
127+
print("\n" + "=" * 60)
128+
print("SCALING BENCHMARK")
129+
print("=" * 60)
130+
131+
point_counts = [50, 100, 200, 500, 1000]
132+
creation_speedups = []
133+
location_speedups = []
134+
135+
for n in point_counts:
136+
_, _, speedup_create = benchmark_triangulation_creation(n, dim=2)
137+
creation_speedups.append(speedup_create)
138+
139+
_, _, speedup_locate = benchmark_point_location(n, n_queries=100)
140+
location_speedups.append(speedup_locate)
141+
142+
print("\n" + "=" * 60)
143+
print("SUMMARY")
144+
print("=" * 60)
145+
print("\nCreation/Addition Speedups:")
146+
for n, s in zip(point_counts, creation_speedups):
147+
print(f" {n:4d} points: {s:.2f}x")
148+
149+
print("\nPoint Location Speedups:")
150+
for n, s in zip(point_counts, location_speedups):
151+
print(f" {n:4d} points: {s:.2f}x")
152+
153+
avg_creation = sum(creation_speedups) / len(creation_speedups)
154+
avg_location = sum(location_speedups) / len(location_speedups)
155+
print("\nAverage speedups:")
156+
print(f" Creation/Addition: {avg_creation:.2f}x")
157+
print(f" Point Location: {avg_location:.2f}x")
158+
159+
160+
def benchmark_high_dimensions():
161+
"""Benchmark performance in higher dimensions."""
162+
print("\n" + "=" * 60)
163+
print("HIGH-DIMENSIONAL BENCHMARK")
164+
print("=" * 60)
165+
166+
for dim in [2, 3, 4, 5]:
167+
n_points = 50 # Fewer points for higher dimensions
168+
print(f"\n{dim}D space with {n_points} points:")
169+
print("-" * 40)
170+
171+
# Generate points
172+
np.random.seed(42)
173+
initial_points = np.random.rand(dim + 1, dim)
174+
additional_points = np.random.rand(n_points, dim)
175+
176+
# Python
177+
start = time.perf_counter()
178+
py_tri = PyTriangulation(initial_points)
179+
for pt in additional_points:
180+
try:
181+
py_tri.add_point(pt)
182+
except ValueError:
183+
pass
184+
py_time = time.perf_counter() - start
185+
print(f" Python: {py_time * 1000:.3f}ms")
186+
187+
# Rust
188+
start = time.perf_counter()
189+
rust_tri = RustTriangulation(initial_points)
190+
for pt in additional_points:
191+
try:
192+
rust_tri.add_point(pt.tolist())
193+
except ValueError:
194+
pass
195+
rust_time = time.perf_counter() - start
196+
print(f" Rust: {rust_time * 1000:.3f}ms")
197+
198+
speedup = py_time / rust_time
199+
print(f" Speedup: {speedup:.2f}x")
200+
201+
202+
if __name__ == "__main__":
203+
print("=" * 60)
204+
print("DIRECT TRIANGULATION PERFORMANCE COMPARISON")
205+
print("=" * 60)
206+
207+
# Basic benchmarks
208+
benchmark_triangulation_creation(100, dim=2)
209+
benchmark_point_location(100, 1000)
210+
211+
# Scaling benchmark
212+
benchmark_scaling()
213+
214+
# High-dimensional benchmark
215+
benchmark_high_dimensions()
216+
217+
print("\n" + "=" * 60)
218+
print("ALL BENCHMARKS COMPLETED!")
219+
print("=" * 60)

0 commit comments

Comments
 (0)