# Delaunayn - trivial triangle removal criteria

and then 4D just adds a p5, and using one level of laplace expansion is a sum of 4 dot.cross’s with the right signs. Could continue this way to dim-n, with expected diminishing returns and continually increasing number of terms (as n! I think).

of course it’s possible that dot and cross have enough overhead that just explicitly writing out the a11a22a33 + ... rather than dot(p, cross(...)) would save a significant time. I’ve done this before and it is generally much faster than a for loop with det, and the n-dimensional vectorization can produce a huge speed increase. but it’s easiest to write with smaller determinants.

at some point there’s a memory issue though. the for loop is slow but low memory. vectorizations with large datasets can hit upper memory bounds fairly easily if we’re sloppy about it, especially as n grows.

Is there any established way of using faster methods until input size is ‘too large’ at which point you switch to a looped code path?

I’m not sure I understand what you are trying to do.
But something like this might work for the 3-d case (using the triple product you are proposing):

diff -r a48151f59b69 scripts/geometry/delaunayn.m
--- a/scripts/geometry/delaunayn.m	Mon Jun 28 10:36:25 2021 +0200
+++ b/scripts/geometry/delaunayn.m	Mon Jun 28 16:04:27 2021 +0200
@@ -108,6 +108,12 @@
det = cross (p12, p23, 2);
idx = abs (det (:,3) ./ sqrt (sumsq (p12, 2))) < tol & ...
abs (det (:,3) ./ sqrt (sumsq (p23, 2))) < tol;
+  elseif (nd == 4)
+    ## 3-D case
+    np = rows (pts);
+    X = reshape (pts(T(:,1:end-1).',:) - pts(T(:,2:end).',:), 3, nt, 3);
+    det_X = abs (dot (X(:,:,1), cross (X(:,:,2), X(:,:,3), 1), 1));
+    idx = det_X ./ sqrt (sum (sumsq (X, 1), 3)) < sqrt(3) * tol;
else
## FIXME: Vectorize this for loop or convert delaunayn to .oct function
idx = [];

Absolutely untested still.

I’m not exactly sure. But the 2-D case might not find all “slim” triangles.

diff -r a48151f59b69 scripts/geometry/delaunayn.m
--- a/scripts/geometry/delaunayn.m	Mon Jun 28 10:36:25 2021 +0200
+++ b/scripts/geometry/delaunayn.m	Mon Jun 28 16:37:28 2021 +0200
@@ -91,29 +91,38 @@
## (reference http://en.wikipedia.org/wiki/Simplex).  Any simplex with a
## relative volume less than some arbitrary criteria is rejected.  The
## criteria we use is the volume of the simplex corresponding to an
-  ## orthogonal simplex is equal edge length all equal to the edge length of
-  ## the original simplex.  If the relative volume is 1e3*eps then the simplex
-  ## is rejected.  Note division of the two volumes means that the factor
-  ## factorial(ndim+1) is dropped.
+  ## orthogonal simplex with (ndim-1) edge lengths equal to the edge lengths of
+  ## the original simplex.  If the relative volume is smaller than 1e3*eps, the
+  ## simplex is rejected.  Note division of the two volumes means that the
+  ## factor factorial(ndim+1) is dropped.
[nt, nd] = size (T);
if (nd == 3)
## 2-D case
np = rows (pts);
+    ## expand to 3-D
ptsz = [pts, zeros(np, 1)];
p1 = ptsz(T(:,1), :);
p2 = ptsz(T(:,2), :);
p3 = ptsz(T(:,3), :);
p12 = p1 - p2;
p23 = p2 - p3;
-    det = cross (p12, p23, 2);
-    idx = abs (det (:,3) ./ sqrt (sumsq (p12, 2))) < tol & ...
-          abs (det (:,3) ./ sqrt (sumsq (p23, 2))) < tol;
+    ## Use triple product to calculate volume
+    area = cross (p12, p23, 2)(:,3);  # "3rd edge" is [0,0,1]
+    idx = abs (area ./ (sqrt (sumsq (p12, 2)) .* sqrt (sumsq (p23, 2))/2)) < tol;
+  elseif (nd == 4)
+    ## 3-D case
+    np = rows (pts);
+    ptsz = [pts, zeros(np, 1)];
+    X = reshape (pts(T(:,1:end-1).',:) - pts(T(:,2:end).',:), 3, nt, 3);
+    det_X = abs (dot (X(:,:,1), cross (X(:,:,2), X(:,:,3), 1), 1));
+    idx = det_X ./ sqrt (sum (sumsq (X, 1), 3)) < tol;
else
## FIXME: Vectorize this for loop or convert delaunayn to .oct function
idx = [];
+    # tol = sqrt (nd-1) * tol;
for i = 1:nt
X = pts(T(i,1:end-1),:) - pts(T(i,2:end),:);
-      if (abs (det (X)) / sqrt (sumsq (X, 2)) < tol)
+      if (abs (det (X)) / sqrt (sum (sumsq (X, 2), 1)) < tol)
idx(end+1) = i;
endif
endfor

