cleanup

Dockerfile

@@ -10,9 +10,8 @@ RUN apt-get update && apt-get install -y \
     # Testing dependencies
     python3-pytest-cov jupyter \
     # Optional NIFTy dependencies
-    libfftw3-dev python3-mpi4py python3-matplotlib python3-pynfft \
+    python3-mpi4py python3-matplotlib python3-pynfft \
   # more optional NIFTy dependencies
-  && pip3 install pyfftw \
   && pip3 install git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git \
   && pip3 install git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git \
   && pip3 install git+https://gitlab.mpcdf.mpg.de/mtr/pypocketfft.git \

README.md

@@ -47,9 +47,9 @@ Installation
 
 - [Python 3](https://www.python.org/) (3.5.x or later)
 - [SciPy](https://www.scipy.org/)
+- [pypocketfft](https://gitlab.mpcdf.mpg.de/mtr/pypocketfft)
 
 Optional dependencies:
-- [pyFFTW](https://pypi.python.org/pypi/pyFFTW) for faster Fourier transforms
 - [pyHealpix](https://gitlab.mpcdf.mpg.de/ift/pyHealpix) (for harmonic
     transforms involving domains on the sphere)
 - [nifty_gridder](https://gitlab.mpcdf.mpg.de/ift/nifty_gridder) (for radio
@@ -73,28 +73,12 @@ NIFTy5 and its mandatory dependencies can be installed via:
 
     sudo apt-get install git python3 python3-pip python3-dev
     pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/nifty.git@NIFTy_5
+    pip3 install --user git+https://gitlab.mpcdf.mpg.de/mtr/pypocketfft
 
 Plotting support is added via:
 
     sudo apt-get install python3-matplotlib
 
-NIFTy uses Numpy's FFT implementation by default. For large problems FFTW may be
-used because of its higher performance. It can be installed via:
-
-    sudo apt-get install libfftw3-dev
-    pip3 install --user pyfftw
-
-To enable FFTW usage in NIFTy, call
-
-    nifty5.fft.enable_fftw()
-
-at the beginning of your code.
-
-(Note: If you encounter problems related to `pyFFTW`, make sure that you are
-using a pip-installed `pyFFTW` package. Unfortunately, some distributions are
-shipping an incorrectly configured `pyFFTW` package, which does not cooperate
-with the installed `FFTW3` libraries.)
-
 Support for spherical harmonic transforms is added via:
 
     pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git

demos/bench_gridder.py

@@ -5,7 +5,6 @@ import numpy as np
 
 import nifty5 as ift
 
-ift.fft.enable_fftw()
 np.random.seed(40)
 
 N0s, a0s, b0s, c0s = [], [], [], []

docs/source/installation.rst

@@ -9,28 +9,12 @@ NIFTy5 and its mandatory dependencies can be installed via::
 
     sudo apt-get install git python3 python3-pip python3-dev
     pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/nifty.git@NIFTy_5
+    pip3 install --user git+https://gitlab.mpcdf.mpg.de/mtr/pypocketfft
 
 Plotting support is added via::
 
     sudo apt-get install python3-matplotlib
 
-NIFTy uses Numpy's FFT implementation by default. For large problems FFTW may be
-used because of its higher performance. It can be installed via::
-
-    sudo apt-get install libfftw3-dev
-    pip3 install --user pyfftw
-
-To enable FFTW usage in NIFTy, call::
-
-    nifty5.fft.enable_fftw()
-
-at the beginning of your code.
-
-(Note: If you encounter problems related to `pyFFTW`, make sure that you are
-using a pip-installed `pyFFTW` package. Unfortunately, some distributions are
-shipping an incorrectly configured `pyFFTW` package, which does not cooperate
-with the installed `FFTW3` libraries.)
-
 Support for spherical harmonic transforms is added via::
 
     pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git

nifty5/fft.py

@@ -20,38 +20,6 @@ import numpy as np
 import pypocketfft
 
 
-_use_fftw = False
-_fftw_prepped = False
-_fft_extra_args = {}
-
-
-def enable_fftw():
-    global _use_fftw
-    _use_fftw = True
-
-
-def disable_fftw():
-    global _use_fftw
-    _use_fftw = False
-
-
-def _init_pyfftw():
-    global _fft_extra_args, _fftw_prepped
-    if not _fftw_prepped:
-        import pyfftw
-        from pyfftw.interfaces.numpy_fft import fftn, rfftn, ifftn
-        pyfftw.interfaces.cache.enable()
-        pyfftw.interfaces.cache.set_keepalive_time(1000.)
-        # Optional extra arguments for the FFT calls
-        # if exact reproducibility is needed,
-        # set "planner_effort" to "FFTW_ESTIMATE"
-        import os
-        nthreads = int(os.getenv("OMP_NUM_THREADS", "1"))
-        _fft_extra_args = dict(planner_effort='FFTW_ESTIMATE',
-                               threads=nthreads)
-        _fftw_prepped = True
-
-
 # FIXME this should not be necessary ... no one should call a complex FFT
 # with a float array.
 def _make_complex(a):
@@ -65,77 +33,24 @@ def _make_complex(a):
 
 
 def fftn(a, axes=None):
-    if _use_fftw:
-        from pyfftw.interfaces.numpy_fft import fftn
-        _init_pyfftw()
-        return fftn(a, axes=axes, **_fft_extra_args)
-    else:
-        return pypocketfft.fftn(_make_complex(a), axes=axes)
+    return pypocketfft.fftn(_make_complex(a), axes=axes)
 
 
 def rfftn(a, axes=None):
-    if _use_fftw:
-        from pyfftw.interfaces.numpy_fft import rfftn
-        _init_pyfftw()
-        return rfftn(a, axes=axes, **_fft_extra_args)
-    else:
-        return pypocketfft.rfftn(a, axes=axes)
+    return pypocketfft.rfftn(a, axes=axes)
 
 
 def ifftn(a, axes=None):
-    if _use_fftw:
-        from pyfftw.interfaces.numpy_fft import ifftn
-        _init_pyfftw()
-        return ifftn(a, axes=axes, **_fft_extra_args)
+    # FIXME this is a temporary fix and can be done more elegantly
+    if axes is None:
+        fct = 1./a.size
     else:
-        # FIXME this is a temporary fix and can be done more elegantly
-        if axes is None:
-            fct = 1./a.size
-        else:
-            fct = 1./np.prod(np.take(a.shape, axes))
-        return pypocketfft.ifftn(_make_complex(a), axes=axes, fct=fct)
+        fct = 1./np.prod(np.take(a.shape, axes))
+    return pypocketfft.ifftn(_make_complex(a), axes=axes, fct=fct)
 
 
 def hartley(a, axes=None):
-    # Check if the axes provided are valid given the shape
-    if axes is not None and \
-            not all(axis < len(a.shape) for axis in axes):
-        raise ValueError("Provided axes do not match array shape")
-    if iscomplextype(a.dtype):
-        raise TypeError("Hartley transform requires real-valued arrays.")
-
-    tmp = rfftn(a, axes=axes)
-
-    def _fill_array(tmp, res, axes):
-        if axes is None:
-            axes = tuple(range(tmp.ndim))
-        lastaxis = axes[-1]
-        ntmplast = tmp.shape[lastaxis]
-        slice1 = (slice(None),)*lastaxis + (slice(0, ntmplast),)
-        np.add(tmp.real, tmp.imag, out=res[slice1])
-
-        def _fill_upper_half(tmp, res, axes):
-            lastaxis = axes[-1]
-            nlast = res.shape[lastaxis]
-            ntmplast = tmp.shape[lastaxis]
-            nrem = nlast - ntmplast
-            slice1 = [slice(None)]*lastaxis + [slice(ntmplast, None)]
-            slice2 = [slice(None)]*lastaxis + [slice(nrem, 0, -1)]
-            for i in axes[:-1]:
-                slice1[i] = slice(1, None)
-                slice2[i] = slice(None, 0, -1)
-            slice1 = tuple(slice1)
-            slice2 = tuple(slice2)
-            np.subtract(tmp[slice2].real, tmp[slice2].imag, out=res[slice1])
-            for i, ax in enumerate(axes[:-1]):
-                dim1 = (slice(None),)*ax + (slice(0, 1),)
-                axes2 = axes[:i] + axes[i+1:]
-                _fill_upper_half(tmp[dim1], res[dim1], axes2)
-
-        _fill_upper_half(tmp, res, axes)
-        return res
-
-    return _fill_array(tmp, np.empty_like(a), axes)
+    return pypocketfft.hartley2(a, axes=axes)
 
 
 # Do a real-to-complex forward FFT and return the _full_ output array

test/test_operators/test_fft_operator.py

@@ -34,22 +34,9 @@ def _get_rtol(tp):
 pmp = pytest.mark.parametrize
 dtype = list2fixture([np.float64, np.float32, np.complex64, np.complex128])
 op = list2fixture([ift.HartleyOperator, ift.FFTOperator])
-fftw = list2fixture([False, True])
-
-
-def test_switch():
-    ift.fft.enable_fftw()
-    assert_(ift.fft._use_fftw is True)
-    ift.fft.disable_fftw()
-    assert_(ift.fft._use_fftw is False)
-    ift.fft.enable_fftw()
-    assert_(ift.fft._use_fftw is True)
-
 
 @pmp('d', [0.1, 1, 3.7])
-def test_fft1D(d, dtype, op, fftw):
-    if fftw:
-        ift.fft.enable_fftw()
+def test_fft1D(d, dtype, op):
     dim1 = 16
     tol = _get_rtol(dtype)
     a = ift.RGSpace(dim1, distances=d)
@@ -69,16 +56,13 @@ def test_fft1D(d, dtype, op, fftw):
         domain=a, random_type='normal', std=7, mean=3, dtype=dtype)
     out = fft.inverse_times(fft.times(inp))
     assert_allclose(inp.local_data, out.local_data, rtol=tol, atol=tol)
-    ift.fft.disable_fftw()
 
 
 @pmp('dim1', [12, 15])
 @pmp('dim2', [9, 12])
 @pmp('d1', [0.1, 1, 3.7])
 @pmp('d2', [0.4, 1, 2.7])
-def test_fft2D(dim1, dim2, d1, d2, dtype, op, fftw):
-    if fftw:
-        ift.fft.enable_fftw()
+def test_fft2D(dim1, dim2, d1, d2, dtype, op):
     tol = _get_rtol(dtype)
     a = ift.RGSpace([dim1, dim2], distances=[d1, d2])
     b = ift.RGSpace(
@@ -97,13 +81,10 @@ def test_fft2D(dim1, dim2, d1, d2, dtype, op, fftw):
         domain=a, random_type='normal', std=7, mean=3, dtype=dtype)
     out = fft.inverse_times(fft.times(inp))
     assert_allclose(inp.local_data, out.local_data, rtol=tol, atol=tol)
-    ift.fft.disable_fftw()
 
 
 @pmp('index', [0, 1, 2])
-def test_composed_fft(index, dtype, op, fftw):
-    if fftw:
-        ift.fft.enable_fftw()
+def test_composed_fft(index, dtype, op):
     tol = _get_rtol(dtype)
     a = [a1, a2,
          a3] = [ift.RGSpace((32,)),
@@ -115,7 +96,6 @@ def test_composed_fft(index, dtype, op, fftw):
         domain=(a1, a2, a3), random_type='normal', std=7, mean=3, dtype=dtype)
     out = fft.inverse_times(fft.times(inp))
     assert_allclose(inp.local_data, out.local_data, rtol=tol, atol=tol)
-    ift.fft.disable_fftw()
 
 
 @pmp('space', [
@@ -123,9 +103,7 @@ def test_composed_fft(index, dtype, op, fftw):
     ift.RGSpace((15, 27), distances=(.7, .33), harmonic=True),
     ift.RGSpace(73, distances=0.5643)
 ])
-def test_normalisation(space, dtype, op, fftw):
-    if fftw:
-        ift.fft.enable_fftw()
+def test_normalisation(space, dtype, op):
     tol = 10*_get_rtol(dtype)
     cospace = space.get_default_codomain()
     fft = op(space, cospace)
@@ -138,4 +116,3 @@ def test_normalisation(space, dtype, op, fftw):
     assert_allclose(
         inp.to_global_data()[zero_idx], out.integrate(), rtol=tol, atol=tol)
     assert_allclose(out.local_data, out2.local_data, rtol=tol, atol=tol)
-    ift.fft.disable_fftw()