cwt with fft convolution support

[alsauve]

This pull request contains the support for fast convolution support with fast fourier transform (FFT) as discussed in issue #479 .
The speed improvement is visible for signals over signal/scales over 100/50, and becomes highly significant by more than a factor 10 for EEG like signals with 100K+ samples.

From a user point of view the only change is the new keyword method='fft', the behaviour of the function has been kept identical to the 1.0.2 release for both signature and dtype of output. Support for compliance between input and output dtype has been prepared.

• The patch passes the 1500+ test suite and provides additionals tests which compare the numerical results of both methods.
• documentation has been added to the cwt() function header

[alsauve]

This is my first pull request on GitHub, please tell me if everything is in order

[grlee77]

Great, thanks @alsauve! I will take a look.

The initial failures in the automated Windows testing on Appveyor are related to lack of np.float128 (and np.complex256) on that platform. Probably it is best to explicitly check for the existence of these attributes before adding them to the dictionary of test cases (e.g. using hasattr(np, 'float128')):

[alsauve]

Hi @grlee77,

For this typing issue, after some digging I found in numpy.core.numeric that a dict is provided for this kind of issue.
Example :

np.sctypes['complex']
Out[18]: [numpy.complex64, numpy.complex128, numpy.complex256]
np.sctypes['float']
Out[19]: [numpy.float16, numpy.float32, numpy.float64, numpy.float128]

Before upgrading to a more robust type handling in cwt() I suggest to comment out for the time being the test function draft test_cwt_dtype()

[grlee77]

This looks pretty good. I will push some minor stylistic changes to your branch and then leave comments on any remaining issues.

I found in numpy.core.numeric

Nice, I was not aware of that one. I would just delete test_cwt_dtype for now. We can add dtype-related tests in a separate PR addressing dtype handling.

[grlee77]

Congrats on starting your first PR!

I am in favor of the functionality added in this PR, but please see the review comments for some issues to consider.

[grlee77]

Also, my changes in a767fc6 are mostly to code that predates this PR. They do not change any functionality, but are just things I noticed and updated a bit when making the other style fixes. The refactored code should be a bit more efficient, but it won't make a big difference because the computation time is dominated by the convolutions.

[alsauve]

A merged version has been pushed with the comments taken into account.

Also the next_fast_len() in cunjunction with caching of the fft of data buffer improved significantly the performances (see figures in change request), so it is included in this version.

[alsauve]

The appveyor build has failed due to the new requirement on module scipy for next_fast_len
This can be an issue if the goal for pywt is to keep a dependency list as small as possible.
Otherwise the required function is short an can be included as it depends only on standard python functions.

[grlee77]

Yeah, we don't want to require scipy, but could make it an optional dependency and just fall back to the simple power of two solution when it is unavailable. (try/except ImportError)

In the future it may be beneficial to allow use of the fft from scipy rather than numpy because we have plans on the scipy end to allow user configuration of the fft library used under the hood. Hopefully we will have a student working on that this summer as part of Google Summer of Code, but that has not been finalized yet. I would just keep using np.fft.fft for now.

[grlee77]

It used to be that NumPy and SciPy were both using FFTPack to do the FFTs. Very recently (for the upcoming 0.17 release), NumPy switched to pocketfft instead. Based on its readme:

Efficient codelets are available for the factors:

2, 3, 4, 5, 7, 11 for complex-valued FFTs
2, 3, 4, 5 for real-valued FFTs

The SciPy function uses factors 2, 3 and 5 so it will still be applicable (although neglecting potential factors of 7 and 11).

[grlee77]

minor typos:

numbe -> number
"less efficient as" -> "less efficient than"

[alsauve]

Hi, thanks for the typos, and btw, congrats for having taken your project to google summer camp.
I looked for pywavelet on the google camp site but could not find it.

In which section is it?

[grlee77]

The GSOC project I was referring to is for SciPy itself rather than PyWavelets. The various Python projects planning to participate in GSOC are listed at:
https://blogs.python-gsoc.org/en/project-ideas/

and the specific project I mentioned is the fftpack one described here:
https://github.com/scipy/scipy/wiki/GSoC-2019-project-ideas

At this point SciPy has received final proposals from potential students, but we won't know until May 6th whether our proposal has been accepted by Google.

[alsauve]

BTW I made plots of a benchmark reference with a jupiter notebook here:
CWT method benchmark

The first run was somewhat biased by a cron job so I'll rerun rerunning the bench on my laptop. But the general idea is there, the graph helps understanding where the method change become interresting.

[grlee77]

BTW I made plots of a benchmark reference with a jupiter notebook here:
CWT method benchmark

Great. I will try to run that with base numpy (nomkl) and then one with (mkl). I think both variants are now available via conda-forge.

We also have a benchmark suite for some things in PyWavelets, although we have not set up any mechanism to host the results online (I think NumPy does that for theirs). If you look in the benchmarks subfolder there is a README.rst that describes it a bit.

There is a basic benchmark for cwt here:
https://github.com/PyWavelets/pywt/blob/master/benchmarks/benchmarks/cwt_benchmarks.py

We could add a case to also loop over methods. If you want to try that we could add it in this PR or we can leave it to a future PR. If you do so, you have to test the benchmarks locally to make sure they run. They do not get automatically run by the CI services.

[grlee77]

Hi @alsauve: sorry this PR was neglected for some time. I think it looks pretty good. I updated the import logic a bit and added a cases running the FFT-based variant to the CWT benchmarks.

Benchmark results (with SciPy 1.3, so using scipy.fftpack for the FFTs):

As expected, cases with larger size and larger number of scales show large benefit for the FFT-based approach.

