To speed up circumsphere, I directly imported the relevant functions, and I also removed the fast_det calls, as these end up defaulting to numpy's linalg.det anyways. I also avoid array allocations by getting rid of np.delete() calls. I made a few more tweaks, particularly in dealing with applying the (-1) multiplier and the constant factor a.
I also added general optimizations by directly importing all numpy functions - for hot loops, this means that Python doesn't have to look up the function repeatedly.
... | ... |
@@ -1,20 +1,22 @@ |
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 | 4 |
from math import factorial |
6 |
- |
|
7 |
-import numpy as np |
|
8 | 5 |
import scipy.spatial |
9 | 6 |
|
7 |
+from numpy import square, zeros, subtract, array, ones, dot, asarray, concatenate, average, eye, mean, abs, sqrt |
|
8 |
+from numpy import sum as nsum |
|
9 |
+from numpy.linalg import det as ndet |
|
10 |
+from numpy.linalg import slogdet, solve, matrix_rank, norm |
|
11 |
+ |
|
10 | 12 |
|
11 | 13 |
def fast_norm(v): |
12 | 14 |
# notice this method can be even more optimised |
13 | 15 |
if len(v) == 2: |
14 |
- return math.sqrt(v[0] * v[0] + v[1] * v[1]) |
|
16 |
+ return sqrt(v[0] * v[0] + v[1] * v[1]) |
|
15 | 17 |
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)) |
|
18 |
+ return sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]) |
|
19 |
+ return sqrt(dot(v, v)) |
|
18 | 20 |
|
19 | 21 |
|
20 | 22 |
def fast_2d_point_in_simplex(point, simplex, eps=1e-8): |
... | ... |
@@ -32,12 +34,13 @@ def fast_2d_point_in_simplex(point, simplex, eps=1e-8): |
32 | 34 |
|
33 | 35 |
|
34 | 36 |
def point_in_simplex(point, simplex, eps=1e-8): |
37 |
+ # simplex is list |
|
35 | 38 |
if len(point) == 2: |
36 | 39 |
return fast_2d_point_in_simplex(point, simplex, eps) |
37 | 40 |
|
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) |
|
41 |
+ x0 = array(simplex[0], dtype=float) |
|
42 |
+ vectors = array(simplex[1:], dtype=float) - x0 |
|
43 |
+ alpha = solve(vectors.T, point - x0) |
|
41 | 44 |
|
42 | 45 |
return all(alpha > -eps) and sum(alpha) < 1 + eps |
43 | 46 |
|
... | ... |
@@ -55,7 +58,7 @@ def fast_2d_circumcircle(points): |
55 | 58 |
tuple |
56 | 59 |
(center point : tuple(int), radius: int) |
57 | 60 |
""" |
58 |
- points = np.array(points) |
|
61 |
+ points = array(points) |
|
59 | 62 |
# transform to relative coordinates |
60 | 63 |
pts = points[1:] - points[0] |
61 | 64 |
|
... | ... |
@@ -73,7 +76,7 @@ def fast_2d_circumcircle(points): |
73 | 76 |
# compute center |
74 | 77 |
x = dx / a |
75 | 78 |
y = dy / a |
76 |
- radius = math.sqrt(x * x + y * y) # radius = norm([x, y]) |
|
79 |
+ radius = sqrt(x * x + y * y) # radius = norm([x, y]) |
|
77 | 80 |
|
78 | 81 |
return (x + points[0][0], y + points[0][1]), radius |
79 | 82 |
|
... | ... |
@@ -91,7 +94,7 @@ def fast_3d_circumcircle(points): |
91 | 94 |
tuple |
92 | 95 |
(center point : tuple(int), radius: int) |
93 | 96 |
""" |
94 |
- points = np.array(points) |
|
97 |
+ points = array(points) |
|
95 | 98 |
pts = points[1:] - points[0] |
96 | 99 |
|
97 | 100 |
(x1, y1, z1), (x2, y2, z2), (x3, y3, z3) = pts |
... | ... |
@@ -119,14 +122,14 @@ def fast_3d_circumcircle(points): |
119 | 122 |
|
120 | 123 |
|
121 | 124 |
def fast_det(matrix): |
122 |
- matrix = np.asarray(matrix, dtype=float) |
|
125 |
+ matrix = asarray(matrix, dtype=float) |
|
123 | 126 |
if matrix.shape == (2, 2): |
124 | 127 |
return matrix[0][0] * matrix[1][1] - matrix[1][0] * matrix[0][1] |
125 | 128 |
elif matrix.shape == (3, 3): |
126 | 129 |
a, b, c, d, e, f, g, h, i = matrix.ravel() |
127 | 130 |
return a * (e * i - f * h) - b * (d * i - f * g) + c * (d * h - e * g) |
128 | 131 |
else: |
129 |
- return np.linalg.det(matrix) |
|
132 |
+ return ndet(matrix) |
|
130 | 133 |
|
131 | 134 |
|
132 | 135 |
def circumsphere(pts): |
... | ... |
@@ -137,20 +140,20 @@ def circumsphere(pts): |
137 | 140 |
return fast_3d_circumcircle(pts) |
138 | 141 |
|
139 | 142 |
# 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 = [] |
|
143 |
+ mat = array([[nsum(square(pt)), *pt, 1] for pt in pts]) |
|
144 |
+ center = zeros(dim) |
|
145 |
+ a = 1 / (2 * ndet(mat[:, 1:])) |
|
146 |
+ factor = a |
|
147 |
+ ind = ones((dim + 2,), bool) |
|
143 | 148 |
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] |
|
149 |
+ ind[i - 1] = True |
|
150 |
+ ind[i] = False |
|
151 |
+ center[i - 1] = factor * ndet(mat[:, ind]) |
|
152 |
+ factor *= -1 |
|
150 | 153 |
|
151 | 154 |
x0 = pts[0] |
152 |
- vec = np.subtract(center, x0) |
|
153 |
- radius = fast_norm(vec) |
|
155 |
+ vec = subtract(center, x0) |
|
156 |
+ radius = sqrt(dot(vec, vec)) |
|
154 | 157 |
|
155 | 158 |
return tuple(center), radius |
156 | 159 |
|
... | ... |
@@ -174,8 +177,8 @@ def orientation(face, origin): |
174 | 177 |
If two points lie on the same side of the face, the orientation will |
175 | 178 |
be equal, if they lie on the other side of the face, it will be negated. |
176 | 179 |
""" |
177 |
- vectors = np.array(face) |
|
178 |
- sign, logdet = np.linalg.slogdet(vectors - origin) |
|
180 |
+ vectors = array(face) |
|
181 |
+ sign, logdet = slogdet(vectors - origin) |
|
179 | 182 |
if logdet < -50: # assume it to be zero when it's close to zero |
180 | 183 |
return 0 |
181 | 184 |
return sign |
... | ... |
@@ -210,20 +213,20 @@ def simplex_volume_in_embedding(vertices) -> float: |
210 | 213 |
# Implements http://mathworld.wolfram.com/Cayley-MengerDeterminant.html |
211 | 214 |
# Modified from https://codereview.stackexchange.com/questions/77593/calculating-the-volume-of-a-tetrahedron |
212 | 215 |
|
213 |
- vertices = np.asarray(vertices, dtype=float) |
|
216 |
+ vertices = asarray(vertices, dtype=float) |
|
214 | 217 |
dim = len(vertices[0]) |
215 | 218 |
if dim == 2: |
216 | 219 |
# Heron's formula |
217 | 220 |
a, b, c = scipy.spatial.distance.pdist(vertices, metric="euclidean") |
218 | 221 |
s = 0.5 * (a + b + c) |
219 |
- return math.sqrt(s * (s - a) * (s - b) * (s - c)) |
|
222 |
+ return sqrt(s * (s - a) * (s - b) * (s - c)) |
|
220 | 223 |
|
221 | 224 |
# β_ij = |v_i - v_k|² |
222 | 225 |
sq_dists = scipy.spatial.distance.pdist(vertices, metric="sqeuclidean") |
223 | 226 |
|
224 | 227 |
# Add border while compressed |
225 | 228 |
num_verts = scipy.spatial.distance.num_obs_y(sq_dists) |
226 |
- bordered = np.concatenate((np.ones(num_verts), sq_dists)) |
|
229 |
+ bordered = concatenate((ones(num_verts), sq_dists)) |
|
227 | 230 |
|
228 | 231 |
# Make matrix and find volume |
229 | 232 |
sq_dists_mat = scipy.spatial.distance.squareform(bordered) |
... | ... |
@@ -236,7 +239,7 @@ def simplex_volume_in_embedding(vertices) -> float: |
236 | 239 |
return 0 |
237 | 240 |
raise ValueError("Provided vertices do not form a simplex") |
238 | 241 |
|
239 |
- return np.sqrt(vol_square) |
|
242 |
+ return sqrt(vol_square) |
|
240 | 243 |
|
241 | 244 |
|
242 | 245 |
class Triangulation: |
... | ... |
@@ -287,8 +290,8 @@ class Triangulation: |
287 | 290 |
raise ValueError("Please provide at least one simplex") |
288 | 291 |
|
289 | 292 |
coords = list(map(tuple, coords)) |
290 |
- vectors = np.subtract(coords[1:], coords[0]) |
|
291 |
- if np.linalg.matrix_rank(vectors) < dim: |
|
293 |
+ vectors = subtract(coords[1:], coords[0]) |
|
294 |
+ if matrix_rank(vectors) < dim: |
|
292 | 295 |
raise ValueError( |
293 | 296 |
"Initial simplex has zero volumes " |
294 | 297 |
"(the points are linearly dependent)" |
... | ... |
@@ -338,9 +341,9 @@ class Triangulation: |
338 | 341 |
if len(simplex) != self.dim + 1: |
339 | 342 |
# We are checking whether point belongs to a face. |
340 | 343 |
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) |
|
344 |
+ x0 = array(self.vertices[simplex[0]]) |
|
345 |
+ vectors = array(self.get_vertices(simplex[1:])) - x0 |
|
346 |
+ alpha = solve(vectors.T, point - x0) |
|
344 | 347 |
if any(alpha < -eps) or sum(alpha) > 1 + eps: |
345 | 348 |
return [] |
346 | 349 |
|
... | ... |
@@ -403,7 +406,7 @@ class Triangulation: |
403 | 406 |
# we do not really need the center, we only need a point that is |
404 | 407 |
# guaranteed to lie strictly within the hull |
405 | 408 |
hull_points = self.get_vertices(self.hull) |
406 |
- pt_center = np.average(hull_points, axis=0) |
|
409 |
+ pt_center = average(hull_points, axis=0) |
|
407 | 410 |
|
408 | 411 |
pt_index = len(self.vertices) |
409 | 412 |
self.vertices.append(new_vertex) |
... | ... |
@@ -447,7 +450,7 @@ class Triangulation: |
447 | 450 |
tuple (center point, radius) |
448 | 451 |
The center and radius of the circumscribed circle |
449 | 452 |
""" |
450 |
- pts = np.dot(self.get_vertices(simplex), transform) |
|
453 |
+ pts = dot(self.get_vertices(simplex), transform) |
|
451 | 454 |
return circumsphere(pts) |
452 | 455 |
|
453 | 456 |
def point_in_cicumcircle(self, pt_index, simplex, transform): |
... | ... |
@@ -455,13 +458,13 @@ class Triangulation: |
455 | 458 |
eps = 1e-8 |
456 | 459 |
|
457 | 460 |
center, radius = self.circumscribed_circle(simplex, transform) |
458 |
- pt = np.dot(self.get_vertices([pt_index]), transform)[0] |
|
461 |
+ pt = dot(self.get_vertices([pt_index]), transform)[0] |
|
459 | 462 |
|
460 |
- return np.linalg.norm(center - pt) < (radius * (1 + eps)) |
|
463 |
+ return norm(center - pt) < (radius * (1 + eps)) |
|
461 | 464 |
|
462 | 465 |
@property |
463 | 466 |
def default_transform(self): |
464 |
- return np.eye(self.dim) |
|
467 |
+ return eye(self.dim) |
|
465 | 468 |
|
466 | 469 |
def bowyer_watson(self, pt_index, containing_simplex=None, transform=None): |
467 | 470 |
"""Modified Bowyer-Watson point adding algorithm. |
... | ... |
@@ -532,9 +535,9 @@ class Triangulation: |
532 | 535 |
volume is only dependent on the shape of the simplex and not on the |
533 | 536 |
absolute size. Due to the weird scaling, the only use of this method |
534 | 537 |
is to check that a simplex is almost flat.""" |
535 |
- vertices = np.array(self.get_vertices(simplex)) |
|
538 |
+ vertices = array(self.get_vertices(simplex)) |
|
536 | 539 |
vectors = vertices[1:] - vertices[0] |
537 |
- average_edge_length = np.mean(np.abs(vectors)) |
|
540 |
+ average_edge_length = mean(abs(vectors)) |
|
538 | 541 |
return self.volume(simplex) / (average_edge_length ** self.dim) |
539 | 542 |
|
540 | 543 |
def add_point(self, point, simplex=None, transform=None): |
... | ... |
@@ -587,8 +590,8 @@ class Triangulation: |
587 | 590 |
return self.bowyer_watson(pt_index, actual_simplex, transform) |
588 | 591 |
|
589 | 592 |
def volume(self, simplex): |
590 |
- prefactor = np.math.factorial(self.dim) |
|
591 |
- vertices = np.array(self.get_vertices(simplex)) |
|
593 |
+ prefactor = factorial(self.dim) |
|
594 |
+ vertices = array(self.get_vertices(simplex)) |
|
592 | 595 |
vectors = vertices[1:] - vertices[0] |
593 | 596 |
return float(abs(fast_det(vectors)) / prefactor) |
594 | 597 |
|