ENH: Add pocketfft sources to numpy for testing, benchmarks, etc.

[mreineck]

This pull request was opened at the request of @matthew-brett and is related to issues #11885, #11886 and #10625.

[ncgeib]

I ran the same scripts as I did for #10625 to compare the accuracy and speed. The accuracy issue for real-valued FFTs is fixed and in general pocketfft is faster than the current implementation in numpy but still slower than scipy's fftpack (not taking the additional twiddle factor creation into account). For prime sizes it is vastly faster due to using an algorithm with a lower complexity.

All seems good to me. Only thing I could think of is that maybe the implementation of pocketfft could be sped up further.

EDIT: Removed speed benchmarks as they did not reflect the current status of the PR. See below.

[mreineck]

Thanks a lot for these tests, @ncgeib!

Concerning performance, I don't fully understand yet why scipy's implementation is faster for regular sizes ... the algorithms are basically the same.

Hmm, one difference is that pocketfft does its own memory management, i.e. for every transform it calls malloc() and free() a handful of times. This could cause a slowdown, especially for short transforms. Maybe that's what we are seeing here.

[mreineck]

Sorry, I'm being stupid ... of course I know why pocketfft is currently slower than scipy!

At the moment, pocketfft is not creating plans beforehand and caching them for re-use. That means that for every repetition of the transform it recomputes the twiddle factors all over again.

This can be fixed of course, but I'll probably need help writing a Python class that caches the pocketfft plans.

[matthew-brett]

I wonder whether some small Cython classes would help here ...

[ncgeib]

@mreineck Yeah, sorry if that was not clear. The comparison is not really fair as the timings exclude the twiddle factor calculation and cache creation for numpy and scipy and include the twiddle factor calculation for pocketfft.
Still, I am not sure if this explains the whole runtime difference. numpy's C port of fftpack is also slower than the Fortran version. We would have to measure either the twiddle computation or use a cached interface of pocketfft.

[mreineck]

The code now does some caching.

If someone could tell me the canonical way to pass a C pointer to Python, I'd be grateful. Currently I'm casting its value to a "long" and feel really dirty about it.

[charris]

canonical way to pass a C pointer to Python

NpyCapsule_*, see numpy/core/include/numpy/npy_3kcompat.h and grep for usage.

[rgommers]

I ran the same scripts as I did for #10625 to compare the accuracy and speed.

@ncgeib nice. It may be useful to contribute the accuracy benchmarks as tests and the speed ones as asv benchmarks to the various packages, if you have time. I doubt our current tests and benchmarks are as comprehensive.

[ncgeib]

I re-ran the speed tests. The comparison should now be fair (and more or less accurate). We see that pocketfft is still just a little bit slower than scipy's fftpack for most test cases. Maybe there is some overhead that could be reduced?

@rgommers If I find the time I will do that. The benchmarks basically copy the ones FFTW does. For the meantime, they can be found at https://github.com/ncgeib/benchmark-ffts

[rgommers]

@rgommers If I find the time I will do that. The benchmarks basically copy the ones FFTW does. For the meantime, they can be found at https://github.com/ncgeib/benchmark-ffts

oh okay, if you mean an actual copy then never mind - FFTW is GPL licensed so we cannot use that code. standalone repo is still useful though, we can then run it again in case of other PRs with major changes like this one

[ncgeib]

@rgommers Sorry, I meant they copy the methodology the FFTW authors use. No code was copied in any way - I simply found the way they test the accuracy of the FFT implementations to be good (comparing against a high-precision calculation of the DFT). This allows to compare the graphs to the ones on their webpage.

I specifically wrote the benchmark scripts while hoping that at one point they could be used when numpy finally gains a better FFT implementation. In that sense I am quite happy now ;)

[mreineck]

@ncgeib I feel much more comfortable with the new benchmarks! The remaining performance difference between scipy and pocketfft is largest for small transform sizes and almost vanishes for longer ones, so my hope is that we can still gain a little bit by optimizing the wrapper code on the Python side.

[mreineck]

Hmmm ... I don't see any locking/unlocking going on in the scipy code when the caches are updated. Could this be the reason for the performance difference? And why is scipy not subject to race conditions?

[rgommers]

Hmmm ... I don't see any locking/unlocking going on in the scipy code when the caches are updated. Could this be the reason for the performance difference? And why is scipy not subject to race conditions?

From scipy/scipy#1419: ... This means that SciPy's FFTs are not threadsafe, and thus the GIL cannot be released. Thus, SciPy has to lock up the interpreter while doing FFTs. ...

[mreineck]

Maybe there's our performance overhead.

Also, scipy is doing cache management on the C side, which could also help.

But anyway I think this is a detail we can look into once we are convinced that everything else works as intended.

[mreineck]

