summaryrefslogtreecommitdiff
path: root/silx/math/fft
diff options
context:
space:
mode:
Diffstat (limited to 'silx/math/fft')
-rw-r--r--silx/math/fft/__init__.py8
-rw-r--r--silx/math/fft/basefft.py146
-rw-r--r--silx/math/fft/clfft.py286
-rw-r--r--silx/math/fft/cufft.py253
-rw-r--r--silx/math/fft/fft.py96
-rw-r--r--silx/math/fft/fftw.py214
-rw-r--r--silx/math/fft/npfft.py124
-rw-r--r--silx/math/fft/setup.py41
-rw-r--r--silx/math/fft/test/__init__.py25
-rw-r--r--silx/math/fft/test/test_fft.py270
10 files changed, 0 insertions, 1463 deletions
diff --git a/silx/math/fft/__init__.py b/silx/math/fft/__init__.py
deleted file mode 100644
index ea12cd6..0000000
--- a/silx/math/fft/__init__.py
+++ /dev/null
@@ -1,8 +0,0 @@
-#!/usr/bin/env python
-# coding: utf-8
-
-__authors__ = ["P. Paleo"]
-__license__ = "MIT"
-__date__ = "12/12/2018"
-
-from .fft import FFT
diff --git a/silx/math/fft/basefft.py b/silx/math/fft/basefft.py
deleted file mode 100644
index 854ca37..0000000
--- a/silx/math/fft/basefft.py
+++ /dev/null
@@ -1,146 +0,0 @@
-#!/usr/bin/env python
-# coding: utf-8
-# /*##########################################################################
-#
-# Copyright (c) 2018 European Synchrotron Radiation Facility
-#
-# Permission is hereby granted, free of charge, to any person obtaining a copy
-# of this software and associated documentation files (the "Software"), to deal
-# in the Software without restriction, including without limitation the rights
-# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-# copies of the Software, and to permit persons to whom the Software is
-# furnished to do so, subject to the following conditions:
-#
-# The above copyright notice and this permission notice shall be included in
-# all copies or substantial portions of the Software.
-#
-# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
-# THE SOFTWARE.
-#
-# ###########################################################################*/
-import numpy as np
-from pkg_resources import parse_version
-
-
-def check_version(package, required_version):
- """
- Check whether a given package version is superior or equal to required_version.
- """
- try:
- ver = getattr(package, "__version__")
- except AttributeError:
- try:
- ver = getattr(package, "version")
- except Exception:
- return False
- req_v = parse_version(required_version)
- ver_v = parse_version(ver)
- return ver_v >= req_v
-
-
-class BaseFFT(object):
- """
- Base class for all FFT backends.
- """
- def __init__(self, **kwargs):
- self.__get_args(**kwargs)
-
- if self.shape is None and self.dtype is None and self.template is None:
- raise ValueError("Please provide either (shape and dtype) or template")
- if self.template is not None:
- self.shape = self.template.shape
- self.dtype = self.template.dtype
- self.user_data = self.template
- self.data_allocated = False
- self.__calc_axes()
- self.__set_dtypes()
- self.__calc_shape()
-
- def __get_args(self, **kwargs):
- expected_args = {
- "shape": None,
- "dtype": None,
- "template": None,
- "shape_out": None,
- "axes": None,
- "normalize": "rescale",
- }
- for arg_name, default_val in expected_args.items():
- if arg_name not in kwargs:
- # Base class was not instantiated properly
- raise ValueError("Please provide argument %s" % arg_name)
- setattr(self, arg_name, default_val)
- for arg_name, arg_val in kwargs.items():
- setattr(self, arg_name, arg_val)
-
- def __set_dtypes(self):
- dtypes_mapping = {
- np.dtype("float32"): np.complex64,
- np.dtype("float64"): np.complex128,
- np.dtype("complex64"): np.complex64,
- np.dtype("complex128"): np.complex128
- }
- dp = {
- np.dtype("float32"): np.float64,
- np.dtype("complex64"): np.complex128
- }
- self.dtype_in = np.dtype(self.dtype)
- if self.dtype_in not in dtypes_mapping:
- raise ValueError("Invalid input data type: got %s" %
- self.dtype_in
- )
- self.dtype_out = dtypes_mapping[self.dtype_in]
-
- def __calc_shape(self):
- # TODO allow for C2C even for real input data (?)
- if self.dtype_in in [np.float32, np.float64]:
- last_dim = self.shape[-1]//2 + 1
- # FFTW convention
- self.shape_out = self.shape[:-1] + (self.shape[-1]//2 + 1,)
- else:
- self.shape_out = self.shape
-
- def __calc_axes(self):
- default_axes = tuple(range(len(self.shape)))
- if self.axes is None:
- self.axes = default_axes
- self.user_axes = None
- else:
- self.user_axes = self.axes
- # Handle possibly negative axes
- self.axes = tuple(np.array(default_axes)[np.array(self.user_axes)])
-
- def _allocate(self, shape, dtype):
- raise ValueError("This should be implemented by back-end FFT")
-
- def set_data(self, dst, src, shape, dtype, copy=True):
- raise ValueError("This should be implemented by back-end FFT")
-
- def allocate_arrays(self):
- if not(self.data_allocated):
- self.data_in = self._allocate(self.shape, self.dtype_in)
- self.data_out = self._allocate(self.shape_out, self.dtype_out)
- self.data_allocated = True
-
- def set_input_data(self, data, copy=True):
- if data is None:
- return self.data_in
- else:
- return self.set_data(self.data_in, data, self.shape, self.dtype_in, copy=copy, name="data_in")
-
- def set_output_data(self, data, copy=True):
- if data is None:
- return self.data_out
- else:
- return self.set_data(self.data_out, data, self.shape_out, self.dtype_out, copy=copy, name="data_out")
-
- def fft(self, array, **kwargs):
- raise ValueError("This should be implemented by back-end FFT")
-
- def ifft(self, array, **kwargs):
- raise ValueError("This should be implemented by back-end FFT")
diff --git a/silx/math/fft/clfft.py b/silx/math/fft/clfft.py
deleted file mode 100644
index dad8ec1..0000000
--- a/silx/math/fft/clfft.py
+++ /dev/null
@@ -1,286 +0,0 @@
-#!/usr/bin/env python
-# coding: utf-8
-# /*##########################################################################
-#
-# Copyright (c) 2018-2019 European Synchrotron Radiation Facility
-#
-# Permission is hereby granted, free of charge, to any person obtaining a copy
-# of this software and associated documentation files (the "Software"), to deal
-# in the Software without restriction, including without limitation the rights
-# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-# copies of the Software, and to permit persons to whom the Software is
-# furnished to do so, subject to the following conditions:
-#
-# The above copyright notice and this permission notice shall be included in
-# all copies or substantial portions of the Software.
-#
-# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
-# THE SOFTWARE.
-#
-# ###########################################################################*/
-import numpy as np
-
-from .basefft import BaseFFT, check_version
-try:
- import pyopencl as cl
- import pyopencl.array as parray
- import gpyfft
- from gpyfft.fft import FFT as cl_fft
- from ...opencl.common import ocl
- __have_clfft__ = True
-except ImportError:
- __have_clfft__ = False
-
-
-# Check gpyfft version
-__required_gpyfft_version__ = "0.3.0"
-if __have_clfft__:
- __have_clfft__ = check_version(gpyfft, __required_gpyfft_version__)
-
-
-class CLFFT(BaseFFT):
- """Initialize a clfft plan.
-
- Please see FFT class for parameters help.
-
- CLFFT-specific parameters
- --------------------------
-
- :param pyopencl.Context ctx:
- If set to other than None, an existing pyopencl context is used.
- :param bool fast_math:
- If set to True, computations will be done with "fast math" mode,
- i.e., more speed but less accuracy.
- :param bool choose_best_device:
- Whether to automatically choose the best available OpenCL device.
- """
- def __init__(
- self,
- shape=None,
- dtype=None,
- template=None,
- shape_out=None,
- axes=None,
- normalize="rescale",
- ctx=None,
- fast_math=False,
- choose_best_device=True,
- ):
- if not(__have_clfft__) or not(__have_clfft__):
- raise ImportError("Please install pyopencl and gpyfft >= %s to use the OpenCL back-end" % __required_gpyfft_version__)
-
- super(CLFFT, self).__init__(
- shape=shape,
- dtype=dtype,
- template=template,
- shape_out=shape_out,
- axes=axes,
- normalize=normalize,
- )
- self.ctx = ctx
- self.choose_best_device = choose_best_device
- self.fast_math = fast_math
- self.backend = "clfft"
-
- self.fix_axes()
- self.init_context_queue()
- self.allocate_arrays()
- self.real_transform = np.isrealobj(self.data_in)
- self.compute_forward_plan()
- self.compute_inverse_plan()
- self.refs = {
- "data_in": self.data_in,
- "data_out": self.data_out,
- }
- # TODO
- # Either pyopencl ElementWiseKernel, or built-in clfft callbacks
- if self.normalize != "rescale":
- raise NotImplementedError(
- "Normalization modes other than rescale are not implemented with OpenCL backend yet."
- )
-
- def fix_axes(self):
- """
- "Fix" axes.
-
- clfft does not have the same convention as FFTW/cuda/numpy.
- """
- self.axes = self.axes[::-1]
-
- def _allocate(self, shape, dtype):
- ary = parray.empty(self.queue, shape, dtype=dtype)
- ary.fill(0)
- return ary
-
-
- def check_array(self, array, shape, dtype, copy=True):
- if array.shape != shape:
- raise ValueError("Invalid data shape: expected %s, got %s" %
- (shape, array.shape)
- )
- if array.dtype != dtype:
- raise ValueError("Invalid data type: expected %s, got %s" %
- (dtype, array.dtype)
- )
-
-
- def set_data(self, dst, src, shape, dtype, copy=True, name=None):
- """
- dst is a device array owned by the current instance
- (either self.data_in or self.data_out).
-
- copy is ignored for device<-> arrays.
- """
- self.check_array(src, shape, dtype)
- if isinstance(src, np.ndarray):
- if name == "data_out":
- # Makes little sense to provide output=numpy_array
- return dst
- if not(src.flags["C_CONTIGUOUS"]):
- src = np.ascontiguousarray(src, dtype=dtype)
- # working on underlying buffer is notably faster
- #~ dst[:] = src[:]
- evt = cl.enqueue_copy(self.queue, dst.data, src)
- evt.wait()
- elif isinstance(src, parray.Array):
- # No copy, use the data as self.d_input or self.d_output
- # (this prevents the use of in-place transforms, however).
- # We have to keep their old references.
- if name is None:
- # This should not happen
- raise ValueError("Please provide either copy=True or name != None")
- assert id(self.refs[name]) == id(dst) # DEBUG
- setattr(self, name, src)
- return src
- else:
- raise ValueError(
- "Invalid array type %s, expected numpy.ndarray or pyopencl.array" %
- type(src)
- )
- return dst
-
-
- def recover_array_references(self):
- self.data_in = self.refs["data_in"]
- self.data_out = self.refs["data_out"]
-
-
- def init_context_queue(self):
- if self.ctx is None:
- if self.choose_best_device:
- self.ctx = ocl.create_context()
- else:
- self.ctx = cl.create_some_context()
- self.queue = cl.CommandQueue(self.ctx)
-
-
- def compute_forward_plan(self):
- self.plan_forward = cl_fft(
- self.ctx,
- self.queue,
- self.data_in,
- out_array=self.data_out,
- axes=self.axes,
- fast_math=self.fast_math,
- real=self.real_transform,
- )
-
-
- def compute_inverse_plan(self):
- self.plan_inverse = cl_fft(
- self.ctx,
- self.queue,
- self.data_out,
- out_array=self.data_in,
- axes=self.axes,
- fast_math=self.fast_math,
- real=self.real_transform,
- )
-
-
- def update_forward_plan_arrays(self):
- self.plan_forward.data = self.data_in
- self.plan_forward.result = self.data_out
-
-
- def update_inverse_plan_arrays(self):
- self.plan_inverse.data = self.data_out
- self.plan_inverse.result = self.data_in
-
-
- def copy_output_if_numpy(self, dst, src):
- if isinstance(dst, parray.Array):
- return
- # working on underlying buffer is notably faster
- #~ dst[:] = src[:]
- evt = cl.enqueue_copy(self.queue, dst, src.data)
- evt.wait()
-
-
- def fft(self, array, output=None, do_async=False):
- """
- Perform a (forward) Fast Fourier Transform.
-
- :param Union[numpy.ndarray,pyopencl.array] array:
- Input data. Must be consistent with the current context.
- :param Union[numpy.ndarray,pyopencl.array] output:
- Output data. By default, output is a numpy.ndarray.
- :param bool do_async:
- Whether to perform operation in asynchronous mode.
- Default is False, meaning that we wait for transform to complete.
- """
- self.set_input_data(array, copy=False)
- self.set_output_data(output, copy=False)
- self.update_forward_plan_arrays()
- event, = self.plan_forward.enqueue()
- if not(do_async):
- event.wait()
- if output is not None:
- self.copy_output_if_numpy(output, self.data_out)
- res = output
- else:
- res = self.data_out.get()
- self.recover_array_references()
- return res
-
-
- def ifft(self, array, output=None, do_async=False):
- """
- Perform a (inverse) Fast Fourier Transform.
-
- :param Union[numpy.ndarray,pyopencl.array] array:
- Input data. Must be consistent with the current context.
- :param Union[numpy.ndarray,pyopencl.array] output:
- Output data. By default, output is a numpy.ndarray.
- :param bool do_async:
- Whether to perform operation in asynchronous mode.
- Default is False, meaning that we wait for transform to complete.
- """
- self.set_output_data(array, copy=False)
- self.set_input_data(output, copy=False)
- self.update_inverse_plan_arrays()
- event, = self.plan_inverse.enqueue(forward=False)
- if not(do_async):
- event.wait()
- if output is not None:
- self.copy_output_if_numpy(output, self.data_in)
- res = output
- else:
- res = self.data_in.get()
- self.recover_array_references()
- return res
-
-
- def __del__(self):
- # It seems that gpyfft underlying clFFT destructors are not called.
- # This results in the following warning:
- # Warning: Program terminating, but clFFT resources not freed.
- # Please consider explicitly calling clfftTeardown( )
- del self.plan_forward
- del self.plan_inverse
-
diff --git a/silx/math/fft/cufft.py b/silx/math/fft/cufft.py
deleted file mode 100644
index 848f3e6..0000000
--- a/silx/math/fft/cufft.py
+++ /dev/null
@@ -1,253 +0,0 @@
-#!/usr/bin/env python
-# coding: utf-8
-# /*##########################################################################
-#
-# Copyright (c) 2018-2019 European Synchrotron Radiation Facility
-#
-# Permission is hereby granted, free of charge, to any person obtaining a copy
-# of this software and associated documentation files (the "Software"), to deal
-# in the Software without restriction, including without limitation the rights
-# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-# copies of the Software, and to permit persons to whom the Software is
-# furnished to do so, subject to the following conditions:
-#
-# The above copyright notice and this permission notice shall be included in
-# all copies or substantial portions of the Software.
-#
-# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
-# THE SOFTWARE.
-#
-# ###########################################################################*/
-import numpy as np
-
-from .basefft import BaseFFT
-try:
- import pycuda.gpuarray as gpuarray
- from skcuda.fft import Plan
- from skcuda.fft import fft as cu_fft
- from skcuda.fft import ifft as cu_ifft
- __have_cufft__ = True
-except ImportError:
- __have_cufft__ = False
-
-
-class CUFFT(BaseFFT):
- """Initialize a cufft plan
-
- Please see FFT class for parameters help.
-
- CUFFT-specific parameters
- --------------------------
-
- :param pycuda.driver.Stream stream:
- Stream with which to associate the plan. If no stream is specified,
- the default stream is used.
- """
- def __init__(
- self,
- shape=None,
- dtype=None,
- template=None,
- shape_out=None,
- axes=None,
- normalize="rescale",
- stream=None,
- ):
- if not(__have_cufft__) or not(__have_cufft__):
- raise ImportError("Please install pycuda and scikit-cuda to use the CUDA back-end")
-
- super(CUFFT, self).__init__(
- shape=shape,
- dtype=dtype,
- template=template,
- shape_out=shape_out,
- axes=axes,
- normalize=normalize,
- )
- self.cufft_stream = stream
- self.backend = "cufft"
-
- self.configure_batched_transform()
- self.allocate_arrays()
- self.real_transform = np.isrealobj(self.data_in)
- self.compute_forward_plan()
- self.compute_inverse_plan()
- self.refs = {
- "data_in": self.data_in,
- "data_out": self.data_out,
- }
- self.configure_normalization()
-
- def _allocate(self, shape, dtype):
- return gpuarray.zeros(shape, dtype)
-
- # TODO support batched transform where batch is other than dimension 0
- def configure_batched_transform(self):
- self.cufft_batch_size = 1
- self.cufft_shape = self.shape
- if (self.axes is not None) and (len(self.axes) < len(self.shape)):
- # In the easiest case, the transform is computed along the fastest dimensions:
- # - 1D transforms of lines of 2D data
- # - 2D transforms of images of 3D data (stacked along slow dim)
- # - 1D transforms of 3D data along fastest dim
- # Otherwise, we have to configure cuda "advanced memory layout",
- # which is not implemented yet.
-
- data_ndims = len(self.shape)
- supported_axes = {
- 2: [(1,)],
- 3: [(1, 2), (2, 1), (1,), (2,)],
- }
- if self.axes not in supported_axes[data_ndims]:
- raise NotImplementedError("With the CUDA backend, batched transform is only supported along fastest dimensions")
- self.cufft_batch_size = self.shape[0]
- self.cufft_shape = self.shape[1:]
- if data_ndims == 3 and len(self.axes) == 1:
- # 1D transform on 3D data: here only supported along fast dim,
- # so batch_size is Nx*Ny
- self.cufft_batch_size = np.prod(self.shape[:2])
- self.cufft_shape = (self.shape[-1],)
- if len(self.cufft_shape) == 1:
- self.cufft_shape = self.cufft_shape[0]
-
- def configure_normalization(self):
- # TODO
- if self.normalize == "ortho":
- raise NotImplementedError(
- "Normalization mode 'ortho' is not implemented with CUDA backend yet."
- )
- self.cufft_scale_inverse = (self.normalize == "rescale")
-
- def check_array(self, array, shape, dtype, copy=True):
- if array.shape != shape:
- raise ValueError("Invalid data shape: expected %s, got %s" %
- (shape, array.shape))
- if array.dtype != dtype:
- raise ValueError("Invalid data type: expected %s, got %s" %
- (dtype, array.dtype))
-
- def set_data(self, dst, src, shape, dtype, copy=True, name=None):
- """
- dst is a device array owned by the current instance
- (either self.data_in or self.data_out).
-
- copy is ignored for device<-> arrays.
- """
- self.check_array(src, shape, dtype)
- if isinstance(src, np.ndarray):
- if name == "data_out":
- # Makes little sense to provide output=numpy_array
- return dst
- if not(src.flags["C_CONTIGUOUS"]):
- src = np.ascontiguousarray(src, dtype=dtype)
- dst[:] = src[:]
- elif isinstance(src, gpuarray.GPUArray):
- # No copy, use the data as self.d_input or self.d_output
- # (this prevents the use of in-place transforms, however).
- # We have to keep their old references.
- if name is None:
- # This should not happen
- raise ValueError("Please provide either copy=True or name != None")
- assert id(self.refs[name]) == id(dst) # DEBUG
- setattr(self, name, src)
- return src
- else:
- raise ValueError(
- "Invalid array type %s, expected numpy.ndarray or pycuda.gpuarray" %
- type(src)
- )
- return dst
-
- def recover_array_references(self):
- self.data_in = self.refs["data_in"]
- self.data_out = self.refs["data_out"]
-
- def compute_forward_plan(self):
- self.plan_forward = Plan(
- self.cufft_shape,
- self.dtype,
- self.dtype_out,
- batch=self.cufft_batch_size,
- stream=self.cufft_stream,
- # cufft extensible plan API is only supported after 0.5.1
- # (commit 65288d28ca0b93e1234133f8d460dc6becb65121)
- # but there is still no official 0.5.2
- #~ auto_allocate=True # cufft extensible plan API
- )
-
- def compute_inverse_plan(self):
- self.plan_inverse = Plan(
- self.cufft_shape, # not shape_out
- self.dtype_out,
- self.dtype,
- batch=self.cufft_batch_size,
- stream=self.cufft_stream,
- # cufft extensible plan API is only supported after 0.5.1
- # (commit 65288d28ca0b93e1234133f8d460dc6becb65121)
- # but there is still no official 0.5.2
- #~ auto_allocate=True
- )
-
- def copy_output_if_numpy(self, dst, src):
- if isinstance(dst, gpuarray.GPUArray):
- return
- dst[:] = src[:]
-
- def fft(self, array, output=None):
- """
- Perform a (forward) Fast Fourier Transform.
-
- :param Union[numpy.ndarray,pycuda.gpuarray] array:
- Input data. Must be consistent with the current context.
- :param Union[numpy.ndarray,pycuda.gpuarray] output:
- Output data. By default, output is a numpy.ndarray.
- """
- data_in = self.set_input_data(array, copy=False)
- data_out = self.set_output_data(output, copy=False)
-
- cu_fft(
- data_in,
- data_out,
- self.plan_forward,
- scale=False
- )
-
- if output is not None:
- self.copy_output_if_numpy(output, self.data_out)
- res = output
- else:
- res = self.data_out.get()
- self.recover_array_references()
- return res
-
- def ifft(self, array, output=None):
- """
- Perform a (inverse) Fast Fourier Transform.
-
- :param Union[numpy.ndarray,pycuda.gpuarray] array:
- Input data. Must be consistent with the current context.
- :param Union[numpy.ndarray,pycuda.gpuarray] output:
- Output data. By default, output is a numpy.ndarray.
- """
- data_in = self.set_output_data(array, copy=False)
- data_out = self.set_input_data(output, copy=False)
-
- cu_ifft(
- data_in,
- data_out,
- self.plan_inverse,
- scale=self.cufft_scale_inverse,
- )
-
- if output is not None:
- self.copy_output_if_numpy(output, self.data_in)
- res = output
- else:
- res = self.data_in.get()
- self.recover_array_references()
- return res
diff --git a/silx/math/fft/fft.py b/silx/math/fft/fft.py
deleted file mode 100644
index eb0d73b..0000000
--- a/silx/math/fft/fft.py
+++ /dev/null
@@ -1,96 +0,0 @@
-#!/usr/bin/env python
-# coding: utf-8
-# /*##########################################################################
-#
-# Copyright (c) 2018-2019 European Synchrotron Radiation Facility
-#
-# Permission is hereby granted, free of charge, to any person obtaining a copy
-# of this software and associated documentation files (the "Software"), to deal
-# in the Software without restriction, including without limitation the rights
-# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-# copies of the Software, and to permit persons to whom the Software is
-# furnished to do so, subject to the following conditions:
-#
-# The above copyright notice and this permission notice shall be included in
-# all copies or substantial portions of the Software.
-#
-# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
-# THE SOFTWARE.
-#
-# ###########################################################################*/
-from .fftw import FFTW
-from .clfft import CLFFT
-from .npfft import NPFFT
-from .cufft import CUFFT
-
-
-def FFT(
- shape=None,
- dtype=None,
- template=None,
- shape_out=None,
- axes=None,
- normalize="rescale",
- backend="numpy",
- **kwargs
-):
- """
- Initialize a FFT plan.
-
- :param List[int] shape:
- Shape of the input data.
- :param numpy.dtype dtype:
- Data type of the input data.
- :param numpy.ndarray template:
- Optional data, replacement for "shape" and "dtype".
- If provided, the arguments "shape" and "dtype" are ignored,
- and are instead inferred from it.
- :param List[int] shape_out:
- Optional shape of output data.
- By default, the data has the same shape as the input
- data (in case of C2C transform), or a shape with the last dimension halved
- (in case of R2C transform). If shape_out is provided, it must be greater
- or equal than the shape of input data. In this case, FFT is performed
- with zero-padding.
- :param List[int] axes:
- Axes along which FFT is computed.
- * For 2D transform: axes=(1,0)
- * For batched 1D transform of 2D image: axes=(0,)
- :param str normalize:
- Whether to normalize FFT and IFFT. Possible values are:
- * "rescale": in this case, Fourier data is divided by "N"
- before IFFT, so that (FFT(data)) = data
- * "ortho": in this case, FFT and IFFT are adjoint of eachother,
- the transform is unitary. Both FFT and IFFT are scaled with 1/sqrt(N).
- * "none": no normalizatio is done : IFFT(FFT(data)) = data*N
- :param str backend:
- FFT Backend to use. Value can be "numpy", "fftw", "opencl", "cuda".
- """
- backends = {
- "numpy": NPFFT,
- "np": NPFFT,
- "fftw": FFTW,
- "opencl": CLFFT,
- "clfft": CLFFT,
- "cuda": CUFFT,
- "cufft": CUFFT,
- }
-
- backend = backend.lower()
- if backend not in backends:
- raise ValueError("Unknown backend %s, available are %s" % (backend, backends))
- F = backends[backend](
- shape=shape,
- dtype=dtype,
- template=template,
- shape_out=shape_out,
- axes=axes,
- normalize=normalize,
- **kwargs
- )
- return F
diff --git a/silx/math/fft/fftw.py b/silx/math/fft/fftw.py
deleted file mode 100644
index ff6966c..0000000
--- a/silx/math/fft/fftw.py
+++ /dev/null
@@ -1,214 +0,0 @@
-#!/usr/bin/env python
-# coding: utf-8
-# /*##########################################################################
-#
-# Copyright (c) 2018-2019 European Synchrotron Radiation Facility
-#
-# Permission is hereby granted, free of charge, to any person obtaining a copy
-# of this software and associated documentation files (the "Software"), to deal
-# in the Software without restriction, including without limitation the rights
-# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-# copies of the Software, and to permit persons to whom the Software is
-# furnished to do so, subject to the following conditions:
-#
-# The above copyright notice and this permission notice shall be included in
-# all copies or substantial portions of the Software.
-#
-# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
-# THE SOFTWARE.
-#
-# ###########################################################################*/
-import numpy as np
-
-from .basefft import BaseFFT, check_version
-try:
- import pyfftw
- __have_fftw__ = True
-except ImportError:
- __have_fftw__ = False
-
-
-# Check pyfftw version
-__required_pyfftw_version__ = "0.10.0"
-if __have_fftw__:
- __have_fftw__ = check_version(pyfftw, __required_pyfftw_version__)
-
-
-class FFTW(BaseFFT):
- """Initialize a FFTW plan.
-
- Please see FFT class for parameters help.
-
- FFTW-specific parameters
- -------------------------
-
- :param bool check_alignment:
- If set to True and "data" is provided, this will enforce the input data
- to be "byte aligned", which might imply extra memory usage.
- :param int num_threads:
- Number of threads for computing FFT.
- """
- def __init__(
- self,
- shape=None,
- dtype=None,
- template=None,
- shape_out=None,
- axes=None,
- normalize="rescale",
- check_alignment=False,
- num_threads=1,
- ):
- if not(__have_fftw__):
- raise ImportError("Please install pyfftw >= %s to use the FFTW back-end" % __required_pyfftw_version__)
- super(FFTW, self).__init__(
- shape=shape,
- dtype=dtype,
- template=template,
- shape_out=shape_out,
- axes=axes,
- normalize=normalize,
- )
- self.check_alignment = check_alignment
- self.num_threads = num_threads
- self.backend = "fftw"
-
- self.allocate_arrays()
- self.set_fftw_flags()
- self.compute_forward_plan()
- self.compute_inverse_plan()
- self.refs = {
- "data_in": self.data_in,
- "data_out": self.data_out,
- }
-
- def set_fftw_flags(self):
- self.fftw_flags = ('FFTW_MEASURE', ) # TODO
- self.fftw_planning_timelimit = None # TODO
- self.fftw_norm_modes = {
- "rescale": {"ortho": False, "normalize": True},
- "ortho": {"ortho": True, "normalize": False},
- "none": {"ortho": False, "normalize": False},
- }
- if self.normalize not in self.fftw_norm_modes:
- raise ValueError("Unknown normalization mode %s. Possible values are %s" %
- (self.normalize, self.fftw_norm_modes.keys())
- )
- self.fftw_norm_mode = self.fftw_norm_modes[self.normalize]
-
- def _allocate(self, shape, dtype):
- return pyfftw.zeros_aligned(shape, dtype=dtype)
-
- def check_array(self, array, shape, dtype, copy=True):
- if array.shape != shape:
- raise ValueError("Invalid data shape: expected %s, got %s" %
- (shape, array.shape)
- )
- if array.dtype != dtype:
- raise ValueError("Invalid data type: expected %s, got %s" %
- (dtype, array.dtype)
- )
-
- def set_data(self, self_array, array, shape, dtype, copy=True, name=None):
- """
- :param self_array: array owned by the current instance
- (either self.data_in or self.data_out).
- :type: numpy.ndarray
- :param self_array: data to set
- :type: numpy.ndarray
- :type tuple shape: shape of the array
- :param dtype: type of the array
- :type: numpy.dtype
- :param bool copy: should we copy the array
- :param str name: name of the array
-
- Copies are avoided when possible.
- """
- self.check_array(array, shape, dtype)
- if id(self.refs[name]) == id(array):
- # nothing to do: fft is performed on self.data_in or self.data_out
- arr_to_use = self.refs[name]
- if self.check_alignment and not(pyfftw.is_byte_aligned(array)):
- # If the array is not properly aligned,
- # create a temp. array copy it to self.data_in or self.data_out
- self_array[:] = array[:]
- arr_to_use = self_array
- else:
- # If the array is properly aligned, use it directly
- if copy:
- arr_to_use = np.copy(array)
- else:
- arr_to_use = array
- return arr_to_use
-
- def compute_forward_plan(self):
- self.plan_forward = pyfftw.FFTW(
- self.data_in,
- self.data_out,
- axes=self.axes,
- direction='FFTW_FORWARD',
- flags=self.fftw_flags,
- threads=self.num_threads,
- planning_timelimit=self.fftw_planning_timelimit,
- # the following seems to be taken into account only when using __call__
- ortho=self.fftw_norm_mode["ortho"],
- normalise_idft=self.fftw_norm_mode["normalize"],
- )
-
- def compute_inverse_plan(self):
- self.plan_inverse = pyfftw.FFTW(
- self.data_out,
- self.data_in,
- axes=self.axes,
- direction='FFTW_BACKWARD',
- flags=self.fftw_flags,
- threads=self.num_threads,
- planning_timelimit=self.fftw_planning_timelimit,
- # the following seem to be taken into account only when using __call__
- ortho=self.fftw_norm_mode["ortho"],
- normalise_idft=self.fftw_norm_mode["normalize"],
- )
-
- def fft(self, array, output=None):
- """
- Perform a (forward) Fast Fourier Transform.
-
- :param numpy.ndarray array:
- Input data. Must be consistent with the current context.
- :param numpy.ndarray output:
- Optional output data.
- """
- data_in = self.set_input_data(array, copy=False)
- data_out = self.set_output_data(output, copy=False)
- self.plan_forward.update_arrays(data_in, data_out)
- # execute.__call__ does both update_arrays() and normalization
- self.plan_forward(
- ortho=self.fftw_norm_mode["ortho"],
- )
- self.plan_forward.update_arrays(self.refs["data_in"], self.refs["data_out"])
- return data_out
-
- def ifft(self, array, output=None):
- """
- Perform a (inverse) Fast Fourier Transform.
-
- :param numpy.ndarray array:
- Input data. Must be consistent with the current context.
- :param numpy.ndarray output:
- Optional output data.
- """
- data_in = self.set_output_data(array, copy=False)
- data_out = self.set_input_data(output, copy=False)
- self.plan_inverse.update_arrays(data_in, data_out)
- # execute.__call__ does both update_arrays() and normalization
- self.plan_inverse(
- ortho=self.fftw_norm_mode["ortho"],
- normalise_idft=self.fftw_norm_mode["normalize"]
- )
- self.plan_inverse.update_arrays(self.refs["data_out"], self.refs["data_in"])
- return data_out
diff --git a/silx/math/fft/npfft.py b/silx/math/fft/npfft.py
deleted file mode 100644
index 20351de..0000000
--- a/silx/math/fft/npfft.py
+++ /dev/null
@@ -1,124 +0,0 @@
-#!/usr/bin/env python
-# coding: utf-8
-# /*##########################################################################
-#
-# Copyright (c) 2018-2019 European Synchrotron Radiation Facility
-#
-# Permission is hereby granted, free of charge, to any person obtaining a copy
-# of this software and associated documentation files (the "Software"), to deal
-# in the Software without restriction, including without limitation the rights
-# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-# copies of the Software, and to permit persons to whom the Software is
-# furnished to do so, subject to the following conditions:
-#
-# The above copyright notice and this permission notice shall be included in
-# all copies or substantial portions of the Software.
-#
-# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
-# THE SOFTWARE.
-#
-# ###########################################################################*/
-import numpy as np
-
-from .basefft import BaseFFT
-
-
-class NPFFT(BaseFFT):
- """Initialize a numpy plan.
-
- Please see FFT class for parameters help.
- """
- def __init__(
- self,
- shape=None,
- dtype=None,
- template=None,
- shape_out=None,
- axes=None,
- normalize="rescale",
- ):
- super(NPFFT, self).__init__(
- shape=shape,
- dtype=dtype,
- template=template,
- shape_out=shape_out,
- axes=axes,
- normalize=normalize,
- )
- self.backend = "numpy"
- self.real_transform = False
- if template is not None and np.isrealobj(template):
- self.real_transform = True
- # For numpy functions.
- # TODO Issue warning if user wants ifft(fft(data)) = N*data ?
- if normalize != "ortho":
- self.normalize = None
- self.set_fft_functions()
- #~ self.allocate_arrays() # not needed for this backend
- self.compute_plans()
-
-
- def set_fft_functions(self):
- # (fwd, inv) = _fft_functions[is_real][ndim]
- self._fft_functions = {
- True: {
- 1: (np.fft.rfft, np.fft.irfft),
- 2: (np.fft.rfft2, np.fft.irfft2),
- 3: (np.fft.rfftn, np.fft.irfftn),
- },
- False: {
- 1: (np.fft.fft, np.fft.ifft),
- 2: (np.fft.fft2, np.fft.ifft2),
- 3: (np.fft.fftn, np.fft.ifftn),
- }
- }
-
-
- def _allocate(self, shape, dtype):
- return np.zeros(self.queue, shape, dtype=dtype)
-
-
- def compute_plans(self):
- ndim = len(self.shape)
- funcs = self._fft_functions[self.real_transform][np.minimum(ndim, 3)]
- if np.version.version[:4] in ["1.8.", "1.9."]:
- # norm keyword was introduced in 1.10 and we support numpy >= 1.8
- self.numpy_args = {}
- else:
- self.numpy_args = {"norm": self.normalize}
- # Batched transform
- if (self.user_axes is not None) and len(self.user_axes) < ndim:
- funcs = self._fft_functions[self.real_transform][np.minimum(ndim-1, 3)]
- self.numpy_args["axes"] = self.user_axes
- # Special case of batched 1D transform on 2D data
- if ndim == 2:
- assert len(self.user_axes) == 1
- self.numpy_args["axis"] = self.user_axes[0]
- self.numpy_args.pop("axes")
- self.numpy_funcs = funcs
-
-
- def fft(self, array):
- """
- Perform a (forward) Fast Fourier Transform.
-
- :param numpy.ndarray array:
- Input data. Must be consistent with the current context.
- """
- return self.numpy_funcs[0](array, **self.numpy_args)
-
-
- def ifft(self, array):
- """
- Perform a (inverse) Fast Fourier Transform.
-
- :param numpy.ndarray array:
- Input data. Must be consistent with the current context.
- """
- return self.numpy_funcs[1](array, **self.numpy_args)
-
diff --git a/silx/math/fft/setup.py b/silx/math/fft/setup.py
deleted file mode 100644
index 76bb864..0000000
--- a/silx/math/fft/setup.py
+++ /dev/null
@@ -1,41 +0,0 @@
-# coding: utf-8
-# /*##########################################################################
-# Copyright (C) 2016-2017 European Synchrotron Radiation Facility
-#
-# Permission is hereby granted, free of charge, to any person obtaining a copy
-# of this software and associated documentation files (the "Software"), to deal
-# in the Software without restriction, including without limitation the rights
-# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-# copies of the Software, and to permit persons to whom the Software is
-# furnished to do so, subject to the following conditions:
-#
-# The above copyright notice and this permission notice shall be included in
-# all copies or substantial portions of the Software.
-#
-# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
-# THE SOFTWARE.
-#
-# ############################################################################*/
-
-__authors__ = ["P. Naudet"]
-__license__ = "MIT"
-__date__ = "12/12/2018"
-
-import numpy
-from numpy.distutils.misc_util import Configuration
-
-
-def configuration(parent_package='', top_path=None):
- config = Configuration('fft', parent_package, top_path)
- config.add_subpackage('test')
- return config
-
-
-if __name__ == "__main__":
- from numpy.distutils.core import setup
- setup(configuration=configuration)
diff --git a/silx/math/fft/test/__init__.py b/silx/math/fft/test/__init__.py
deleted file mode 100644
index 83f8926..0000000
--- a/silx/math/fft/test/__init__.py
+++ /dev/null
@@ -1,25 +0,0 @@
-# coding: utf-8
-# /*##########################################################################
-# Copyright (C) 2016-2019 European Synchrotron Radiation Facility
-#
-# Permission is hereby granted, free of charge, to any person obtaining a copy
-# of this software and associated documentation files (the "Software"), to deal
-# in the Software without restriction, including without limitation the rights
-# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-# copies of the Software, and to permit persons to whom the Software is
-# furnished to do so, subject to the following conditions:
-#
-# The above copyright notice and this permission notice shall be included in
-# all copies or substantial portions of the Software.
-#
-# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
-# THE SOFTWARE.
-#
-# ############################################################################*/
-
-from .test_fft import suite
diff --git a/silx/math/fft/test/test_fft.py b/silx/math/fft/test/test_fft.py
deleted file mode 100644
index 9ef2fd2..0000000
--- a/silx/math/fft/test/test_fft.py
+++ /dev/null
@@ -1,270 +0,0 @@
-#!/usr/bin/env python
-# coding: utf-8
-# /*##########################################################################
-#
-# Copyright (c) 2018-2019 European Synchrotron Radiation Facility
-#
-# Permission is hereby granted, free of charge, to any person obtaining a copy
-# of this software and associated documentation files (the "Software"), to deal
-# in the Software without restriction, including without limitation the rights
-# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-# copies of the Software, and to permit persons to whom the Software is
-# furnished to do so, subject to the following conditions:
-#
-# The above copyright notice and this permission notice shall be included in
-# all copies or substantial portions of the Software.
-#
-# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
-# THE SOFTWARE.
-#
-# ###########################################################################*/
-"""Test of the FFT module"""
-
-import numpy as np
-import unittest
-import logging
-try:
- from scipy.misc import ascent
- __have_scipy = True
-except ImportError:
- __have_scipy = False
-from silx.utils.testutils import ParametricTestCase
-from silx.math.fft.fft import FFT
-from silx.math.fft.clfft import __have_clfft__
-from silx.math.fft.cufft import __have_cufft__
-from silx.math.fft.fftw import __have_fftw__
-
-from silx.test.utils import test_options
-
-logger = logging.getLogger(__name__)
-
-
-class TransformInfos(object):
- def __init__(self):
- self.dimensions = [
- "1D",
- "batched_1D",
- "2D",
- "batched_2D",
- "3D",
- ]
- self.modes = {
- "R2C": np.float32,
- "R2C_double": np.float64,
- "C2C": np.complex64,
- "C2C_double": np.complex128,
- }
- self.sizes = {
- "1D": [(128,), (127,)],
- "2D": [(128, 128), (128, 127), (127, 128), (127, 127)],
- "3D": [(64, 64, 64), (64, 64, 63), (64, 63, 64), (63, 64, 64),
- (64, 63, 63), (63, 64, 63), (63, 63, 64), (63, 63, 63)]
- }
- self.axes = {
- "1D": None,
- "batched_1D": (-1,),
- "2D": None,
- "batched_2D": (-2, -1),
- "3D": None,
- }
- self.sizes["batched_1D"] = self.sizes["2D"]
- self.sizes["batched_2D"] = self.sizes["3D"]
-
-
-class TestData(object):
- def __init__(self):
- self.data = ascent().astype("float32")
- self.data1d = self.data[:, 0] # non-contiguous data
- self.data3d = np.tile(self.data[:64, :64], (64, 1, 1))
- self.data_refs = {
- 1: self.data1d,
- 2: self.data,
- 3: self.data3d,
- }
-
-
-@unittest.skipUnless(__have_scipy, "scipy is missing")
-class TestFFT(ParametricTestCase):
- """Test cuda/opencl/fftw backends of FFT"""
-
- def setUp(self):
- self.tol = {
- np.dtype("float32"): 1e-3,
- np.dtype("float64"): 1e-9,
- np.dtype("complex64"): 1e-3,
- np.dtype("complex128"): 1e-9,
- }
- self.transform_infos = TransformInfos()
- self.test_data = TestData()
-
- @staticmethod
- def calc_mae(arr1, arr2):
- """
- Compute the Max Absolute Error between two arrays
- """
- return np.max(np.abs(arr1 - arr2))
-
- @unittest.skipIf(not __have_cufft__,
- "cuda back-end requires pycuda and scikit-cuda")
- def test_cuda(self):
- import pycuda.autoinit
-
- # Error is higher when using cuda. fast_math mode ?
- self.tol[np.dtype("float32")] *= 2
-
- self.__run_tests(backend="cuda")
-
- @unittest.skipIf(not __have_clfft__,
- "opencl back-end requires pyopencl and gpyfft")
- def test_opencl(self):
- from silx.opencl.common import ocl
- if ocl is not None:
- self.__run_tests(backend="opencl", ctx=ocl.create_context())
-
- @unittest.skipIf(not __have_fftw__,
- "fftw back-end requires pyfftw")
- def test_fftw(self):
- self.__run_tests(backend="fftw")
-
- def __run_tests(self, backend, **extra_args):
- """Run all tests with the given backend
-
- :param str backend:
- :param dict extra_args: Additional arguments to provide to FFT
- """
- for trdim in self.transform_infos.dimensions:
- for mode in self.transform_infos.modes:
- for size in self.transform_infos.sizes[trdim]:
- with self.subTest(trdim=trdim, mode=mode, size=size):
- self.__test(backend, trdim, mode, size, **extra_args)
-
- def __test(self, backend, trdim, mode, size, **extra_args):
- """Compare given backend with numpy for given conditions"""
- logger.debug("backend: %s, trdim: %s, mode: %s, size: %s",
- backend, trdim, mode, str(size))
- if size == "3D" and test_options.TEST_LOW_MEM:
- self.skipTest("low mem")
-
- ndim = len(size)
- input_data = self.test_data.data_refs[ndim].astype(
- self.transform_infos.modes[mode])
- tol = self.tol[np.dtype(input_data.dtype)]
- if trdim == "3D":
- tol *= 10 # Error is relatively high in high dimensions
- # It seems that cuda has problems with C2D batched 1D
- if trdim == "batched_1D" and backend == "cuda" and mode == "C2C":
- tol *= 10
-
- # Python < 3.5 does not want to mix **extra_args with existing kwargs
- fft_args = {
- "template": input_data,
- "axes": self.transform_infos.axes[trdim],
- "backend": backend,
- }
- fft_args.update(extra_args)
- F = FFT(
- **fft_args
- )
- F_np = FFT(
- template=input_data,
- axes=self.transform_infos.axes[trdim],
- backend="numpy"
- )
-
- # Forward FFT
- res = F.fft(input_data)
- res_np = F_np.fft(input_data)
- mae = self.calc_mae(res, res_np)
- all_close = np.allclose(res, res_np, atol=tol, rtol=tol),
- self.assertTrue(
- all_close,
- "FFT %s:%s, MAE(%s, numpy) = %f (tol = %.2e)" % (mode, trdim, backend, mae, tol)
- )
-
- # Inverse FFT
- res2 = F.ifft(res)
- mae = self.calc_mae(res2, input_data)
- self.assertTrue(
- mae < tol,
- "IFFT %s:%s, MAE(%s, numpy) = %f" % (mode, trdim, backend, mae)
- )
-
-
-@unittest.skipUnless(__have_scipy, "scipy is missing")
-class TestNumpyFFT(ParametricTestCase):
- """
- Test the Numpy backend individually.
- """
-
- def setUp(self):
- transforms = {
- "1D": {
- True: (np.fft.rfft, np.fft.irfft),
- False: (np.fft.fft, np.fft.ifft),
- },
- "2D": {
- True: (np.fft.rfft2, np.fft.irfft2),
- False: (np.fft.fft2, np.fft.ifft2),
- },
- "3D": {
- True: (np.fft.rfftn, np.fft.irfftn),
- False: (np.fft.fftn, np.fft.ifftn),
- },
- }
- transforms["batched_1D"] = transforms["1D"]
- transforms["batched_2D"] = transforms["2D"]
- self.transforms = transforms
- self.transform_infos = TransformInfos()
- self.test_data = TestData()
-
- def test(self):
- """Test the numpy backend against native fft.
-
- Results should be exactly the same.
- """
- for trdim in self.transform_infos.dimensions:
- for mode in self.transform_infos.modes:
- for size in self.transform_infos.sizes[trdim]:
- with self.subTest(trdim=trdim, mode=mode, size=size):
- self.__test(trdim, mode, size)
-
- def __test(self, trdim, mode, size):
- logger.debug("trdim: %s, mode: %s, size: %s", trdim, mode, str(size))
- ndim = len(size)
- input_data = self.test_data.data_refs[ndim].astype(
- self.transform_infos.modes[mode])
- np_fft, np_ifft = self.transforms[trdim][np.isrealobj(input_data)]
-
- F = FFT(
- template=input_data,
- axes=self.transform_infos.axes[trdim],
- backend="numpy"
- )
- # Test FFT
- res = F.fft(input_data)
- ref = np_fft(input_data)
- self.assertTrue(np.allclose(res, ref))
-
- # Test IFFT
- res2 = F.ifft(res)
- ref2 = np_ifft(ref)
- self.assertTrue(np.allclose(res2, ref2))
-
-
-def suite():
- suite = unittest.TestSuite()
- for cls in (TestNumpyFFT, TestFFT):
- suite.addTest(
- unittest.defaultTestLoader.loadTestsFromTestCase(cls))
- return suite
-
-
-if __name__ == '__main__':
- unittest.main(defaultTest="suite")
-
-