[WIP] Numba accelerated alpha filtrations

[ulupo]

A ridiculously preliminary version of a numba-accelerated version of alpha.py. The numba dependency has not yet been added to the installation requirements.

There are two families of changes (besides some minor code cleanup):

1. Changes due to numba's stricter requirements on types. One of the effects of this is that some class methods have to become separate functions. Furthermore, the warnings in the circumcenter function currently break numba compilation so I commented them out.
2. Miscellaneous algorithmic changes. Chiefly, I removed some redundant conditionals and the else clause in the for loop over the sigmas. I also comment out some code which was unused, e.g. the variable assignment currently in line 198. And the commented out code in lines 112-127 are an alternative "global" version of the Gabriel condition check which I believe can shave off some further runtime due to fewer calls to the circumcenter function (as mentioned briefly to @ctralie).

[sauln]

Sweet! I'd love to see some numbers on how this improves the speed! I assume it's big time.

[ulupo]

Indeed there is a speedup, at least anecdotally. Sorry I don't have time right now to give you a comprehensive benchmark test, but here are some very quick numbers from my MacBook:

• 100 pts in 2D: for X = 2 * tadasets.dsphere(n=100, d=1, r=1, noise=0.2), from 41.1 ms to 16.5 ms;
• 1000 pts in 3D: for X = 2 * tadasets.dsphere(n=1000, d=2, r=1, noise=0.2), from 7.68 s to 2.2 s or even 1.76 s when checking the Gabriel condition globally (i.e. using the commented-out alternative for faces).

And this is without properly coding things up for maximum use of that sweet numba sauce.

[ulupo]

Update. The latest incarnation improves upon the version on master as follows on my machine (noise is set to 0.2 in all experiments):

4D

• 10000 pts: from ~21min 5s to ~1min 50s.
• 1000 pts: from ~1min 45s ± 137ms to 8.88s ± 105ms.
• 100 pts: from 6.18s ± 56.5ms to 523ms ± 11.6ms.

This is a ~11.8x linear speedup.

3D

• 10000 pts: from ~1min 26s to ~11.6s.
• 1000 pts: from 8.13s ± 137ms to 1.1s ± 16.6ms.
• 100 pts: from 660ms ± 21.2ms to 88.3ms ± 2.23ms.

This is a ~7.5x linear speedup.

2D

• 10000 pts: from 4.71s ± 37.5ms to 1.18s ± 34ms.
• 1000 pts: from 464ms ± 7.28ms to 114ms ± 3.13ms.
• 100 pts: from 44.9ms ± 2ms to 11.1ms ± 204µs.

This is a ~4x linear speedup.

The extra performance comes mainly from the following changes:

1. New numba-accelerated code for the computation of the circumradius and/or the circumcentre, leveraging Miroslav Fiedler's method and using a jitted implementation of the pairwise distance function (considerably faster than scipy counterparts such as pdist + squareform or cdist, for small matrices).
2. Avoiding duplicate computations: all simplices in a non-top dimension dim can be stored as dictionary keys while traversing the set of simplices of dimension dim + 1. Henceforth, they can be efficiently enumerated later on without repetition. Currently, circumcircle calculations are duplicated because a simplex typically has several cofaces.

One notable missing feature is any sort of treatment of numerical issues coming from lack of general position. This is still WIP with @ctralie. And correctness even in the presence of general position still needs to be checked more thoroughly!

Another frontier is jitting not just the workhorse circumcircle functions but also the main algorithm. There is a technical issue here in that numba does not currently accept variable-length tuples in nopython mode, but I might be able to find a workaround.

[codecov]

Codecov Report

Merging #23 (765f724) into master (eb20753) will decrease coverage by 18.12%.
The diff coverage is 27.67%.

❗ Current head 765f724 differs from pull request most recent head 3d8b9ef. Consider uploading reports for the commit 3d8b9ef to get more accurate results

@@             Coverage Diff             @@
##           master      #23       +/-   ##
===========================================
- Coverage   94.47%   76.34%   -18.13%     
===========================================
  Files          12       12               
  Lines         398      427       +29     
===========================================
- Hits          376      326       -50     
- Misses         22      101       +79     

Impacted Files	Coverage Δ
cechmate/filtrations/alpha.py	37.87% <27.67%> (-59.21%)	⬇️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update eb20753...3d8b9ef. Read the comment docs.

[ulupo]

I was finally able to njit everything! As an example of the latest figures, I can now run 10000 points in 4 dimensions in about 20 s (from the original 21 mins, i.e. a 60x speedup!). The main bottleneck is in fact a "stupid" list comprehension doing little more than dictionary unpacking and happening outside of the Numba code. At the moment, this can take up to almost 50% of the whole execution time! <- Solved!