There seem to be unnecessary array copies in fftpack.py (and via copy/paste also in pocketfft.py):

The function fft starts with:

    a = asarray(a).astype(complex, copy=False)

(i.e. potentially there is no copy), while ifft starts with:

    # The copy may be required for multithreading.
    a = array(a, copy=True, dtype=complex)

This looks conspicuously inconsistent...

[mreineck]

BTW, I think that most of numpy.fft's slowness compared to scipy is explained by the many calls to _check_size() in helper.py. This is seriously suboptimal code.

[mreineck]

After some more measurements I'm convinced that the remaining overhead is caused by cache lookup and management. If I use a primitive cache which stores just a single transform, I get speeds that lie practically on top of the scipy data points everywhere.
Not sure what to make of that...

[mreineck]

This is what I get with the latest patch. Please keep in mind that this is not directly comparable to earlier performance plots, because they were measured on a different computer!

[ncgeib]

I can confirm these on my machine. In the latest state of the PR the runtime of pocketfft is basically the same as scipy's fftpack for regular numbers.

[mreineck]

Great!
So, what's the procedure from here on?

All Travis builds fail (even those for Python 3), because C89 is enforced. Do we need to request a change in policy here?

[charris]

We will try to move to C99 after Python 2.7 is dropped. Note that other parts of NumPy will need some minor fixes, or at least that was my impression last time I tried with the C99 flag.

[charris]

I'm thinking of branching 1.16.0 around the end of November.

[charris]

You can try changing the build/test environment in this PR to see what needs doing, but it would be best to keep those changes as separate commits as we might want to split the resulting PR when we move up. In particular, you can drop the 2.7 tests and look to making the rest of NumPy C99 safe. If some of that can be done without breaking 2.7 on windows, so much the better.

[charris]

I've put up #11905 to start testing with c99 on travisCI.

[charris]

Looks like we need to modify lines 344-352 in runtests.py to remove -Werror=declaration-after-statement and add -std=c99. I'll make a PR.

[charris]

See #12610.

[mhvk]

This all looks great, but does need a release note!

Also, having not looked in any detail: the old fft could not do single precision; can pocketfft? Important for large but noisy data sets. (But obviously no reason to hold up the PR; if it were possible to implement easily, however, maybe raise an issue so it can be done later.)

[charris]

@mreineck I've made a separate PR for the travis-test.sh testing fix, after that goes in you will want to rebase, probably as part of adding a release note.

[tacaswell]

I think matplotlib did that way back when and had some problems, athough I'm not really familiar with that package. @tacaswell Are you folks still using C++.

Yes, we still have a bunch of c++ (mostly because Agg is templated c++).

[charris]

@mreineck Might also want to squash the 84 commits :)

[mreineck]

@mreineck I've made a separate PR for the travis-test.sh testing fix, after that goes in you will want to rebase, probably as part of adding a release note.

What would that look like? Some addition in doc/changelog/1.17.0-changelog.rst?

I'll fetch the latest changes from master and squash as soon as the tester runs through.

[mreineck]

Also, having not looked in any detail: the old fft could not do single precision; can pocketfft?

Not at the moment. I think it could be added, but I don't really have the experience where we have to keep double precision to avoid losing too much accuracy. If memory usage is the concern, people will be better off with an in-place implementation that computes twiddle factors on the fly. Performance-wise I don't expect any speedup potential unless explicit vectorization support is added.

[charris]

@mreineck Put the notes doc/release/1.17.0-notes.rst under "Improvements". I think a one line list entry in "Highlights" would also be appropriate, i.e., * Blah, blah....

[mhvk]

@mreineck - with single precision, all calculations would be single precision; in the old fftpack, there was

#define DOUBLE
#ifdef DOUBLE
#define Treal double
#else
#define Treal float
#endif

and then all code is defined in terms of Treal. I think if one were to add this, one would write the code duplication using the machinery also used for ufuncs. But that is for another PR! (Perhaps as part of making the fft routines generalized ufuncs...)

[mreineck]

Simply replacing double by float will have very undesirable consequences for twiddle factor computation. But I agree, that is for another time...

[mreineck]

OK, I think I'm done.

[charris]

OK, let's get this in. I haven't reviewed the sources, but I think the FFT is not difficult to test, so will rely on the tests and user feedback leading up to the 1.17.0 release. Thanks @mreineck .

[charris]

@mreineck We should also add some benchmarks for this. There are currently some in scypy that can probably be modified to cover things like large(ish) prime sizes.

[mreineck]

@charris Thanks a lot for merging!

I agree that correctness testing for FFTs is not hard to do. The previously existing unit tests have been kept, and I added identity tests (real and complex) for all lengths between 1 and 512. That should catch most potential problems.

Concerning benchmarks, I think @ncgeib's work (#11888 (comment)) would be the best start.