haven’t read through your version yet, but I should have put my full code block in. The following has 2d, 3d, 4d cases broken out but not optimized, and it’s written out for readability, not memory use. I did get rid of the ptsz data duplication for the 2D case, put all of the p = … p12 = …, etc., data duplication is still in there. I had just repeated the <tol checks for 3d and 4d, which knowing isn’t ‘right’ began this whole post. (also not sure if I follow where your sqrt(3)*tol came from.)

...
[nt, nd] = size (T);
##  if (nd == 3)
if (any (nd == [3, 4, 5]))
## 2,3,4D cases
p1 = pts(T(:,1), :);
p2 = pts(T(:,2), :);
p3 = pts(T(:,3), :);
p12 = p1 - p2;
p23 = p2 - p3;

if (nd == 3)
# 2D
p12 = [p12, zeros(nt, 1)];
p23 = [p23, zeros(nt, 1)];
tri_areas = abs (cross (p12, p23, 2));
idx = (tri_areas(:,3) ./ sqrt (sumsq (p12, 2))) < tol & ...
(tri_areas(:,3) ./ sqrt (sumsq (p23, 2))) < tol;

else
## 3D or 4D, calculate additional vectors, do higher lev volume calc.
p4 = pts(T(:,4), :);
p34 = p3 - p4;
if (nd == 4)
## 3D
tet_vols = abs (dot (p12, cross (p23, p34, 2), 2));
idx = (tet_vols ./ sqrt (sumsq (p12, 2))) < tol & ...
(tet_vols ./ sqrt (sumsq (p23, 2))) < tol & ...
(tet_vols ./ sqrt (sumsq (p34, 2))) < tol;

else
## 4D
p5 = pts(T(:,5), :);
p45 = p4 - p5;
pentachoron_vols = ...  ## 1 level laplace expansion
p12(:,1).* dot (p23(:,2:4), cross (p34(:,2:4), p45(:,2:4), 2), 2) ...
- p12(:,2).* dot (p23(:,[1,3:4]), cross (p34(:,[1,3:4]), p45(:,[1,3:4]), 2), 2) ...
+ p12(:,3).* dot (p23(:,[1:2,4]), cross (p34(:,[1:2,4]), p45(:,[1:2,4]), 2), 2) ...
- p12(:,4).* dot (p23(:,1:3), cross (p34(:,1:3), p45(:,1:3), 2), 2);

idx = (pentachoron_vols./ sqrt (sumsq (p12, 2))) < tol & ...
(pentachoron_vols ./ sqrt (sumsq (p23, 2))) < tol & ...
(pentachoron_vols ./ sqrt (sumsq (p34, 2))) < tol & ...
(pentachoron_vols ./ sqrt (sumsq (p45, 2))) < tol;

endif
endif

