[Octave-bug-tracker] [bug #60818] delaunayn - 2D code path vectorization doesn't match nD algorithm

[Nicholas Jankowski]

URL:
  <https://savannah.gnu.org/bugs/?60818>

                 Summary: delaunayn - 2D code path vectorization doesn't match
nD algorithm
                 Project: GNU Octave
            Submitted by: nrjank
            Submitted on: Thu 24 Jun 2021 02:38:13 PM EDT
                Category: Octave Function
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: None
                  Status: None
             Assigned to: None
         Originator Name: Nicholas Jankowski
        Originator Email: 
             Open/Closed: Open
                 Release: dev
         Discussion Lock: Any
        Operating System: Any

    _______________________________________________________

Details:

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)




    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?60818>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/