introduce 'fast_det' for 2x2 and 3x3 matrices

Author: basnijholt

adaptive/learner/learnerND.py

@@ -13,18 +13,19 @@
 
 from ..notebook_integration import ensure_holoviews, ensure_plotly
 from .triangulation import (Triangulation, point_in_simplex,
-                            circumsphere, simplex_volume_in_embedding)
+                            circumsphere, simplex_volume_in_embedding,
+                            fast_det)
 from ..utils import restore, cache_latest
 
 
 def volume(simplex, ys=None):
     # Notice the parameter ys is there so you can use this volume method as
     # as loss function
-    matrix = np.array(np.subtract(simplex[:-1], simplex[-1]), dtype=float)
-    dim = len(simplex) - 1
+    matrix = np.subtract(simplex[:-1], simplex[-1], dtype=float)
 
     # See https://www.jstor.org/stable/2315353
-    vol = np.abs(np.linalg.det(matrix)) / np.math.factorial(dim)
+    dim = len(simplex) - 1
+    vol = np.abs(fast_det(matrix)) / np.math.factorial(dim)
     return vol
 
 

adaptive/learner/triangulation.py

@@ -126,6 +126,17 @@ def fast_3d_circumcircle(points):
     return center, radius
 
 
+def fast_det(matrix):
+    matrix = np.asarray(matrix, dtype=float)
+    if matrix.shape == (2, 2):
+        return matrix[0][0] * matrix[1][1] - matrix[1][0] * matrix[0][1]
+    elif matrix.shape == (3, 3):
+        a, b, c, d, e, f, g, h, i = matrix.ravel()
+        return a * (e*i - f*h) - b * (d*i - f*g) + c * (d*h - e*g)
+    else:
+        return np.linalg.det(matrix)
+
+
 def circumsphere(pts):
     dim = len(pts) - 1
     if dim == 2:
@@ -140,9 +151,9 @@ def circumsphere(pts):
     for i in range(1, len(pts)):
         r = np.delete(mat, i, 1)
         factor = (-1) ** (i + 1)
-        center.append(factor * np.linalg.det(r))
+        center.append(factor * fast_det(r))
 
-    a = np.linalg.det(np.delete(mat, 0, 1))
+    a = fast_det(np.delete(mat, 0, 1))
     center = [x / (2 * a) for x in center]
 
     x0 = pts[0]
@@ -227,7 +238,7 @@ def simplex_volume_in_embedding(vertices) -> float:
     sq_dists_mat = scipy.spatial.distance.squareform(bordered)
 
     coeff = - (-2) ** (num_verts-1) * factorial(num_verts-1) ** 2
-    vol_square = np.linalg.det(sq_dists_mat) / coeff
+    vol_square = fast_det(sq_dists_mat) / coeff
 
     if vol_square < 0:
         if abs(vol_square) < 1e-15:
@@ -587,7 +598,7 @@ def volume(self, simplex):
         prefactor = np.math.factorial(self.dim)
         vertices = np.array(self.get_vertices(simplex))
         vectors = vertices[1:] - vertices[0]
-        return float(abs(np.linalg.det(vectors)) / prefactor)
+        return float(abs(fast_det(vectors)) / prefactor)
 
     def volumes(self):
         return [self.volume(sim) for sim in self.simplices]