else
keyboard
## FIXME: Vectorize this for loop or convert delaunayn to .oct function
idx = [];
for i = 1:nt
X = pts(T(i,1:end-1),:) - pts(T(i,2:end),:);
if (abs (det (X)) / sqrt (sumsq (X, 2)) < tol)
idx(end+1) = i;
endif
endfor
endif
T(idx,:) = [];

endfunction

I had started looking at more generalized vectorization such as you did with the X = reshape(... part. Because what you did for 3D would work for nD with the right variables for reshape. Then I got stuck on generalizing the determinant for ndim>3. Using 1 level of laplace expansion isn’t too cumbersome to just write out 4D as above. 5D would have 5x as many statements as that, etc, etc. (even more efficient full cofactor expansion goes as O(n!*n). Looking at reshaping so individuals volumes projected in 3rd dimension,
I looked at using LU decomposition instead (O(n^3), since if A = LU, abs(det(A)) = abs( diag(U)) (or something along those lines). I never got back to figuring out if there’s a decent, parallel way of calculating a whole list of U’s.

timing test cases:

rand("state", [1:625]');
x = rand(1000,1)*4-2; y = rand(1000,1)*4-2; z = rand(1000,1)*4-2; p = rand(1000,1)*4-2;  q = rand(1000,1)*4-2;

everything pushed to the for loop (even 2D)

>> tic;for idx = 1:10, delaunayn([x,y]);endfor;toc
Elapsed time is 1.29099 seconds.
>> tic;for idx = 1:10, delaunayn([x,y,z]);endfor;toc
Elapsed time is 4.1812 seconds.
>> tic;for idx = 1:10, delaunayn([x,y,z,p]);endfor;toc
Elapsed time is 18.694 seconds.
>> tic;for idx = 1:10, delaunayn([x,y,z,p,q]);endfor;toc
Elapsed time is 90.6858 seconds.

my script above (only vectorizing up to 4D)

>> tic;for idx = 1:10, delaunayn([x,y]);endfor;toc
Elapsed time is 0.063765 seconds.
>> tic;for idx = 1:10, delaunayn([x,y,z]);endfor;toc
Elapsed time is 0.206687 seconds.
>> tic;for idx = 1:10, delaunayn([x,y,z,p]);endfor;toc
Elapsed time is 1.55906 seconds.

comparing improved 3D’s:

yours:
>> tic;for idx = 1:1000, delaunayn([x,y,z]);endfor;toc
Elapsed time is 18.942 seconds.

mine:
Ȁ>> tic;for idx = 1:1000, delaunayn([x,y,z]);endfor;toc
Elapsed time is 18.3565 seconds.

so (absent the questionable validity of the all of the different >tol checks) both scripts seems to make a huge difference on timing.

playing a little more. I turned the streamlined 3D determinant into a subfunction, then made 1-level recursive laplace expansion functions for up to 7D determinants. See attached.
delaunayn.m (9.9 KB)

timing results compared to current function:

dim, n, before, after, imp.factor
2D, 1000, 4.91613, 4.89016, 1.03x
3D, 100, 42.1211, 1.8576, 22.675x
4D, 10, 18.2257, 1.34021, 13.599x
5D, 10, 91.6353, 15.0109, 6.1046x
6D, 1, 49.0621, 25.0106, 1.9617x
7D, 1, 284.357, 159.944, 1.7779x

I thought we’d see diminished returns by 6. but guess not. If i substracted the actual delaunay calculation time and just timed the removal algorithm it would look even better. Its really just a copy/paste job and some quick repetitive editing to go to higher dims, but i think 7D is more than enough.

If making the vectors was changed to something like your X() and reshaped, it could probably be done in one common command for all nDims, and then the nested if’s could probably be replaced with a switch (nd).

Similarly, with a common X the <tol check (I didn’t bother to ‘fix’ the ones in there right now) the code can probably be made more general and performed after the volume calculations.

Neither would change speed but might clean up the code a bit.

just now seeing an interesting general approach on the bug tracker that works vectorized for nDim. Uses the LU decomposition after restructuring the problem as a sparse block-diagonal matrix so the linear algebra can be done in parallel.