ENH: Add new scipy.fft module using pocketfft

[peterbell10]

Closes #10175

Enhancements

• Support for long double transforms
• Multi-threading using OpenMP (disabled for now)
• NumPy-like real to complex rffts including n-dimensional transforms rfftn
• Orthonormal transforms (norm = 'ortho')
• Bluestein's algorithm
• GIL is released during transforms

Issues

• The build definitions don't select the extra_compile_args based on the compiler which would be needed to compile with OpenMP in a portable way.
• FFTPACK optimized for the case when real sequences were passed to fft which is missing from pocketfft currently.
• Normalization factors are only passed to c++ with double precision, even for long double transforms.
• This adds a build time dependency on pybind11 and I'm not sure if it is required at run-time as well.
• Docstrings need to be written
• The pocketfft backend is passing all of the numpy.fft tests and most of the scipy.fftpack tests after some modifications. I have 1 remaining failure locally which is a precision difference for float32 transforms. Needs investigation.

Open questions

• How portable is OpenMP? (scipy.fft: Can OpenMP be used portably #10239)
• How should it be reported if multi-threading isn't available. @rgommers suggested we should raise when the user specifies a number of threads that isn't available. (Response when number of workers cannot be satisfied #10242)
• NumPy's hermitian transforms (hfft, ihfft) are missing. This is trivial to add since it is only a slight variant on rfft and irfft. We could also extend beyond numpy and give n-dimensional versions hfftn and ihfftn. (scipy.fft: What functions should be included #10240)
• Some people find the argument overwrite_x to be confusing (The overwrite_x argument of FFT functions is deceptive in certain conditions cupy/cupy#1988). It has been suggested to remove it entirely but could also be renamed. FFTW calls this FFTW_DESTROY_INPUT which makes it clearer that this is not a guarunteed inplace transform/out-parameter. (fft(pack): The overwrite_x parameter can be confusing #10241)

[rgommers]

Thanks @peterbell10, great start!

I created a new scipy.fft label for the new submodule.

For all the things that need significant discussion and deciding, I suggest to create separate issues. And keep the list in this PR as the "master issue". Otherwise things will get unwieldy. Suggestion of separate issues:

• OpenMP portability and build config (incl pybind11).
• what functions to include in the new module (e.g. hfft)
• multithreading behavior
• overwrite_x

[peterbell10]

Thanks @rgommers, that's a good idea.

[peterbell10]

Regarding the failing test, single precision fftpack tests were marked as xfail so fftpack failed for the same test case. Investigating it shows that pocketfft is in fact more accurate than fftpack.

Max abs diff fft:      2.99796178049074e-06
Max abs diff fftpack:  3.06302563664315e-06
Max rel diff fft:      4.73642558242986e-05
Max rel diff fftpack:  2.63469637785951e-04

I've given single precision variants slightly more tolerance so the test should be passing now.

[peterbell10]

I've had a look through the pybind11 repo and the only python code is either tests or the get_include functionality. That, combined with the fact that there is nothing else linked in; I'm now reasonably confident that it's only a build-time dependency.

[peterbell10]

Also, I'm planning to base a lot of the documentation off the numpy fft docs as they seemed to be slightly more consistent. I hope that's not a problem.

[rgommers]

Also, I'm planning to base a lot of the documentation off the numpy fft docs as they seemed to be slightly more consistent. I hope that's not a problem.

Definitely not, please borrow whatever is helpful!

[peterbell10]

I'm struggling to get the documentation to build correctly. The scipy.fft module docstring gives me a numpydoc error, despite it being very similar to every other module docstring

Unknown section Fast Fourier Transforms (ffts) in the docstring of <module 'scipy.fft'[...]

which comes from this section heading in the docstring, which is exactly the same as for fftpack

Fast Fourier Transforms (FFTs)
==============================

If I remove the section headings entirely, the module summary builds. However, there are no .rst files for the fft module in doc/source/generated so the functions aren't documented anywhere. I assume I need to tell the generator about the fft module somewhere but I can't see where the generator is.

Any help would be appreciated.

[rgommers]

Unknown section

Hmm, that error is normally for function docstrings. Perhaps numpydoc gets itself confused and is parsing a module docstring as function docstring. Could be because scipy.fft used to be a function and you didn't do a clean rebuild?

[peterbell10]

CircleCI should be doing a clean build and it's happening on there as well.

Edit: I just tried removing all reference to numpy.fft.fft from the scipy __init__.py and still get the same error.

[peterbell10]

It looks like the root cause is that isinstance(scipy.fft, Callable)==True, so numpydoc is treating it like a function.

[peterbell10]

I've added a workaround for sphinx documenting scipy.fft as a function and the numpydoc error is gone. However, the function documentation still isn't being generated.

I've also been trying to update the setup.py to conditionally add --std=c++11 for gcc-like compilers. The correct way to do this seems to be providing a custom build_ext class but there doesn't seem to be a standard way to update cmdclass in the numpy.distutils Configuration object.

[larsoner]

Okay the FFT people appear to be happy and no other complaints, in it goes. Thanks @peterbell10 !

[WarrenWeckesser]

This adds a build time dependency on pybind11 and I'm not sure if it is required at run-time as well.

This should be documented in INSTALL.rst.txt.

[endolith]

"discrete" typo

[endolith]

"discrete" typo

[endolith]

"multi-dimentional" typo

[endolith]

I guess these docstrings will be expanded to the same level of detail as the old ones in a different PR?

[peterbell10]

The user-facing API is in scipy/fft/_basic.py which has the full docstrings.

[peterbell10]

See the rendered docs

[peterbell10]

I've added a release note to the wiki page:
https://github.com/scipy/scipy/wiki/Release-note-entries-for-SciPy-1.4.0#scipyfft-added

[rgommers]

That looks great, thanks Peter

[WarrenWeckesser]

I was about to update some old code to use scipy.fft instead of scipy.fftpack, but I see that scipy.fft doesn't have a replacement for scipy.fftpack.diff. Was there ever a plan to add diff (or something similar) to scipy.fft? Or is something equivalent already available somewhere else?

[mreineck]

I vaguely remember that we discussed that at the time and think we agreed that no replacement was needed, since diff was not documented and not used. That said, it should be easy to add as a post-processing step, correct?

[WarrenWeckesser]

The docstring that I linked to is not complete, but I wouldn't say the function was not documented. It is definitely part of the public API. It might not be widely used, but:

• I used it years ago in the SciPy cookbook example "Korteweg de Vries equation", and at least two people have used that example as a model for their own code:

  • How to implement boundary value in Python using spectral methods
  • how to use solve_ivp solve Partial Differential Equations with spectral method?

  Presumably others have used the reference.

• It shows up now and then on stackoverflow, e.g.

  • numpy.diff and scipy.fftpack.diff giving different results when differentiating
  • Computing numeric derivative via FFT - SciPy
  • Derivative of an array in python?

• Some other random mentions on the web:

  • Practical notes on selected numerical methods with examples (PDF)
  • "SciPy for Signal Processing"

Having said all that, apparently this is the first time someone has asked about the "missing" functions, so I guess the demand isn't really that high. Certainly the method is widely used, but I suspect the existence of the function is not well known. If you understand the method, it is not hard to implement, but I wouldn't call our current implementation trivial:

scipy/scipy/fftpack/_pseudo_diffs.py

Lines 49 to 74 in 39fa485

	tmp = asarray(x)
	if order == 0:
	return tmp
	if iscomplexobj(tmp):
	return diff(tmp.real,order,period)+1j*diff(tmp.imag,order,period)
	if period is not None:
	c = 2*pi/period
	else:
	c = 1.0
	n = len(x)
	omega = _cache.get((n,order,c))
	if omega is None:
	if len(_cache) > 20:
	while _cache:
	_cache.popitem()
	def kernel(k,order=order,c=c):
	if k:
	return pow(c*k,order)
	return 0
	omega = convolve.init_convolution_kernel(n,kernel,d=order,
	zero_nyquist=1)
	_cache[(n,order,c)] = omega
	overwrite_x = _datacopied(tmp, x)
	return convolve.convolve(tmp,omega,swap_real_imag=order % 2,
	overwrite_x=overwrite_x)