via asv run --python=same --quick in the benchmarks folder:

[  2.63%] ··· cwt_benchmarks.CwtTimeSuite.time_cwt                                                                                                                      ok
[  2.63%] ··· ====== ========= =========== ========== =========== ========== ============ ===========
              --                                         max_scale / method                          
              ---------------- ----------------------------------------------------------------------
                n     wavelet   16 / conv   16 / fft   64 / conv   64 / fft   256 / conv   256 / fft 
              ====== ========= =========== ========== =========== ========== ============ ===========
                32      cmor     6.22±0ms   5.93±0ms    9.79±0ms   17.0±0ms    55.7±0ms     79.2±0ms 
                32     cgau4     10.3±0ms   12.3±0ms    19.3±0ms   28.1±0ms    74.0±0ms     102±0ms  
                32      fbsp     8.94±0ms   12.6±0ms    25.4±0ms   43.2±0ms    140±0ms      211±0ms  
                32     gaus4     6.77±0ms   9.17±0ms    15.7±0ms   26.5±0ms    62.6±0ms     97.6±0ms 
                32      mexh     6.41±0ms   7.90±0ms    16.2±0ms   26.8±0ms    76.3±0ms     117±0ms  
                32      morl     6.54±0ms   8.97±0ms    16.4±0ms   27.2±0ms    76.1±0ms     119±0ms  
                32      shan     9.11±0ms   12.5±0ms    25.5±0ms   42.5±0ms    145±0ms      207±0ms  
               128      cmor     9.02±0ms   11.5±0ms    23.5±0ms   30.2±0ms    131±0ms      131±0ms  
               128     cgau4     10.6±0ms   12.5±0ms    22.0±0ms   29.1±0ms    97.0±0ms     102±0ms  
               128      fbsp     9.58±0ms   13.2±0ms    33.2±0ms   45.9±0ms    217±0ms      218±0ms  
               128     gaus4     6.78±0ms   9.43±0ms    17.6±0ms   25.7±0ms    80.3±0ms     98.6±0ms 
               128      mexh     6.37±0ms   9.89±0ms    18.5±0ms   28.7±0ms    99.8±0ms     121±0ms  
               128      morl     6.42±0ms   9.01±0ms    18.8±0ms   28.1±0ms    78.0±0ms     90.8±0ms 
               128      shan     7.93±0ms   10.2±0ms    25.9±0ms   36.6±0ms    197±0ms      190±0ms  
               512      cmor     8.93±0ms   9.34±0ms    30.3±0ms   26.6±0ms    219±0ms      111±0ms  
               512     cgau4     8.76±0ms   11.2±0ms    24.5±0ms   25.0±0ms    148±0ms      90.2±0ms 
               512      fbsp     10.5±0ms   11.6±0ms    52.4±0ms   36.0±0ms    466±0ms      198±0ms  
               512     gaus4     6.65±0ms   7.95±0ms    24.0±0ms   28.9±0ms    128±0ms      106±0ms  
               512      mexh     7.70±0ms   10.2±0ms    27.4±0ms   32.6±0ms    156±0ms      130±0ms  
               512      morl     7.33±0ms   10.9±0ms    27.7±0ms   31.8±0ms    187±0ms      87.3±0ms 
               512      shan     7.31±0ms   8.12±0ms    37.0±0ms   27.0±0ms    486±0ms      179±0ms  
               2048     cmor     8.36±0ms   8.16±0ms    56.5±0ms   25.9±0ms    949±0ms      122±0ms  
               2048    cgau4     8.47±0ms   8.61±0ms    67.2±0ms   43.5±0ms    600±0ms      138±0ms  
               2048     fbsp     23.6±0ms   20.9±0ms    186±0ms    59.3±0ms    2.19±0s      253±0ms  
               2048    gaus4     5.97±0ms   6.86±0ms    25.1±0ms   25.0±0ms    274±0ms      92.4±0ms 
               2048     mexh     10.6±0ms   11.3±0ms    42.3±0ms   35.8±0ms    423±0ms      149±0ms  
               2048     morl     11.1±0ms   13.2±0ms    56.4±0ms   43.7±0ms    438±0ms      111±0ms  
               2048     shan     12.8±0ms   10.6±0ms    186±0ms    36.7±0ms    1.90±0s      254±0ms  
              ====== ========= =========== ========== =========== ========== ============ ===========

[grlee77]

I will probably have to rebase on current master to get the CI to pass here.

[alsauve]

Hi @grlee77, good to hear about you. I was myself pretty busy and also neglected to restart this thread.

The times of your benchmark seem in line with my own checks here.
And the performance inversion is clearly visible for large scales.

I remember now that the last addition of scipy next_fast_len() had improved significantly the processing times. I consider this path pretty stable and usable as is, so you can proceed for the merge :-)

[grlee77]

@rgommers: Is there any objection to an optional import of scipy here? We fall back to numpy when it is not available so we do not need to add SciPy as a dependency.

I guess importing scipy may add 100-200 ms to the import time of PyWavelets, but I don't see much other drawback to allowing this and I think the potential to use the new functionality in scipy.fft makes it worth it.

@alsauve: The rationale for preferring scipy.fft when available is that this new SciPy module is based on pypocketfft rather than FFTPack and has a number of performance and accuracy improvements. It also has mechansims for the user to offload the FFT computations to third party libraries that define an appropriate backend.

[rgommers]

@rgommers: Is there any objection to an optional import of scipy here? We fall back to numpy when it is not available so we do not need to add SciPy as a dependency.

I agree, useful speedup. Would be good to add SciPy in one of the CI runs though.

[grlee77]

closes #479

[grlee77]

Thanks @alsauve!