Optimize harmonic processing speed with numba

[anangsahroni]

Hi! Thank you very much for the code, I use it heavily on my thesis.

As mentioned in #7 which one of the todo list items is to convert the code to low level language to increase the performance, here I provide a numba-optimized code for rfpy.harmonics. numba translates Python functions to optimized machine code at runtime, it can be an alternative or temporary solution to performance optimization without rewriting the code to lower level language.

I optimized rfpy_harmonics and created a simple performance comparison by running the script for station MMPY for 9 different scenarios depending on the number of receiver functions (nrf), each of the scenario is executed based on two cases: with and without numba.

The device and environment for testing:

ROG Zephyrus G14 GA401QH_GA401QH
AMD Ryzen 7 5800HS with Radeon Graphics (16 CPUs), ~3.2GHz
8192MB RAM
Ubuntu 18.04 WSL2, Windows 11
Python 3.8.13

The input parameters and the result:

Table:

station	start	end	trange	nrf	time_without_numba	time_with_numba
MMPY	2017-01-01	2018-03-30	2-10	137	691.202392578125	47.04966640472412
MMPY	2017-01-01	2017-11-06	2-10	120	630.3968324661255	42.97214984893799
MMPY	2017-01-01	2017-09-23	2-10	100	537.779358625412	33.96166014671326
MMPY	2017-01-01	2017-08-10	2-10	80	420.6219515800476	25.681000471115112
MMPY	2017-01-01	2017-07-14	2-10	60	314.8274037837982	21.694521188735962
MMPY	2017-01-01	2017-06-15	2-10	41	91.12138342857361	8.571951150894165
MMPY	2017-01-01	2017-05-15	2-10	30	68.21855211257935	6.957925796508789
MMPY	2017-01-01	2017-04-09	2-10	21	52.25753951072693	5.868242502212524
MMPY	2017-01-01	2017-02-23	2-10	11	32.92144703865051	4.929212331771851

Graph:

As shown on the table and graph, numba can significantly increase the performance.

Changes

1. Added _dcomp_calculate_harmonics_isolated static method to isolate harmonics calculation based on numba, obspy.Stream object in this function is converted to numpy.array to accomodate numba
2. Added use_numba as argument for rfpy.harmonics.Harmonic.dcomp_find_azim, when this argument is True then the processing will be executed by numba
3. Added progressbar when processing is running using tqdm and numba_progressbar
4. Added --use-numba argument for rfpy_harmonics
5. Added numba and numba_progressbar in setup.py
6. Added documentation in rfpy.rst and tutorial.rst

Notes

1. Because of the current behavior of numba as mentioned here, canceling processing with CTRL+C is currently not supported. Actually I have managed to successfully solve this problem on branch harmonics-numba-sigint but it requires GCC (Linux) or MSVC (Windows) to work.
2. Sometimes the returned azimuth can be 180 degrees shifted compared to no numba because of the floating point behavior of numba as mentioned here and here. For example, when running the first scenario on the table (137 nrf), the variances of C1 for index 83 and 173, corresponding to azimuth 166 and 346 respectively:

• with numba: 0.0008908156080249394 and 0.0008908156080249397
• without numba: 0.0008908156080249395 and 0.0008908156080249395

3. numba will be a little bit slower in the first run because of compilation

Screenshots