Optimize cirumsphere and triangulation.py
Joseph Weston authored on 05/05/2020 19:11:46... | ... |
@@ -1,20 +1,45 @@ |
1 |
-import math |
|
2 | 1 |
from collections import Counter |
3 | 2 |
from collections.abc import Iterable, Sized |
4 | 3 |
from itertools import chain, combinations |
5 |
-from math import factorial |
|
4 |
+from math import factorial, sqrt |
|
6 | 5 |
|
7 |
-import numpy as np |
|
8 | 6 |
import scipy.spatial |
7 |
+from numpy import abs as np_abs |
|
8 |
+from numpy import ( |
|
9 |
+ array, |
|
10 |
+ asarray, |
|
11 |
+ average, |
|
12 |
+ concatenate, |
|
13 |
+ dot, |
|
14 |
+ eye, |
|
15 |
+ mean, |
|
16 |
+ ones, |
|
17 |
+ square, |
|
18 |
+ subtract, |
|
19 |
+) |
|
20 |
+from numpy import sum as np_sum |
|
21 |
+from numpy import zeros |
|
22 |
+from numpy.linalg import det as ndet |
|
23 |
+from numpy.linalg import matrix_rank, norm, slogdet, solve |
|
9 | 24 |
|
10 | 25 |
|
11 | 26 |
def fast_norm(v): |
27 |
+ """ Manually take the vector norm for len 2, 3 vectors. Defaults to a square root of the dot product |
|
28 |
+ for larger vectors. |
|
29 |
+ |
|
30 |
+ Note that for large vectors, it is possible for integer overflow to occur. |
|
31 |
+ For instance: |
|
32 |
+ vec = [49024, 59454, 12599, -63721, 18517, 27961] |
|
33 |
+ dot(vec, vec) = -1602973744 |
|
34 |
+ |
|
35 |
+ """ |
|
36 |
+ len_v = len(v) |
|
12 | 37 |
# notice this method can be even more optimised |
13 |
- if len(v) == 2: |
|
14 |
- return math.sqrt(v[0] * v[0] + v[1] * v[1]) |
|
15 |
- if len(v) == 3: |
|
16 |
- return math.sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]) |
|
17 |
- return math.sqrt(np.dot(v, v)) |
|
38 |
+ if len_v == 2: |
|
39 |
+ return sqrt(v[0] * v[0] + v[1] * v[1]) |
|
40 |
+ if len_v == 3: |
|
41 |
+ return sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]) |
|
42 |
+ return sqrt(dot(v, v)) |
|
18 | 43 |
|
19 | 44 |
|
20 | 45 |
def fast_2d_point_in_simplex(point, simplex, eps=1e-8): |
... | ... |
@@ -35,9 +60,9 @@ def point_in_simplex(point, simplex, eps=1e-8): |
35 | 60 |
if len(point) == 2: |
36 | 61 |
return fast_2d_point_in_simplex(point, simplex, eps) |
37 | 62 |
|
38 |
- x0 = np.array(simplex[0], dtype=float) |
|
39 |
- vectors = np.array(simplex[1:], dtype=float) - x0 |
|
40 |
- alpha = np.linalg.solve(vectors.T, point - x0) |
|
63 |
+ x0 = array(simplex[0], dtype=float) |
|
64 |
+ vectors = array(simplex[1:], dtype=float) - x0 |
|
65 |
+ alpha = solve(vectors.T, point - x0) |
|
41 | 66 |
|
42 | 67 |
return all(alpha > -eps) and sum(alpha) < 1 + eps |
43 | 68 |
|
... | ... |
@@ -53,9 +78,9 @@ def fast_2d_circumcircle(points): |
53 | 78 |
Returns |
54 | 79 |
------- |
55 | 80 |
tuple |
56 |
- (center point : tuple(int), radius: int) |
|
81 |
+ (center point : tuple(float), radius: float) |
|
57 | 82 |
""" |
58 |
- points = np.array(points) |
|
83 |
+ points = array(points) |
|
59 | 84 |
# transform to relative coordinates |
60 | 85 |
pts = points[1:] - points[0] |
61 | 86 |
|
... | ... |
@@ -73,7 +98,7 @@ def fast_2d_circumcircle(points): |
73 | 98 |
# compute center |
74 | 99 |
x = dx / a |
75 | 100 |
y = dy / a |
76 |
- radius = math.sqrt(x * x + y * y) # radius = norm([x, y]) |
|
101 |
+ radius = sqrt(x * x + y * y) # radius = norm([x, y]) |
|
77 | 102 |
|
78 | 103 |
return (x + points[0][0], y + points[0][1]), radius |
79 | 104 |
|
... | ... |
@@ -89,9 +114,9 @@ def fast_3d_circumcircle(points): |
89 | 114 |
Returns |
90 | 115 |
------- |
91 | 116 |
tuple |
92 |
- (center point : tuple(int), radius: int) |
|
117 |
+ (center point : tuple(float), radius: float) |
|
93 | 118 |
""" |
94 |
- points = np.array(points) |
|
119 |
+ points = array(points) |
|
95 | 120 |
pts = points[1:] - points[0] |
96 | 121 |
|
97 | 122 |
(x1, y1, z1), (x2, y2, z2), (x3, y3, z3) = pts |
... | ... |
@@ -119,17 +144,36 @@ def fast_3d_circumcircle(points): |
119 | 144 |
|
120 | 145 |
|
121 | 146 |
def fast_det(matrix): |
122 |
- matrix = np.asarray(matrix, dtype=float) |
|
147 |
+ matrix = asarray(matrix, dtype=float) |
|
123 | 148 |
if matrix.shape == (2, 2): |
124 | 149 |
return matrix[0][0] * matrix[1][1] - matrix[1][0] * matrix[0][1] |
125 | 150 |
elif matrix.shape == (3, 3): |
126 | 151 |
a, b, c, d, e, f, g, h, i = matrix.ravel() |
127 | 152 |
return a * (e * i - f * h) - b * (d * i - f * g) + c * (d * h - e * g) |
128 | 153 |
else: |
129 |
- return np.linalg.det(matrix) |
|
154 |
+ return ndet(matrix) |
|
130 | 155 |
|
131 | 156 |
|
132 | 157 |
def circumsphere(pts): |
158 |
+ """Compute the center and radius of a N dimension sphere which touches each point in pts. |
|
159 |
+ |
|
160 |
+ Parameters |
|
161 |
+ ---------- |
|
162 |
+ pts : array-like, of shape (N-dim + 1, N-dim) |
|
163 |
+ The points for which we would like to compute a circumsphere. |
|
164 |
+ |
|
165 |
+ Returns |
|
166 |
+ ------- |
|
167 |
+ center : tuple of floats of size N-dim |
|
168 |
+ radius : a positive float |
|
169 |
+ A valid center and radius, if a circumsphere is possible, and no points are repeated. |
|
170 |
+ If points are repeated, or a circumsphere is not possible, will return nans, and a |
|
171 |
+ ZeroDivisionError may occur. |
|
172 |
+ Will fail for matrices which are not (N-dim + 1, N-dim) in size due to non-square determinants: |
|
173 |
+ will raise numpy.linalg.LinAlgError. |
|
174 |
+ May fail for points that are integers (due to 32bit integer overflow). |
|
175 |
+ """ |
|
176 |
+ |
|
133 | 177 |
dim = len(pts) - 1 |
134 | 178 |
if dim == 2: |
135 | 179 |
return fast_2d_circumcircle(pts) |
... | ... |
@@ -137,20 +181,23 @@ def circumsphere(pts): |
137 | 181 |
return fast_3d_circumcircle(pts) |
138 | 182 |
|
139 | 183 |
# Modified method from http://mathworld.wolfram.com/Circumsphere.html |
140 |
- mat = [[np.sum(np.square(pt)), *pt, 1] for pt in pts] |
|
141 |
- |
|
142 |
- center = [] |
|
184 |
+ mat = array([[np_sum(square(pt)), *pt, 1] for pt in pts]) |
|
185 |
+ center = zeros(dim) |
|
186 |
+ a = 1 / (2 * ndet(mat[:, 1:])) |
|
187 |
+ factor = a |
|
188 |
+ # Use ind to index into the matrix columns |
|
189 |
+ ind = ones((dim + 2,), bool) |
|
143 | 190 |
for i in range(1, len(pts)): |
144 |
- r = np.delete(mat, i, 1) |
|
145 |
- factor = (-1) ** (i + 1) |
|
146 |
- center.append(factor * fast_det(r)) |
|
147 |
- |
|
148 |
- a = fast_det(np.delete(mat, 0, 1)) |
|
149 |
- center = [x / (2 * a) for x in center] |
|
191 |
+ ind[i - 1] = True |
|
192 |
+ ind[i] = False |
|
193 |
+ center[i - 1] = factor * ndet(mat[:, ind]) |
|
194 |
+ factor *= -1 |
|
150 | 195 |
|
196 |
+ # Use subtract as we don't know the type of x0. |
|
151 | 197 |
x0 = pts[0] |
152 |
- vec = np.subtract(center, x0) |
|
153 |
- radius = fast_norm(vec) |
|
198 |
+ vec = subtract(center, x0) |
|
199 |
+ # Vector norm. |
|
200 |
+ radius = sqrt(dot(vec, vec)) |
|
154 | 201 |
|
155 | 202 |
return tuple(center), radius |
156 | 203 |
|
... | ... |
@@ -174,8 +221,8 @@ def orientation(face, origin): |
174 | 221 |
If two points lie on the same side of the face, the orientation will |
175 | 222 |
be equal, if they lie on the other side of the face, it will be negated. |
176 | 223 |
""" |
177 |
- vectors = np.array(face) |
|
178 |
- sign, logdet = np.linalg.slogdet(vectors - origin) |
|
224 |
+ vectors = array(face) |
|
225 |
+ sign, logdet = slogdet(vectors - origin) |
|
179 | 226 |
if logdet < -50: # assume it to be zero when it's close to zero |
180 | 227 |
return 0 |
181 | 228 |
return sign |
... | ... |
@@ -198,7 +245,7 @@ def simplex_volume_in_embedding(vertices) -> float: |
198 | 245 |
|
199 | 246 |
Returns |
200 | 247 |
------- |
201 |
- volume : int |
|
248 |
+ volume : float |
|
202 | 249 |
the volume of the simplex with given vertices. |
203 | 250 |
|
204 | 251 |
Raises |
... | ... |
@@ -210,20 +257,20 @@ def simplex_volume_in_embedding(vertices) -> float: |
210 | 257 |
# Implements http://mathworld.wolfram.com/Cayley-MengerDeterminant.html |
211 | 258 |
# Modified from https://codereview.stackexchange.com/questions/77593/calculating-the-volume-of-a-tetrahedron |
212 | 259 |
|
213 |
- vertices = np.asarray(vertices, dtype=float) |
|
260 |
+ vertices = asarray(vertices, dtype=float) |
|
214 | 261 |
dim = len(vertices[0]) |
215 | 262 |
if dim == 2: |
216 | 263 |
# Heron's formula |
217 | 264 |
a, b, c = scipy.spatial.distance.pdist(vertices, metric="euclidean") |
218 | 265 |
s = 0.5 * (a + b + c) |
219 |
- return math.sqrt(s * (s - a) * (s - b) * (s - c)) |
|
266 |
+ return sqrt(s * (s - a) * (s - b) * (s - c)) |
|
220 | 267 |
|
221 | 268 |
# β_ij = |v_i - v_k|² |
222 | 269 |
sq_dists = scipy.spatial.distance.pdist(vertices, metric="sqeuclidean") |
223 | 270 |
|
224 | 271 |
# Add border while compressed |
225 | 272 |
num_verts = scipy.spatial.distance.num_obs_y(sq_dists) |
226 |
- bordered = np.concatenate((np.ones(num_verts), sq_dists)) |
|
273 |
+ bordered = concatenate((ones(num_verts), sq_dists)) |
|
227 | 274 |
|
228 | 275 |
# Make matrix and find volume |
229 | 276 |
sq_dists_mat = scipy.spatial.distance.squareform(bordered) |
... | ... |
@@ -232,11 +279,11 @@ def simplex_volume_in_embedding(vertices) -> float: |
232 | 279 |
vol_square = fast_det(sq_dists_mat) / coeff |
233 | 280 |
|
234 | 281 |
if vol_square < 0: |
235 |
- if abs(vol_square) < 1e-15: |
|
282 |
+ if vol_square > -1e-15: |
|
236 | 283 |
return 0 |
237 | 284 |
raise ValueError("Provided vertices do not form a simplex") |
238 | 285 |
|
239 |
- return np.sqrt(vol_square) |
|
286 |
+ return sqrt(vol_square) |
|
240 | 287 |
|
241 | 288 |
|
242 | 289 |
class Triangulation: |
... | ... |
@@ -287,8 +334,8 @@ class Triangulation: |
287 | 334 |
raise ValueError("Please provide at least one simplex") |
288 | 335 |
|
289 | 336 |
coords = list(map(tuple, coords)) |
290 |
- vectors = np.subtract(coords[1:], coords[0]) |
|
291 |
- if np.linalg.matrix_rank(vectors) < dim: |
|
337 |
+ vectors = subtract(coords[1:], coords[0]) |
|
338 |
+ if matrix_rank(vectors) < dim: |
|
292 | 339 |
raise ValueError( |
293 | 340 |
"Initial simplex has zero volumes " |
294 | 341 |
"(the points are linearly dependent)" |
... | ... |
@@ -338,9 +385,9 @@ class Triangulation: |
338 | 385 |
if len(simplex) != self.dim + 1: |
339 | 386 |
# We are checking whether point belongs to a face. |
340 | 387 |
simplex = self.containing(simplex).pop() |
341 |
- x0 = np.array(self.vertices[simplex[0]]) |
|
342 |
- vectors = np.array(self.get_vertices(simplex[1:])) - x0 |
|
343 |
- alpha = np.linalg.solve(vectors.T, point - x0) |
|
388 |
+ x0 = array(self.vertices[simplex[0]]) |
|
389 |
+ vectors = array(self.get_vertices(simplex[1:])) - x0 |
|
390 |
+ alpha = solve(vectors.T, point - x0) |
|
344 | 391 |
if any(alpha < -eps) or sum(alpha) > 1 + eps: |
345 | 392 |
return [] |
346 | 393 |
|
... | ... |
@@ -403,7 +450,7 @@ class Triangulation: |
403 | 450 |
# we do not really need the center, we only need a point that is |
404 | 451 |
# guaranteed to lie strictly within the hull |
405 | 452 |
hull_points = self.get_vertices(self.hull) |
406 |
- pt_center = np.average(hull_points, axis=0) |
|
453 |
+ pt_center = average(hull_points, axis=0) |
|
407 | 454 |
|
408 | 455 |
pt_index = len(self.vertices) |
409 | 456 |
self.vertices.append(new_vertex) |
... | ... |
@@ -447,7 +494,7 @@ class Triangulation: |
447 | 494 |
tuple (center point, radius) |
448 | 495 |
The center and radius of the circumscribed circle |
449 | 496 |
""" |
450 |
- pts = np.dot(self.get_vertices(simplex), transform) |
|
497 |
+ pts = dot(self.get_vertices(simplex), transform) |
|
451 | 498 |
return circumsphere(pts) |
452 | 499 |
|
453 | 500 |
def point_in_cicumcircle(self, pt_index, simplex, transform): |
... | ... |
@@ -455,13 +502,13 @@ class Triangulation: |
455 | 502 |
eps = 1e-8 |
456 | 503 |
|
457 | 504 |
center, radius = self.circumscribed_circle(simplex, transform) |
458 |
- pt = np.dot(self.get_vertices([pt_index]), transform)[0] |
|
505 |
+ pt = dot(self.get_vertices([pt_index]), transform)[0] |
|
459 | 506 |
|
460 |
- return np.linalg.norm(center - pt) < (radius * (1 + eps)) |
|
507 |
+ return norm(center - pt) < (radius * (1 + eps)) |
|
461 | 508 |
|
462 | 509 |
@property |
463 | 510 |
def default_transform(self): |
464 |
- return np.eye(self.dim) |
|
511 |
+ return eye(self.dim) |
|
465 | 512 |
|
466 | 513 |
def bowyer_watson(self, pt_index, containing_simplex=None, transform=None): |
467 | 514 |
"""Modified Bowyer-Watson point adding algorithm. |
... | ... |
@@ -532,9 +579,9 @@ class Triangulation: |
532 | 579 |
volume is only dependent on the shape of the simplex and not on the |
533 | 580 |
absolute size. Due to the weird scaling, the only use of this method |
534 | 581 |
is to check that a simplex is almost flat.""" |
535 |
- vertices = np.array(self.get_vertices(simplex)) |
|
582 |
+ vertices = array(self.get_vertices(simplex)) |
|
536 | 583 |
vectors = vertices[1:] - vertices[0] |
537 |
- average_edge_length = np.mean(np.abs(vectors)) |
|
584 |
+ average_edge_length = mean(np_abs(vectors)) |
|
538 | 585 |
return self.volume(simplex) / (average_edge_length ** self.dim) |
539 | 586 |
|
540 | 587 |
def add_point(self, point, simplex=None, transform=None): |
... | ... |
@@ -587,8 +634,8 @@ class Triangulation: |
587 | 634 |
return self.bowyer_watson(pt_index, actual_simplex, transform) |
588 | 635 |
|
589 | 636 |
def volume(self, simplex): |
590 |
- prefactor = np.math.factorial(self.dim) |
|
591 |
- vertices = np.array(self.get_vertices(simplex)) |
|
637 |
+ prefactor = factorial(self.dim) |
|
638 |
+ vertices = array(self.get_vertices(simplex)) |
|
592 | 639 |
vectors = vertices[1:] - vertices[0] |
593 | 640 |
return float(abs(fast_det(vectors)) / prefactor) |
594 | 641 |
|
... | ... |
@@ -60,3 +60,40 @@ def test_triangulation_find_opposing_vertices_raises_if_simplex_is_invalid(): |
60 | 60 |
|
61 | 61 |
with pytest.raises(ValueError): |
62 | 62 |
tri.get_opposing_vertices((2, 3, 5)) |
63 |
+ |
|
64 |
+ |
|
65 |
+def test_circumsphere(): |
|
66 |
+ from adaptive.learner.triangulation import circumsphere, fast_norm |
|
67 |
+ from numpy import allclose |
|
68 |
+ from numpy.random import normal, uniform |
|
69 |
+ |
|
70 |
+ def generate_random_sphere_points(dim, radius=0): |
|
71 |
+ """ Refer to https://math.stackexchange.com/a/1585996 """ |
|
72 |
+ |
|
73 |
+ vec = [None] * (dim + 1) |
|
74 |
+ center = uniform(-100, 100, dim) |
|
75 |
+ radius = uniform(1.0, 100.0) if radius == 0 else radius |
|
76 |
+ for i in range(dim + 1): |
|
77 |
+ points = normal(0, size=dim) |
|
78 |
+ x = fast_norm(points) |
|
79 |
+ points = points / x * radius |
|
80 |
+ vec[i] = tuple(points + center) |
|
81 |
+ |
|
82 |
+ return radius, center, vec |
|
83 |
+ |
|
84 |
+ center_diff_err = "Calculated center (%s) differs from true center (%s)\n" |
|
85 |
+ for dim in range(2, 10): |
|
86 |
+ radius, center, points = generate_random_sphere_points(dim) |
|
87 |
+ circ_center, circ_radius = circumsphere(points) |
|
88 |
+ err_msg = "" |
|
89 |
+ if not allclose(circ_center, center): |
|
90 |
+ err_msg += center_diff_err % ( |
|
91 |
+ ", ".join([str(x) for x in circ_center]), |
|
92 |
+ ", ".join([str(x) for x in center]), |
|
93 |
+ ) |
|
94 |
+ if not allclose(radius, circ_radius): |
|
95 |
+ err_msg += "Calculated radius {} differs from true radius {}".format( |
|
96 |
+ circ_radius, radius, |
|
97 |
+ ) |
|
98 |
+ if err_msg: |
|
99 |
+ raise AssertionError(err_msg) |