realized i let this sit for some time. I've refreshed the patch, made a few minor tweaks and format corrections, added some input validation and BISTS, and updated NEWS.md.8. passes all tests and a make check, so doesn't seem to break the few other functions that call it.

again, this takes care of the initial bug, and low hanging fruit for performance improvement.  since it's more than just a bugfix, pushed to default for v8 as https://hg.savannah.gnu.org/hgweb/octave/rev/bf8f33249e86

will retitle and point bug #53942 here for future work on >3D path performance and memory use improvement.

Marking as Ready for Test.

minor patch update, fixing an error introduced by some of the simplex check indexing when _delaunayn_ returns T = 0.  now just bypasses the simplex checking in that case.  no other significant changes yet.

(file #51877)

patch of latest version attached.

I had already made some changes to drop memory use a bit. used ~ on lu outputs for all but u and p.  then 'reused' eqs instead of u, got rid the one reordered tri index variable. i think it dropped mem usage by ~10%, but the big memory jumps are still when forming eqs and running lu.



(file #51753)

We could arbitrarily split up the vector equations by simplexes, but in what way? half? quarter? When we set out I was wondering if there was a way to loop over the dimensions. With volume calculation being the main task I figured some aspect of iteration increasing the volume projection would work.  not something I could see with the LU decomp though.

I'm still not sure why the method does so poorly with low point count, but other than that, with the tweaked version that removes calls to 'dot' and 'cross' which is fairly quick anyway, the laplace methed is faster with lower memory up to and including 6D.  So if LU is just used for degree > 6, maybe it's natural to partition by degree number.

A good question. In the extreme each volume’s eqn  diagonal could be independently calculated, that’s what the for loop did. And I’m pretty sure  the blocks would independently diagonalize (since you can use matrix row operations instead of lu decomp to get the diagonal, and row ops to get row echelon form are block independent), but not sure how to partition that in any non-arbitrary way.

I’m half convinced this whole problem stems from not making use of what’s already going on in the initial Delaunay calculation, and we’re doubling effort for an after the fact check.

so all that said:

biggest memory jump seems to happen inside lu.  no real way to work around that.

laplace does scale worse than lu at high dimension and high number of points, but seems laplace still edges it out to about 5D (ignoring hte small 10pt case), after which laplace scaling gets pretty poor.

just based on timing, I'd recommend using laplace expansion up to 5D, then switch to the LU approach.  unfortunately that will run into memory issues for >5D. I don't really know any way around that for those running larger, higher dimension triangulations. Personally, I don't know anyone running such high dimension triangulations.

Attached is the m file which would do that. Haven't worked up a patch yet, as this should still get some BISTs added to it.



(file #51724)

ok. ran a bunch of speed and memory checks for different dimensions (up to 7D) and  numbers of points(10-10000). There were definite points at high dimension where the LU approach was faster than the direct laplace expansion.  However, the memory footprint is 3-6x higher.

looking at the code to see where the memory increases come in, for one example, 4D and 10k points - peak mem LU-919MB, laplace-233MB :

after running main _delaunayn_: working mem 66MB, peak 224MB
after building eqs:  working 185MB, peak 587MB
after running LU:  working 250MB, peak 919MB

So, both building eqs and running LU are significant memory users, doubling and quadrupling the other method's usage. While my machine hit memory limits just running the _delaunayn_ for 7D @ 10k pts, i also ran out of mem with the LU algorithm at 6D, 10k pts.

timing:
I ran _delaunayn_ through the tests without the simplex volume check to get a baseline, then ran the old loop and the two vectorized approaches.  did enough loops and tic/tocs to get a good average, then figured out average time/operation and subtracted out _delaunayn_'s time to compare just the checks.

here was the original loop (i didn't take time getting the really low memory numbers):

old loop
est max memory
dim        10        100        1000        10000
2d                                55
3d                        50        71
4d                        57        153
5d                        98        655
6d                        322        3500
7d                        1572        OOM

old loop
check time, ms
dim  10        100        1000      10000
2d   0.412164  0.462299   0.44162   0.9745
3d   1.341657  38.14379   435.1904  4755.529
4d   1.309238  128.04121  1950.719  22246.64
5d   1.371815  426.4956   8552.84   121692.96
6d   0.861909  1427.137   46292.5   667997
7d   0.646422  4765.347   221883.5  OOM

the LU algorithm:

LU
est max memory
dim        10        100        1000        10000
2d                                64
3d                        66        252
4d                54        123        933
5d                63        528        5880
6d                140        3100        OOM
7d                453        13000        OOM

lu
check time, ms
dim  10      100     1000     10000
2d   0.4971  0.9715  6.6968   69.4045
3d   0.4556  2.3249  36.732   636.099
4d   0.5227  15.078  336.72   4516.34
5d   0.5603  96.882  2586.3   35837.96
6d   0.5480  494.23   18004   OOM
7d   0.5369  2471.3  205843   OOM

lu
check/loop
dim  10           100   1000   10000
2d   1.20  2.10  15.16  71.22
3d   0.33  0.06  0.084  0.133
4d   0.39  0.11  0.172  0.203
5d   0.40  0.22  0.302  0.294
6d   0.63  0.34  0.388        OOM
7d   0.83  0.51  0.927        OOM

and the laplace expansion:

laplace
est max memory
dim        10        100        1000        10000
2d                                55
3d                                78
4d                                233
5d                        85        1100
6d                        590        8400
7d                105        3800        OOM

laplace
check time, ms
dim  10             100     1000      10000
2d   0.2019  0.2276  0.10983   9.0344
3d   0.3230  0.3180  1.2734    33.169
4d   0.9225  1.1493  28.029    533.24
5d   3.9620  13.853  418.83    9204.96
6d   23.342  298.25  18801.6   278337
7d   164.25  7893.5  656063.5  OOM

laplace
check/loop
dim  10     100    1000   10000
2d   0.489  0.490  0.248  -0.28
3d   0.240  0.008  0.002  0.006
4d   0.704  0.008  0.014  0.023
5d   2.888  0.032  0.048  0.074
6d   27.08  0.208  0.406  0.413
7d   254.0  1.656  2.956  OOM

ok, so comparing the two:

memory_lu/memory_laplace
dim  10  100  1000  10000
2d                  1.16
3d                  3.23
4d                  4.00
5d            6.21  5.34
6d            5.25
7d            4.31  3.42

time_lu / time_laplace
dim  10     100   1000   10000
2d   2.462  4.26  60.97  7.68
3d   1.410  7.30  28.84  19.1
4d   0.566  13.1  12.01  8.46
5d   0.141  6.99  6.175  3.89
6d   0.023  1.65  0.957  inf
7d   0.003  0.31  0.313  OOM

@nrjank: Did you come around to running the performance tests and adding some self tests?

and now looking at the geometric alternative of doing the volume comparison for the LU code, and I'm realizing we may be calculating volume wrong in the Laplace expansion code. the simple orthogonal cases may have just masked it.

volume of a triangle, tetragonal, etc., is 1/n! * volume of the parallelogram/parallelepiped, etc., defined by the vector product. but that requires use of a coterminated vectors.  (p21, p31, p41, etc.) not (p21, p32, p43) like we have there.  I was just following the 2D code which used p12 and p23, since they're co-terminated except for a sign. But using p12, p23, p34, ... for higher dims uses vectors not part of the parallelepiped to calculate it's volume, which i suspect is very wrong.

It's a trivial switch in the Laplace expansion determinant code.  will make that change and update the patch after check on the effect of leaving q out or changing it to q in the LU code.

Dmitri, according to the Octave help

"If A is full then subroutines from LAPACK are used, and if A is sparse then UMFPACK is used."

https://octave.sourceforge.io/octave/function/lu.html

I'm not sure if this makes any particular difference?

I opened bug #60867 for `whos` not showing the sparse attribute.

I would have thought `eqs` was a sparse (block diagonal) matrix. But it looks like it is a full matrix in comment #26.
That could also be the reason why the LU decomposition is unexpectedly slow.
I wonder which command/assignment causes the conversion from sparse to full.

Question with the patch should

if (any(nd == [3]))

be

if (nd == 2 || nd == 3)

Regarding the variables at the end

• l and q are never used
• u could be overwritten rather than using a new variable R
• eqs is not needed after the input to lu
• edgvec is not used after it enters eqs but if the check is changed this may also change

I added "whos" at the end of the function. This is what I see for 6D 1000 pt case:

octave:2> whos
Variables visible from the current scope:

variables in scope: top scope

   Attr Name        Size                     Bytes  Class
   ==== ====        ====                     =====  =====
        p        1000x1                       8000  double
        q        1000x1                       8000  double
        r        1000x1                       8000  double
        s        1000x1                       8000  double
        x        1000x1                       8000  double
        y        1000x1                       8000  double
        z        1000x1                       8000  double

Total is 7000 elements using 56000 bytes

octave:3> tic; delaunayn([x,y,z,p,q,r]); toc
Variables visible from the current scope:

variables in scope: delaunayn: /home/dima/scratch/delaunayn.m

   Attr Name                 Size                     Bytes  Class
   ==== ====                 ====                     =====  =====
        R              3654654x1                   58474480  double
    f   T               608426x7                   34071856  double
        edgvec         3654654x6                  175423392  double
        eqs            3654654x3654654            380084024  double
        idx                  1x683                     5464  double
        l              3654654x3654654            233897864  double
        nd                   1x1                          8  double
        nt                   1x1                          8  double
        p              3654654x1                   29237232  double
    f   pts               1000x6                      48000  double
        q              3654654x1                   29237232  double
        reorderdtriidx       1x3654654             29237232  double
        tol                  1x1                          8  double
        u              3654654x3654654            233897864  double
    f   varargin             0x0                          0  cell

Total is 40069528391356 elements using 1203614664 bytes

Elapsed time is 6.89901 seconds.

Do we actually need "l" matrix? I do not see it being used
anywhere in the code.

Dmitri.
--

To fix the ordering issue change the line

reorderdtriidx = kron (1:nt, ones (1,nd-1))(p)(q);

to

reorderdtriidx = kron (1:nt, ones (1,nd-1))(p);

Please check this I normally use lu on full matrixes so not completely across the two outputs of sparse matrices. This change fixed the example below.

Could we work out which line is causing the memory issue. I visually cannot see why it is so high. I thought the line with logical may be making a full array but calling logical on a sparse array keeps it sparse.

The reason for

(nd - 1) * nt

is the same reason as (the above the number of rows or columns of the matrix)

max (size (A))

in the rank check. I set up the test by working out which rows are causing the matrix to be singular rather than working from a geometry perspective. It is normally done using singular value decomposition or a permuted QR factorization so I am probably stretching the applicability of the check with a LU decomposition. Feel free to change the check to what you think is better.

stepping through the nD code, using a 3D example, it seems to be missing one of the zero volume tets:

Matlab:

A = [0 0 0;
1 0 0;
0 1 0;
0 0 1;
1 1 0;
1 0 1;
0 1 1;
1 1 1;
0.5 0.5 0.5];

B = delaunay(A)
B =
     3     4     1     9
     3     7     9     5
     8     9     7     5
     8     6     9     5
     8     6     7     9
     9     7     4     6
     3     7     4     9
     3     9     2     5
     9     6     2     5
     9     4     2     6
     9     4     1     2
     3     9     1     2

size(B)
ans =
    12     4

Octave old code:

>> A = [0 0 0; 1 0 0; 0 1 0; 0 0 1; 1 1 0; 1 0 1; 0 1 1; 1 1 1; 0.5 0.5 0.5]
A =

        0        0        0
   1.0000        0        0
        0   1.0000        0
        0        0   1.0000
   1.0000   1.0000        0
   1.0000        0   1.0000
        0   1.0000   1.0000
   1.0000   1.0000   1.0000
   0.5000   0.5000   0.5000

>> B = delaunayn(A)
B =

   5   9   2   1
   5   9   3   1
   6   9   2   1
   6   9   4   1
   6   5   9   2
   6   5   8   9
   7   9   3   1
   7   9   4   1
   7   5   9   3
   7   5   8   9
   7   6   9   4
   7   6   8   9

>> size(B)
ans =

   12    4

new code, using LU decomp for 3D:

>> B = delaunayn(A)
B =

   5   3   2   1
   5   9   3   1
   6   4   2   1
   6   9   2   1
   6   9   4   1
   6   5   8   2
   6   5   9   2
   7   9   3   1
   7   9   4   1
   7   5   9   3
   7   5   8   9
   7   6   9   4
   7   6   8   9

>> size(B)
ans =

   13    4

>> p12 = A(B(:,1),:)-A(B(:,2),:); p23 = A(B(:,2),:)-A(B(:,3),:); p34 = A(B(:,3),:)-A(B
(:,4),:);
>> vol = dot (p12, cross (p23, p34, 2), 2)
vol =

        0
  -0.5000
        0
  -0.5000
   0.5000
        0
   0.5000
   0.5000
  -0.5000
  -0.5000
  -0.5000
   0.5000
   0.5000

so, it left in the three 0 vol tets but dropped 2 others.  will look through to see if it's something obvious.

looking at timing. It appears that Octave LU calls are just very computational and memory expensive:

using the profiler to look at a 3D case:

Laplace Expansion:
>> rand("state", [1:625]');
>> x = rand(100,1)*4-2;y = rand(100,1)*4-2; z = rand(100,1)*4-2;
>> profile clear; profile on; for idx = 1:10000, delaunayn([x,y,z]);endfor, profile off; profshow

   #      Function Attr     Time (s)   Time (%)        Calls
------------------------------------------------------------
   5 __delaunayn__            13.772      78.08        10000
   1     delaunayn             1.899      10.77        10000
  13         cross             1.245       7.06        10000
  12      binary -             0.088       0.50        60000
  19           cat             0.088       0.50        10000
  18     binary .*             0.087       0.50        80000
  22          sqrt             0.080       0.45        30000
  21         sumsq             0.070       0.40        30000
   6           isa             0.044       0.25        10000
  20           dot             0.035       0.20        10000
   3      binary <             0.030       0.17        50000
   9          size             0.024       0.14        30000
   7           eps             0.023       0.13        10000
   2        nargin             0.022       0.13        50001
  15         ndims             0.022       0.12        30000
  10     binary ==             0.019       0.11        30000
  16          ones             0.018       0.10        10000
  23     binary ./             0.017       0.10        10000
  11           any             0.017       0.10        10000
  24           abs             0.013       0.07        10000


using the LU decomposition:

>> clear all
>> rand("state", [1:625]');
>> x = rand(100,1)*4-2;y = rand(100,1)*4-2; z = rand(100,1)*4-2; p = rand(100,1)*4-2;
q = rand(100,1)*4-2; r = rand(100,1)*4-2; s = rand(100,1)*4-2;
>> profile clear; profile on; for idx = 1:10000, delaunayn([x,y,z]);endfor, profile of
f; profshow
   #      Function Attr     Time (s)   Time (%)        Calls
------------------------------------------------------------
  23            lu            14.207      40.10        10000
   5 __delaunayn__            12.965      36.59        10000
   1     delaunayn             4.423      12.48        10000
  27        unique             0.783       2.21        10000
  15          kron             0.753       2.13        30000
  22       logical             0.576       1.63        10000
  17         speye             0.368       1.04        10000
  16        sparse             0.218       0.62        20000
  24          diag             0.184       0.52        10000
   3      binary <             0.182       0.51        30000
  25           abs             0.095       0.27        10000
  13      binary -             0.078       0.22        70000
  12    postfix .'             0.078       0.22        20000
  26           max             0.073       0.21        10000
  21          true             0.055       0.16        30000
   7           eps             0.051       0.14        20000
   6           isa             0.038       0.11        10000
  36         zeros             0.036       0.10        10000
  30         false             0.028       0.08        20001
  14          ones             0.028       0.08        20000

so at least in 3D, for 100 points the LU approach is spending just as much time in LU as in _delaunayn_, plus ~1/2 of _delaunayn_ again on other ops within delaunayn.  The Laplace expansion avoids most other function calls except for cross, spending a total of ~1/3 of the time spent in _delaunayn_.

jumping out to 6D:

LU:

>> clear all
>> x = rand(100,1)*4-2;y = rand(100,1)*4-2; z = rand(100,1)*4-2; p = rand(100,1)*4-2;q = rand(100,1)*4-2; r = rand(100,1)*4-2;
>> profile clear; profile on; for idx = 1:100, delaunayn([x,y,z,p,q,r]);endfor, profile off; profshow
   #      Function Attr     Time (s)   Time (%)        Calls
------------------------------------------------------------
  23            lu            26.310      47.92          100
   5 __delaunayn__            12.579      22.91          100
   1     delaunayn             9.997      18.21          100
  15          kron             2.593       4.72          300
  22       logical             1.388       2.53          100
  12    postfix .'             0.857       1.56          200
...


Laplace expansion:

>> clear all
>> x = rand(100,1)*4-2;y = rand(100,1)*4-2; z = rand(100,1)*4-2; p = rand(100,1)*4-2;q = rand(100,1)*4-2; r = rand(100,1)*4-2;
>> profile clear; profile on; for idx = 1:100, delaunayn([x,y,z,p,q,r]);endfor, profile off; profshow
   #          Function Attr     Time (s)   Time (%)        Calls
----------------------------------------------------------------
   5     __delaunayn__            14.085      57.97          100
  22         binary .*             1.971       8.11        87600
  17             cross             1.690       6.95        12000
  15 delaunayn>detvec4             1.341       5.52         3000
  12          binary -             1.335       5.49        44100
  14 delaunayn>detvec5             0.910       3.75          600
  24               dot             0.869       3.57        12000
...

not much improvement.


upping to 1000pts:

3D, 1000pts

LU:
   #      Function Attr     Time (s)   Time (%)        Calls
------------------------------------------------------------
  23            lu            28.536      55.89         1000
   5 __delaunayn__            17.235      33.76         1000
   1     delaunayn             2.665       5.22         1000
  15          kron             0.811       1.59         3000
  22       logical             0.674       1.32         1000


Laplace Exp:
   #          Function Attr     Time (s)   Time (%)        Calls
----------------------------------------------------------------
   5     __delaunayn__            17.201      94.61         1000
   1         delaunayn             0.443       2.43         1000
  14             cross             0.143       0.79         1000
  12          binary -             0.081       0.44         6000



6D, 1000pts:

LU decomp:
   #      Function Attr     Time (s)   Time (%)        Calls
------------------------------------------------------------
  23            lu            97.561      40.34           10
   5 __delaunayn__            85.050      35.17           10
   1     delaunayn            35.671      14.75           10
  15          kron             9.233       3.82           30
  22       logical             5.001       2.07           10
  12    postfix .'             3.478       1.44           20

(peak mem about 3GB, total time ~241s)

Laplace expansion:
   #          Function Attr     Time (s)   Time (%)        Calls
----------------------------------------------------------------
   5     __delaunayn__            85.856      35.19           10
  22         binary .*            55.332      22.68         8820
  12          binary -            30.297      12.42         4410
  15 delaunayn>detvec4            20.128       8.25          300
  23               cat            10.416       4.27         1200
  24               dot            10.393       4.26         1200

(peak mem about 700MB, total time ~244s)

so here at the cost of 3x the memory, LU decomp breaks even at 6D.


upping to 10000 points. 

3D, 10kpts:
LU:
   #      Function Attr     Time (s)   Time (%)        Calls
------------------------------------------------------------
  23            lu            41.780      47.47          100
   5 __delaunayn__            27.882      31.68          100
   1     delaunayn            10.659      12.11          100
  15          kron             2.891       3.29          300
  22       logical             1.455       1.65          100


Laplace expansion:
   #      Function Attr     Time (s)   Time (%)        Calls
------------------------------------------------------------
   5 __delaunayn__            27.508      91.32          100
   1     delaunayn             1.558       5.17          100
  12      binary -             0.665       2.21          600
  19           cat             0.099       0.33          100
  22          sqrt             0.076       0.25          300
  18     binary .*             0.073       0.24          800

6D:

Laplace expansion:
   #          Function Attr     Time (s)   Time (%)        Calls
----------------------------------------------------------------
   5     __delaunayn__           195.219      44.54            1
  22         binary .*            87.128      19.88          876
  12          binary -            46.551      10.62          441
  15 delaunayn>detvec4            29.981       6.84           30
  24               dot            16.814       3.84          120

peak memory usage by Octave in the first few minutes hit about 3GB, then 7GB toward the end

LU:
peak memory usage by Octave in the first few minutes hit about 3GB, then 13GB toward the end, then
hit an out of memory limit.

incomplete profshow:
>> profshow
   #      Function Attr     Time (s)   Time (%)        Calls
------------------------------------------------------------
   5 __delaunayn__           198.244      54.59            1
   1     delaunayn           132.159      36.39            1
  15          kron            13.536       3.73            2
  22       logical             8.555       2.36            1
  12    postfix .'             5.262       1.45            2
  13      binary -             3.891       1.07            5
  16        sparse             1.474       0.41            2
  23       profile             0.011       0.00            1
  24     binary !=             0.002       0.00            1

I think it crashed out during the LU decomp, which is just showing up as delaunayn time.  Not sure how well that time can be trusted, since the times I was hitting 97+ memory I saw the disk usage spike as it started trying to diskswap. 

So in summary, seems like with Octave's LU it starts to be faster at higher dimensions for larger arrays, but at a big memory cost.

I think the issue mentioned in comment 15 is a bug with sparse. Simple example:

a=sparse(1);
a(1)=single(20);

To fix change the line with an issue to

eqs(logical (kron (speye (nt, nt), true (nd-1))))=double(edgvec.'(:));

Attached is an mildly corrected patch - the 2D and 3D cases didn't need /2 and /6. (forgot they were dropped with the volume/volume calculation.) and I left a few typos in both 2D and 3D versions. Otherwise it's still the same.

I'll reexamine the timing tomorrow.  When I ran the tests where I manually wrote out the explicit determinants, i still had the incorrect  (and maybe faster) tol checks in there. I'm wondering if a correct tol check would make it and the LU approach comparable in speed. I also excepted to see LU decomp do better after ~5 or 6 dimensions, but the numbers below were showing that to still be 2x faster.

the easiest question - is the small simplex test necessary? for matlab compatibility it is, as they throw out zero-volume simplexes.  with floating point error determining 'zero-volume' = less than some arbitrary distance from zero.  maybe we could make an optional flag to disable it, but it needs to be there by default.

i think one simple case needing small simplex removal.  A simple cube with a center point:

A = [0 0 0;
     1 0 0;
     0 1 0;
     0 0 1;
     1 1 0;
     1 0 1;
     0 1 1;
     1 1 1;
     0.5 0.5 0.5];

delaunayn(A), size(ans)
ans =

   5   9   2   1
   5   9   3   1
   6   9   2   1
   6   9   4   1
   6   5   9   2
   6   5   8   9
   7   9   3   1
   7   9   4   1
   7   5   9   3
   7   5   8   9
   7   6   9   4
   7   6   8   9

ans =

   12    4

the size before small simplexes were removed, however:

T =

   5   3   2   1
   5   9   2   1
   5   9   3   1
   6   4   2   1
   6   9   2   1
   6   9   4   1
   6   5   8   2
   6   5   9   2
   6   5   8   9
   7   9   3   1
   7   9   4   1
   7   5   9   3
   7   5   8   9
   7   6   9   4
   7   6   8   9

debug> size(T)
ans =

   15    4

with the volumes:

vol =

        0
   0.5000
  -0.5000
        0
  -0.5000
   0.5000
        0
   0.5000
   0.5000
   0.5000
  -0.5000
  -0.5000
  -0.5000
   0.5000
   0.5000

it apparently created three zero volume tetrahedra from the cube faces:

>> A([5 3 2 1], :)
ans =

   1   1   0
   0   1   0
   1   0   0
   0   0   0

>> A([6 4 2 1], :)
ans =

   1   0   1
   0   0   1
   1   0   0
   0   0   0

>> A([6 5 8 2], :)
ans =

   1   0   1
   1   1   0
   1   1   1
   1   0   0

(file #51627)

Theoretical basis of the check eps(max(R))
The check I proposed is a way of estimating the rank or finding the rows which cause the matrix to be singular or nearly singular. The rank function default tolerance is

tol = max (size (A)) * sigma(1) * eps;

where sigma(1) is a norm of the matrix found from the singular value decomposition. So max(R) is an estimate of the matrix norm. For this application as long as this estimate is within a few orders of magnitude of the true matrix norm I don't imagine any problems. We could use a 1, inf, Frobenius, or max norm at a small additional computational expense. I don't see the need for making the tolerance a function of every simplex, a deliberate difference in the simplex size of over 10 orders of magnitude would be required to make any difference. This method could also be applied using QR or SVD decompositions but it would be more expensive. In summary, the tolerance check I suggested is based on checking which rows are causing the matrix to be nearly singular rather than any geometrical arguments.

Suggestions

• If idx is not empty we should display a warning to the user. One of the main uses for Delaunay is mesh generation, removing simplexes could cause issues. A warning would potentially get the user to check their points and fix issues with the point locations. An example would be a nearly duplicated point and it would be better to remove the nearly duplicated point and re-mesh than to remove some simplexes.
• The comment in the latest patch about the determinant should say that it will find the determinant of one simplex all of them is more complicated, see the end of comment 1.
• Do we even need to remove zero volume simplexes, they do cause problems when solving partial differential equations on a mesh containing them but changing the location of the points typically allows the Delaunay algorithm to find a mesh without low volume simplexes.

Can we check the time data with what we should expect:

• LU: I expect nt*(nd-1)^3 plus some overhead
• Leibniz rule repeated 2x2 det: nt*((nd-1)!)

So around six dimensions, the LU method should become faster and it should become much faster for higher dimensions.

attached is a version of the patch that removes some unneeded assignments from 2D (no real need for ptsz and separate p1,2,etc), and adds in a similar path for 3d.  the 3D tol check matches 2D calculating the volume of an equivalent right-angled tetrahedron, and tweaked a few comments.

see attached.

(file #51625)

That is a pity. It could probably have saved a lot of work, had we found bug #53942 before. Sorry for not noticing earlier.
There are probably a lot of other good patches on the tracker that went unnoticed for years.

With respect to the patch here. I haven't tested it as of now.
But I am wondering why it is possible to check all simplices against the same `eps(max(R))`. The previous algorithm used a reference that was different for each simplex.
But maybe I am missing something obvious.

The patch did apply for me. I amended it with a few comments and fixed what I think was a minor inaccuracy in the 2-D case. See the attachment.
Is that all OK?

We might want to tweak the introductory comment a bit more once I better understand the part about `eps(max(R))`.

If the number of simplices is very high, would it make sense to break the calculation into subsets of these simplices and use the new algorithm for each subset? Or is this overkill?


(file #51624)

for the record, I thought I had searched for delaunay bugs/patches when I started this.  apologies for missing bug #53942.

As requested a patch file and m file. The last few patches I have sent have not worked. I completely removed my local copy and redownloaded the latest version and hope it works this time.

(file #51621, file #51622)

hah,  that's great.  I had started looking at LU decomp for just that reason but stopped short of setting up the block diagonal matrix to make it all work in parallel.

will look at it and compare with some of the more brute force stuff I played with just breaking out the determinants (see the  Discourse discussion https://octave.discourse.group/t/delaunayn-trivial-triangle-removal-criteria, or the attached modified delaunayn.m)

we were still unsure of the validity of the arbitrary relative volume check and I was looking for something in the literature to justify one approach or another.

(file #51620)

As discussed in some depth on the Octave Discourse [1],  delaunayn.m has a longstanding FIXME to vectorize the check for trivially small simplexes.  It currently loops over every shape. bug# 53689 separated out at least the 2D case and vectorized it, showing significant speedup.

Looking at doing the same for at least 3D if not nD, we noticed the 2D codepath does not actually reproduce the algorithm used by the for loop. Vectorizing higher dimensions requires choosing what algorithm to implement, and whether the 2D case needs to be made to match >2D.

'Brief' summary of the differences:
The original code admits the evaluation criteria is arbitrary.  It compares each volume to 'a reference volume', and discards volumes that are too small. the 'reference volume', or equivalent, is the difference.

The looped nD code takes the shape-defining edge vectors, orthogonalizes them to create a resultant vector, and a matrix division of the volume by that resultant vector produces another vector. (in 2D this new vector * the first vector has the same volume as the simplex). the components of the new vector are compared to 'tol' (=1000*eps), and if all components are smaller than tol, that shape is deemed trivially small and discarded.

2D path:  simplex volume is divided by the length of each triangle edge producing a resultant length (as if it was a rectangle), and the length of that result is compared to tol. if all such vectors are <tol, that volume is discarded. 

so the main difference is that 'tol' is compared either to the calculated vector length (2d), or it's component lengths (nD). a test pushing a 2D example through both code paths shows very different comparisons are made to 'tol'. 

the benefit of the nD case is that for any dimension, tol is compared to 1D lengths. if the 2D approach were applied to nD cases, you would compare tol to lengths, area, volume,... etc., with increasing dimension.

to the 2D code path's advantage, using rdivide instead of mrdivide, is much easier to vectorize and could be extended to >2d fairly easily.

I'd recommend trying to stick with the nD case for dimensional consistency, unless a better argument can be made, although this will make even the 2D case much harder to vectorize.  That said, the code even admits the check is arbitrarily decided, and maybe we can come up with a more self-consistent 'arbitrary' measure.

[1] https://octave.discourse.group/t/delaunayn-trivial-triangle-removal-criteria

(added discourse discussion users and bug# 53689 author to notification)