From 607b253b3f856acc4103ddd7988e2842341a88f6 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Sun, 31 Aug 2025 00:00:56 +0530 Subject: [PATCH 01/19] Add `mkl_fft` in fft --- environment-dev-arm.yml | 1 + environment-dev-gpu.yml | 1 + environment-dev.yml | 1 + examples/plot_fft.py | 24 ++++++ pylops/signalprocessing/fft.py | 123 ++++++++++++++++++++++++++- pylops/signalprocessing/fft2d.py | 141 ++++++++++++++++++++++++++++++- pylops/signalprocessing/fftnd.py | 125 ++++++++++++++++++++++++++- pylops/utils/deps.py | 20 +++++ pytests/test_ffts.py | 114 ++++++++++++++++++++++--- requirements-dev-gpu.txt | 1 + requirements-dev.txt | 1 + requirements-doc.txt | 1 + 12 files changed, 533 insertions(+), 20 deletions(-) diff --git a/environment-dev-arm.yml b/environment-dev-arm.yml index 49ae8824..ae485474 100644 --- a/environment-dev-arm.yml +++ b/environment-dev-arm.yml @@ -23,6 +23,7 @@ dependencies: - autopep8 - isort - black + - mkl_fft - pip: - torch - devito diff --git a/environment-dev-gpu.yml b/environment-dev-gpu.yml index 8436d228..1ae98691 100644 --- a/environment-dev-gpu.yml +++ b/environment-dev-gpu.yml @@ -21,6 +21,7 @@ dependencies: - autopep8 - isort - black + - mkl_fft - pip: - torch - pytest-runner diff --git a/environment-dev.yml b/environment-dev.yml index fccd5663..8ca130c0 100644 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -24,6 +24,7 @@ dependencies: - autopep8 - isort - black + - mkl_fft - pip: - torch - devito diff --git a/examples/plot_fft.py b/examples/plot_fft.py index a8af9da6..5718e252 100644 --- a/examples/plot_fft.py +++ b/examples/plot_fft.py @@ -66,6 +66,30 @@ axs[1].set_xlim([0, 3 * f0]) plt.tight_layout() +############################################################################### +# PyLops also has a third FFT engine (engine='mkl_fft') that uses the well-known +# `Intel MKL FFT `_. This is a Python wrapper around +# the `Intel® oneAPI Math Kernel Library (oneMKL) `_ +# Fourier Transform functions. It lets PyLops run discrete Fourier transforms faster +# by using Intel’s highly optimized math routines. + +FFTop = pylops.signalprocessing.FFT(dims=nt, nfft=nfft, sampling=dt, engine="mkl_fft") +D = FFTop * d + +# Inverse for FFT +dinv = FFTop.H * D +dinv = FFTop / D + +fig, axs = plt.subplots(1, 2, figsize=(10, 4)) +axs[0].plot(t, d, "k", lw=2, label="True") +axs[0].plot(t, dinv.real, "--r", lw=2, label="Inverted") +axs[0].legend() +axs[0].set_title("Signal") +axs[1].plot(FFTop.f[: int(FFTop.nfft / 2)], np.abs(D[: int(FFTop.nfft / 2)]), "k", lw=2) +axs[1].set_title("Fourier Transform with MKL FFT") +axs[1].set_xlim([0, 3 * f0]) +plt.tight_layout() + ############################################################################### # We can also apply the one dimensional FFT to to a two-dimensional # signal (along one of the first axis) diff --git a/pylops/signalprocessing/fft.py b/pylops/signalprocessing/fft.py index fe853771..5df2ed1f 100644 --- a/pylops/signalprocessing/fft.py +++ b/pylops/signalprocessing/fft.py @@ -16,10 +16,20 @@ from pylops.utils.typing import DTypeLike, InputDimsLike, NDArray pyfftw_message = deps.pyfftw_import("the fft module") +mkl_fft_message = deps.mkl_fft_import("the mkl fft module") if pyfftw_message is None: import pyfftw +if mkl_fft_message is None: + import mkl_fft.interfaces.numpy_fft as mkl_backend + + try: + import scipy.fft # noqa: F401 + import mkl_fft.interfaces.scipy_fft as mkl_backend + except ImportError: + pass + logger = logging.getLogger(__name__) @@ -394,6 +404,94 @@ def __truediv__(self, y: npt.ArrayLike) -> npt.ArrayLike: return self._rmatvec(y) / self._scale +class _FFT_mklfft(_BaseFFT): + """One-dimensional Fast-Fourier Transform using mkl_fft""" + + def __init__( + self, + dims: Union[int, InputDimsLike], + axis: int = -1, + nfft: Optional[int] = None, + sampling: float = 1.0, + norm: str = "ortho", + real: bool = False, + ifftshift_before: bool = False, + fftshift_after: bool = False, + dtype: DTypeLike = "complex128", + **kwargs_fft, + ) -> None: + super().__init__( + dims=dims, + axis=axis, + nfft=nfft, + sampling=sampling, + norm=norm, + real=real, + ifftshift_before=ifftshift_before, + fftshift_after=fftshift_after, + dtype=dtype, + ) + self._kwargs_fft = kwargs_fft + self._norm_kwargs = {"norm": None} + if self.norm is _FFTNorms.ORTHO: + self._norm_kwargs["norm"] = "ortho" + self._scale = np.sqrt(1 / self.nfft) + elif self.norm is _FFTNorms.NONE: + self._scale = self.nfft + elif self.norm is _FFTNorms.ONE_OVER_N: + self._scale = 1.0 / self.nfft + + @reshaped + def _matvec(self, x: NDArray) -> NDArray: + if self.ifftshift_before: + x = mkl_backend.ifftshift(x, axes=self.axis) + if not self.clinear: + x = np.real(x) + if self.real: + y = mkl_backend.rfft(x, n=self.nfft, axis=self.axis, **self._norm_kwargs) + y = np.swapaxes(y, -1, self.axis) + y[..., 1 : 1 + (self.nfft - 1) // 2] *= np.sqrt(2) + y = np.swapaxes(y, self.axis, -1) + else: + y = mkl_backend.fft(x, n=self.nfft, axis=self.axis, **self._norm_kwargs) + if self.norm is _FFTNorms.ONE_OVER_N: + y *= self._scale + if self.fftshift_after: + y = mkl_backend.fftshift(y, axes=self.axis) + return y + + @reshaped + def _rmatvec(self, x: NDArray) -> NDArray: + if self.fftshift_after: + x = mkl_backend.ifftshift(x, axes=self.axis) + if self.real: + x = x.copy() + x = np.swapaxes(x, -1, self.axis) + x[..., 1 : 1 + (self.nfft - 1) // 2] /= np.sqrt(2) + x = np.swapaxes(x, self.axis, -1) + y = mkl_backend.irfft(x, n=self.nfft, axis=self.axis, **self._norm_kwargs) + else: + y = mkl_backend.ifft(x, n=self.nfft, axis=self.axis, **self._norm_kwargs) + if self.norm is _FFTNorms.NONE: + y *= self._scale + + if self.nfft > self.dims[self.axis]: + y = np.take(y, range(0, self.dims[self.axis]), axis=self.axis) + elif self.nfft < self.dims[self.axis]: + y = np.pad(y, self.ifftpad) + + if not self.clinear: + y = np.real(y) + if self.ifftshift_before: + y = mkl_backend.fftshift(y, axes=self.axis) + return y + + def __truediv__(self, y): + if self.norm is not _FFTNorms.ORTHO: + return self._rmatvec(y) / self._scale + return self._rmatvec(y) + + def FFT( dims: Union[int, InputDimsLike], axis: int = -1, @@ -481,7 +579,7 @@ def FFT( frequencies are arranged from zero to largest positive, and then from negative Nyquist to the frequency bin before zero. engine : :obj:`str`, optional - Engine used for fft computation (``numpy``, ``fftw``, or ``scipy``). Choose + Engine used for fft computation (``numpy``, ``fftw``, ``scipy`` or ``mkl_fft``). Choose ``numpy`` when working with cupy and jax arrays. .. note:: Since version 1.17.0, accepts "scipy". @@ -534,7 +632,7 @@ def FFT( - If ``dims`` is provided and ``axis`` is bigger than ``len(dims)``. - If ``norm`` is not one of "ortho", "none", or "1/n". NotImplementedError - If ``engine`` is neither ``numpy``, ``fftw``, nor ``scipy``. + If ``engine`` is neither ``numpy``, ``fftw``, ``scipy`` nor ``mkl_fft``. See Also -------- @@ -579,7 +677,24 @@ def FFT( dtype=dtype, **kwargs_fft, ) - elif engine == "numpy" or (engine == "fftw" and pyfftw_message is not None): + elif engine == "mkl_fft" and mkl_fft_message is None: + f = _FFT_mklfft( + dims, + axis=axis, + nfft=nfft, + sampling=sampling, + norm=norm, + real=real, + ifftshift_before=ifftshift_before, + fftshift_after=fftshift_after, + dtype=dtype, + **kwargs_fft, + ) + elif ( + engine == "numpy" + or (engine == "fftw" and pyfftw_message is not None) + or (engine == "mkl_fft" and mkl_fft_message is not None) + ): if engine == "fftw" and pyfftw_message is not None: logger.warning(pyfftw_message) f = _FFT_numpy( @@ -608,6 +723,6 @@ def FFT( **kwargs_fft, ) else: - raise NotImplementedError("engine must be numpy, fftw or scipy") + raise NotImplementedError("engine must be numpy, fftw, scipy or mkl_fft") f.name = name return f diff --git a/pylops/signalprocessing/fft2d.py b/pylops/signalprocessing/fft2d.py index 5df5df34..c21d47ff 100644 --- a/pylops/signalprocessing/fft2d.py +++ b/pylops/signalprocessing/fft2d.py @@ -8,10 +8,22 @@ from pylops import LinearOperator from pylops.signalprocessing._baseffts import _BaseFFTND, _FFTNorms +from pylops.utils import deps from pylops.utils.backend import get_array_module from pylops.utils.decorators import reshaped from pylops.utils.typing import DTypeLike, InputDimsLike +mkl_fft_message = deps.mkl_fft_import("the mkl fft module") + +if mkl_fft_message is None: + import mkl_fft.interfaces.numpy_fft as mkl_backend + + try: + import scipy.fft + import mkl_fft.interfaces.scipy_fft as mkl_backend + except ImportError: + pass + class _FFT2D_numpy(_BaseFFTND): """Two dimensional Fast-Fourier Transform using NumPy""" @@ -235,6 +247,115 @@ def __truediv__(self, y): return self._rmatvec(y) +class _FFT2D_mklfft(_BaseFFTND): + """Two-dimensional Fast-Fourier Transform using mkl_fft""" + + def __init__( + self, + dims: InputDimsLike, + axes: InputDimsLike = (-2, -1), + nffts: Optional[Union[int, InputDimsLike]] = None, + sampling: Union[float, Sequence[float]] = 1.0, + norm: str = "ortho", + real: bool = False, + ifftshift_before: bool = False, + fftshift_after: bool = False, + dtype: DTypeLike = "complex128", + **kwargs_fft, + ) -> None: + super().__init__( + dims=dims, + axes=axes, + nffts=nffts, + sampling=sampling, + norm=norm, + real=real, + ifftshift_before=ifftshift_before, + fftshift_after=fftshift_after, + dtype=dtype, + ) + + # checks + if self.ndim < 2: + raise ValueError("FFT2D requires at least two input dimensions") + if self.naxes != 2: + raise ValueError("FFT2D must be applied along exactly two dimensions") + + self._kwargs_fft = kwargs_fft + self.f1, self.f2 = self.fs + del self.fs + + self._norm_kwargs: Dict[str, Union[None, str]] = {"norm": None} + if self.norm is _FFTNorms.ORTHO: + self._norm_kwargs["norm"] = "ortho" + self._scale = np.sqrt(1 / np.prod(np.sqrt(self.nffts))) + elif self.norm is _FFTNorms.NONE: + self._scale = np.sqrt(np.prod(self.nffts)) + elif self.norm is _FFTNorms.ONE_OVER_N: + self._scale = np.sqrt(1.0 / np.prod(self.nffts)) + + @reshaped + def _matvec(self, x): + if self.ifftshift_before.any(): + x = mkl_backend.ifftshift(x, axes=self.axes[self.ifftshift_before]) + if not self.clinear: + x = np.real(x) + if self.real: + y = mkl_backend.rfft2( + x, s=self.nffts, axes=self.axes, **self._norm_kwargs, **self._kwargs_fft + ) + y = np.swapaxes(y, -1, self.axes[-1]) + y[..., 1 : 1 + (self.nffts[-1] - 1) // 2] *= np.sqrt(2) + y = np.swapaxes(y, self.axes[-1], -1) + else: + y = mkl_backend.fft2( + x, s=self.nffts, axes=self.axes, **self._norm_kwargs, **self._kwargs_fft + ) + if self.norm is _FFTNorms.ONE_OVER_N: + y *= self._scale + y = y.astype(self.cdtype) + if self.fftshift_after.any(): + y = mkl_backend.fftshift(y, axes=self.axes[self.fftshift_after]) + return y + + @reshaped + def _rmatvec(self, x): + if self.fftshift_after.any(): + x = mkl_backend.ifftshift(x, axes=self.axes[self.fftshift_after]) + if self.real: + x = x.copy() + x = np.swapaxes(x, -1, self.axes[-1]) + x[..., 1 : 1 + (self.nffts[-1] - 1) // 2] /= np.sqrt(2) + x = np.swapaxes(x, self.axes[-1], -1) + y = mkl_backend.irfft2( + x, s=self.nffts, axes=self.axes, **self._norm_kwargs, **self._kwargs_fft + ) + else: + y = mkl_backend.ifft2( + x, s=self.nffts, axes=self.axes, **self._norm_kwargs, **self._kwargs_fft + ) + if self.norm is _FFTNorms.NONE: + y *= self._scale + print(y.shape, self.dims[self.axes[0]]) + if self.nffts[0] > self.dims[self.axes[0]]: + y = np.take(y, np.arange(self.dims[self.axes[0]]), axis=self.axes[0]) + if self.nffts[1] > self.dims[self.axes[1]]: + y = np.take(y, np.arange(self.dims[self.axes[1]]), axis=self.axes[1]) + if self.doifftpad: + y = np.pad(y, self.ifftpad) + if not self.clinear: + y = np.real(y) + y = y.astype(self.rdtype) + if self.ifftshift_before.any(): + y = np.fft.fftshift(y, axes=self.axes[self.ifftshift_before]) + return y + + def __truediv__(self, y): + if self.norm is not _FFTNorms.ORTHO: + return self._rmatvec(y) / self._scale / self._scale + return self._rmatvec(y) + + def FFT2D( dims: InputDimsLike, axes: InputDimsLike = (-2, -1), @@ -331,7 +452,7 @@ def FFT2D( engine : :obj:`str`, optional .. versionadded:: 1.17.0 - Engine used for fft computation (``numpy`` or ``scipy``). Choose + Engine used for fft computation (``numpy`` or ``scipy`` or ``mkl_fft``). Choose ``numpy`` when working with cupy and jax arrays. dtype : :obj:`str`, optional Type of elements in input array. Note that the ``dtype`` of the operator @@ -387,7 +508,7 @@ def FFT2D( two elements. - If ``norm`` is not one of "ortho", "none", or "1/n". NotImplementedError - If ``engine`` is neither ``numpy``, nor ``scipy``. + If ``engine`` is neither ``numpy``, ``scipy`` nor ``mkl_fft``. See Also -------- @@ -420,7 +541,19 @@ def FFT2D( signals. """ - if engine == "numpy": + if engine == "mkl_fft" and mkl_fft_message is None: + f = _FFT2D_mklfft( + dims=dims, + axes=axes, + nffts=nffts, + sampling=sampling, + norm=norm, + real=real, + ifftshift_before=ifftshift_before, + fftshift_after=fftshift_after, + dtype=dtype, + ) + elif engine == "numpy" or (engine == "mkl_fft" and mkl_fft_message is not None): f = _FFT2D_numpy( dims=dims, axes=axes, @@ -447,6 +580,6 @@ def FFT2D( **kwargs_fft, ) else: - raise NotImplementedError("engine must be numpy or scipy") + raise NotImplementedError("engine must be numpy, scipy or mkl_fft") f.name = name return f diff --git a/pylops/signalprocessing/fftnd.py b/pylops/signalprocessing/fftnd.py index 5608b804..2a47e505 100644 --- a/pylops/signalprocessing/fftnd.py +++ b/pylops/signalprocessing/fftnd.py @@ -8,10 +8,22 @@ from pylops import LinearOperator from pylops.signalprocessing._baseffts import _BaseFFTND, _FFTNorms +from pylops.utils import deps from pylops.utils.backend import get_array_module, get_sp_fft from pylops.utils.decorators import reshaped from pylops.utils.typing import DTypeLike, InputDimsLike, NDArray +mkl_fft_message = deps.mkl_fft_import("the mkl fft module") + +if mkl_fft_message is None: + import mkl_fft.interfaces.numpy_fft as mkl_backend + + try: + import scipy.fft + import mkl_fft.interfaces.scipy_fft as mkl_backend + except ImportError: + pass + class _FFTND_numpy(_BaseFFTND): """N-dimensional Fast-Fourier Transform using NumPy""" @@ -216,6 +228,102 @@ def __truediv__(self, y: npt.ArrayLike) -> npt.ArrayLike: return self._rmatvec(y) +class _FFTND_mklfft(_BaseFFTND): + """N-dimensional Fast-Fourier Transform using MKL FFT""" + + def __init__( + self, + dims: Union[int, InputDimsLike], + axes: Union[int, InputDimsLike] = (-3, -2, -1), + nffts: Optional[Union[int, InputDimsLike]] = None, + sampling: Union[float, Sequence[float]] = 1.0, + norm: str = "ortho", + real: bool = False, + ifftshift_before: bool = False, + fftshift_after: bool = False, + dtype: DTypeLike = "complex128", + **kwargs_fft, + ) -> None: + super().__init__( + dims=dims, + axes=axes, + nffts=nffts, + sampling=sampling, + norm=norm, + real=real, + ifftshift_before=ifftshift_before, + fftshift_after=fftshift_after, + dtype=dtype, + ) + self._kwargs_fft = kwargs_fft + self._norm_kwargs = {"norm": None} # equivalent to "backward" in Numpy/Scipy + if self.norm is _FFTNorms.ORTHO: + self._norm_kwargs["norm"] = "ortho" + elif self.norm is _FFTNorms.NONE: + self._scale = np.prod(self.nffts) + elif self.norm is _FFTNorms.ONE_OVER_N: + self._scale = 1.0 / np.prod(self.nffts) + + @reshaped + def _matvec(self, x: NDArray) -> NDArray: + if self.ifftshift_before.any(): + x = mkl_backend.ifftshift(x, axes=self.axes[self.ifftshift_before]) + if not self.clinear: + x = np.real(x) + if self.real: + y = mkl_backend.rfftn( + x, s=self.nffts, axes=self.axes, **self._norm_kwargs, **self._kwargs_fft + ) + # Apply scaling to obtain a correct adjoint for this operator + y = np.swapaxes(y, -1, self.axes[-1]) + y[..., 1 : 1 + (self.nffts[-1] - 1) // 2] *= np.sqrt(2) + y = np.swapaxes(y, self.axes[-1], -1) + else: + y = mkl_backend.fftn( + x, s=self.nffts, axes=self.axes, **self._norm_kwargs, **self._kwargs_fft + ) + if self.norm is _FFTNorms.ONE_OVER_N: + y *= self._scale + if self.fftshift_after.any(): + y = mkl_backend.fftshift(y, axes=self.axes[self.fftshift_after]) + return y + + @reshaped + def _rmatvec(self, x: NDArray) -> NDArray: + if self.fftshift_after.any(): + x = mkl_backend.ifftshift(x, axes=self.axes[self.fftshift_after]) + if self.real: + # Apply scaling to obtain a correct adjoint for this operator + x = x.copy() + x = np.swapaxes(x, -1, self.axes[-1]) + x[..., 1 : 1 + (self.nffts[-1] - 1) // 2] /= np.sqrt(2) + x = np.swapaxes(x, self.axes[-1], -1) + y = mkl_backend.irfftn( + x, s=self.nffts, axes=self.axes, **self._norm_kwargs, **self._kwargs_fft + ) + else: + y = mkl_backend.ifftn( + x, s=self.nffts, axes=self.axes, **self._norm_kwargs, **self._kwargs_fft + ) + if self.norm is _FFTNorms.NONE: + y *= self._scale + for ax, nfft in zip(self.axes, self.nffts): + if nfft > self.dims[ax]: + y = np.take(y, range(self.dims[ax]), axis=ax) + if self.doifftpad: + y = np.pad(y, self.ifftpad) + if not self.clinear: + y = np.real(y) + if self.ifftshift_before.any(): + y = mkl_backend.fftshift(y, axes=self.axes[self.ifftshift_before]) + return y + + def __truediv__(self, y): + if self.norm is not _FFTNorms.ORTHO: + return self._rmatvec(y) / self._scale + return self._rmatvec(y) + + def FFTND( dims: Union[int, InputDimsLike], axes: Union[int, InputDimsLike] = (-3, -2, -1), @@ -410,7 +518,20 @@ def FFTND( for real input signals. """ - if engine == "numpy": + if engine == "mkl_fft" and mkl_fft_message is None: + f = _FFTND_mklfft( + dims=dims, + axes=axes, + nffts=nffts, + sampling=sampling, + norm=norm, + real=real, + ifftshift_before=ifftshift_before, + fftshift_after=fftshift_after, + dtype=dtype, + **kwargs_fft, + ) + elif engine == "numpy" or (engine == "mkl_fft" and mkl_fft_message is not None): f = _FFTND_numpy( dims=dims, axes=axes, @@ -437,6 +558,6 @@ def FFTND( **kwargs_fft, ) else: - raise NotImplementedError("engine must be numpy or scipy") + raise NotImplementedError("engine must be numpy, scipy or mkl_fft") f.name = name return f diff --git a/pylops/utils/deps.py b/pylops/utils/deps.py index df4a0d9e..4232cc5c 100644 --- a/pylops/utils/deps.py +++ b/pylops/utils/deps.py @@ -11,6 +11,7 @@ "sympy_enabled", "torch_enabled", "pytensor_enabled", + "mkl_fft_enabled", ] import os @@ -241,6 +242,24 @@ def pytensor_import(message: Optional[str] = None) -> str: return pytensor_message +def mkl_fft_import(message): + if mkl_fft_enabled: + try: + import_module("mkl_fft") # noqa: F401 + mkl_fft_message = None + except Exception as e: + mkl_fft_message = f"Failed to import mkl_fft (error:{e}), use numpy." + else: + mkl_fft_message = ( + "mkl_fft not available, reverting to numpy. " + "In order to be able to use " + f"{message} run " + f'"pip install mkl_fft" or ' + f'"conda install -c conda-forge mkl_fft".' + ) + return mkl_fft_message + + # Set package availability booleans # cupy and jax: the package is imported to check everything is working correctly, # if not the package is disabled. We do this here as these libraries are used as drop-in @@ -264,3 +283,4 @@ def pytensor_import(message: Optional[str] = None) -> str: sympy_enabled = util.find_spec("sympy") is not None torch_enabled = util.find_spec("torch") is not None pytensor_enabled = util.find_spec("pytensor") is not None +mkl_fft_enabled = util.find_spec("mkl_fft") is not None diff --git a/pytests/test_ffts.py b/pytests/test_ffts.py index 2c8d13c8..72a38360 100644 --- a/pytests/test_ffts.py +++ b/pytests/test_ffts.py @@ -227,11 +227,66 @@ def _choose_random_axes(ndim, n_choices=2): "dtype": np.complex128, "kwargs": {}, } # nfftnt, complex input, mkl-fft engine +par3t = { + "nt": 41, + "nx": 31, + "ny": 10, + "nfft": None, + "real": True, + "engine": "mkl_fft", + "ifftshift_before": False, + "dtype": np.float64, + "kwargs": {}, +} # nfft=nt, real input, mkl-fft engine +par4t = { + "nt": 41, + "nx": 31, + "ny": 10, + "nfft": 64, + "real": True, + "engine": "mkl_fft", + "ifftshift_before": False, + "dtype": np.float64, + "kwargs": {}, +} # nfft>nt, real input, mkl-fft engine +par5t = { + "nt": 41, + "nx": 31, + "ny": 10, + "nfft": 16, + "real": False, + "engine": "mkl_fft", + "ifftshift_before": False, + "dtype": np.complex128, + "kwargs": {}, +} # nfft=2.28.0 +mkl-fft diff --git a/requirements-doc.txt b/requirements-doc.txt index beccb183..0c275860 100644 --- a/requirements-doc.txt +++ b/requirements-doc.txt @@ -31,3 +31,4 @@ flake8 mypy pytensor>=2.28.0 pymc>=5.21.0 +mkl-fft From 5443b44833e39ac83c2cec8cb2c6a0831cfe282a Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Sun, 31 Aug 2025 00:16:56 +0530 Subject: [PATCH 02/19] Ignore flake8 checks --- pylops/signalprocessing/fft.py | 4 ++-- pylops/signalprocessing/fft2d.py | 4 ++-- pylops/signalprocessing/fftnd.py | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/pylops/signalprocessing/fft.py b/pylops/signalprocessing/fft.py index 5df2ed1f..cbbbcc7d 100644 --- a/pylops/signalprocessing/fft.py +++ b/pylops/signalprocessing/fft.py @@ -25,8 +25,8 @@ import mkl_fft.interfaces.numpy_fft as mkl_backend try: - import scipy.fft # noqa: F401 - import mkl_fft.interfaces.scipy_fft as mkl_backend + import scipy.fft # noqa: F811 + import mkl_fft.interfaces.scipy_fft as mkl_backend # noqa: F811 except ImportError: pass diff --git a/pylops/signalprocessing/fft2d.py b/pylops/signalprocessing/fft2d.py index c21d47ff..16881499 100644 --- a/pylops/signalprocessing/fft2d.py +++ b/pylops/signalprocessing/fft2d.py @@ -19,8 +19,8 @@ import mkl_fft.interfaces.numpy_fft as mkl_backend try: - import scipy.fft - import mkl_fft.interfaces.scipy_fft as mkl_backend + import scipy.fft # noqa: F811 + import mkl_fft.interfaces.scipy_fft as mkl_backend # noqa: F811 except ImportError: pass diff --git a/pylops/signalprocessing/fftnd.py b/pylops/signalprocessing/fftnd.py index 2a47e505..e57910e9 100644 --- a/pylops/signalprocessing/fftnd.py +++ b/pylops/signalprocessing/fftnd.py @@ -19,8 +19,8 @@ import mkl_fft.interfaces.numpy_fft as mkl_backend try: - import scipy.fft - import mkl_fft.interfaces.scipy_fft as mkl_backend + import scipy.fft # noqa: F401 + import mkl_fft.interfaces.scipy_fft as mkl_backend # noqa: F811 except ImportError: pass From fa06f03279b8239dd439d7dde59a52743cc0c9b5 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Sun, 31 Aug 2025 00:32:19 +0530 Subject: [PATCH 03/19] Add sys_platform to ignore mac --- requirements-dev-gpu.txt | 3 ++- requirements-dev.txt | 2 +- requirements-doc.txt | 2 +- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/requirements-dev-gpu.txt b/requirements-dev-gpu.txt index 396b52b4..94005a2c 100644 --- a/requirements-dev-gpu.txt +++ b/requirements-dev-gpu.txt @@ -24,4 +24,5 @@ isort black flake8 mypy -mkl_fft +mkl-fft; sys_platform != "darwin" + diff --git a/requirements-dev.txt b/requirements-dev.txt index cb959fc3..c1b49522 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -30,4 +30,4 @@ black flake8 mypy pytensor>=2.28.0 -mkl-fft +mkl-fft; sys_platform != "darwin" diff --git a/requirements-doc.txt b/requirements-doc.txt index 0c275860..26a4403d 100644 --- a/requirements-doc.txt +++ b/requirements-doc.txt @@ -31,4 +31,4 @@ flake8 mypy pytensor>=2.28.0 pymc>=5.21.0 -mkl-fft +mkl-fft; sys_platform != "darwin" From dad316a77af57f0d3ee64289240922006f323358 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Sun, 31 Aug 2025 01:28:09 +0530 Subject: [PATCH 04/19] Update fft code to use scipy backend(we import) and handle dtypes in code --- pylops/signalprocessing/fft.py | 35 ++++++++++++++++++++++++------ pylops/signalprocessing/fft2d.py | 35 ++++++++++++++++++++++++------ pylops/signalprocessing/fftnd.py | 37 +++++++++++++++++++++++++------- pytests/test_ffts.py | 13 ----------- 4 files changed, 85 insertions(+), 35 deletions(-) diff --git a/pylops/signalprocessing/fft.py b/pylops/signalprocessing/fft.py index cbbbcc7d..57ebcfc5 100644 --- a/pylops/signalprocessing/fft.py +++ b/pylops/signalprocessing/fft.py @@ -22,13 +22,7 @@ import pyfftw if mkl_fft_message is None: - import mkl_fft.interfaces.numpy_fft as mkl_backend - - try: - import scipy.fft # noqa: F811 - import mkl_fft.interfaces.scipy_fft as mkl_backend # noqa: F811 - except ImportError: - pass + import mkl_fft.interfaces.scipy_fft as mkl_backend logger = logging.getLogger(__name__) @@ -420,6 +414,21 @@ def __init__( dtype: DTypeLike = "complex128", **kwargs_fft, ) -> None: + if np.dtype(dtype) == np.float16: + warnings.warn( + "mkl_fft backend is unavailable with float16 dtype. Will use float32." + ) + dtype = np.float32 + elif np.dtype(dtype) == np.complex256: + warnings.warn( + "mkl_fft backend is unavailable with complex256 dtype. Will use complex128." + ) + dtype = np.complex128 + elif np.dtype(dtype) == np.float128: + warnings.warn( + "mkl_fft backend is unavailable with float128 dtype. Will use float64." + ) + dtype = np.float64 super().__init__( dims=dims, axis=axis, @@ -443,6 +452,12 @@ def __init__( @reshaped def _matvec(self, x: NDArray) -> NDArray: + if x.dtype == np.float16: + x = x.astype(np.float32) + elif x.dtype == np.float128: + x = x.astype(np.float64) + elif x.dtype == np.complex256: + x = x.astype(np.complex128) if self.ifftshift_before: x = mkl_backend.ifftshift(x, axes=self.axis) if not self.clinear: @@ -462,6 +477,12 @@ def _matvec(self, x: NDArray) -> NDArray: @reshaped def _rmatvec(self, x: NDArray) -> NDArray: + if x.dtype == np.float16: + x = x.astype(np.float32) + elif x.dtype == np.float128: + x = x.astype(np.float64) + elif x.dtype == np.complex256: + x = x.astype(np.complex128) if self.fftshift_after: x = mkl_backend.ifftshift(x, axes=self.axis) if self.real: diff --git a/pylops/signalprocessing/fft2d.py b/pylops/signalprocessing/fft2d.py index 16881499..915ceffc 100644 --- a/pylops/signalprocessing/fft2d.py +++ b/pylops/signalprocessing/fft2d.py @@ -16,13 +16,7 @@ mkl_fft_message = deps.mkl_fft_import("the mkl fft module") if mkl_fft_message is None: - import mkl_fft.interfaces.numpy_fft as mkl_backend - - try: - import scipy.fft # noqa: F811 - import mkl_fft.interfaces.scipy_fft as mkl_backend # noqa: F811 - except ImportError: - pass + import mkl_fft.interfaces.scipy_fft as mkl_backend class _FFT2D_numpy(_BaseFFTND): @@ -263,6 +257,21 @@ def __init__( dtype: DTypeLike = "complex128", **kwargs_fft, ) -> None: + if np.dtype(dtype) == np.float16: + warnings.warn( + "mkl_fft backend is unavailable with float16 dtype. Will use float32." + ) + dtype = np.float32 + elif np.dtype(dtype) == np.complex256: + warnings.warn( + "mkl_fft backend is unavailable with complex256 dtype. Will use complex128." + ) + dtype = np.complex128 + elif np.dtype(dtype) == np.float128: + warnings.warn( + "mkl_fft backend is unavailable with float128 dtype. Will use float64." + ) + dtype = np.float64 super().__init__( dims=dims, axes=axes, @@ -296,6 +305,12 @@ def __init__( @reshaped def _matvec(self, x): + if x.dtype == np.float16: + x = x.astype(np.float32) + elif x.dtype == np.float128: + x = x.astype(np.float64) + elif x.dtype == np.complex256: + x = x.astype(np.complex128) if self.ifftshift_before.any(): x = mkl_backend.ifftshift(x, axes=self.axes[self.ifftshift_before]) if not self.clinear: @@ -320,6 +335,12 @@ def _matvec(self, x): @reshaped def _rmatvec(self, x): + if x.dtype == np.float16: + x = x.astype(np.float32) + elif x.dtype == np.float128: + x = x.astype(np.float64) + elif x.dtype == np.complex256: + x = x.astype(np.complex128) if self.fftshift_after.any(): x = mkl_backend.ifftshift(x, axes=self.axes[self.fftshift_after]) if self.real: diff --git a/pylops/signalprocessing/fftnd.py b/pylops/signalprocessing/fftnd.py index e57910e9..06f95275 100644 --- a/pylops/signalprocessing/fftnd.py +++ b/pylops/signalprocessing/fftnd.py @@ -16,13 +16,7 @@ mkl_fft_message = deps.mkl_fft_import("the mkl fft module") if mkl_fft_message is None: - import mkl_fft.interfaces.numpy_fft as mkl_backend - - try: - import scipy.fft # noqa: F401 - import mkl_fft.interfaces.scipy_fft as mkl_backend # noqa: F811 - except ImportError: - pass + import mkl_fft.interfaces.scipy_fft as mkl_backend class _FFTND_numpy(_BaseFFTND): @@ -244,6 +238,21 @@ def __init__( dtype: DTypeLike = "complex128", **kwargs_fft, ) -> None: + if np.dtype(dtype) == np.float16: + warnings.warn( + "mkl_fft backend is unavailable with float16 dtype. Will use float32." + ) + dtype = np.float32 + elif np.dtype(dtype) == np.complex256: + warnings.warn( + "mkl_fft backend is unavailable with complex256 dtype. Will use complex128." + ) + dtype = np.complex128 + elif np.dtype(dtype) == np.float128: + warnings.warn( + "mkl_fft backend is unavailable with float128 dtype. Will use float64." + ) + dtype = np.float64 super().__init__( dims=dims, axes=axes, @@ -266,6 +275,12 @@ def __init__( @reshaped def _matvec(self, x: NDArray) -> NDArray: + if x.dtype == np.float16: + x = x.astype(np.float32) + elif x.dtype == np.float128: + x = x.astype(np.float64) + elif x.dtype == np.complex256: + x = x.astype(np.complex128) if self.ifftshift_before.any(): x = mkl_backend.ifftshift(x, axes=self.axes[self.ifftshift_before]) if not self.clinear: @@ -290,6 +305,12 @@ def _matvec(self, x: NDArray) -> NDArray: @reshaped def _rmatvec(self, x: NDArray) -> NDArray: + if x.dtype == np.float16: + x = x.astype(np.float32) + elif x.dtype == np.float128: + x = x.astype(np.float64) + elif x.dtype == np.complex256: + x = x.astype(np.complex128) if self.fftshift_after.any(): x = mkl_backend.ifftshift(x, axes=self.axes[self.fftshift_after]) if self.real: @@ -318,7 +339,7 @@ def _rmatvec(self, x: NDArray) -> NDArray: y = mkl_backend.fftshift(y, axes=self.axes[self.ifftshift_before]) return y - def __truediv__(self, y): + def __truediv__(self, y: npt.ArrayLike) -> npt.ArrayLike: if self.norm is not _FFTNorms.ORTHO: return self._rmatvec(y) / self._scale return self._rmatvec(y) diff --git a/pytests/test_ffts.py b/pytests/test_ffts.py index 72a38360..f6d7c5bf 100644 --- a/pytests/test_ffts.py +++ b/pytests/test_ffts.py @@ -314,11 +314,9 @@ def test_unknown_engine(par): engine=["numpy", "fftw", "scipy", "mkl_fft"], ) # Generate all combinations of the above parameters -# Intel mkl fft raises NotImplemented Error on np.float16, npfloat128 and np.complex256 pars_fft_small_real = [ dict(zip(par_lists_fft_small_real.keys(), value)) for value in itertools.product(*par_lists_fft_small_real.values()) - if (value[3] != "mkl_fft") or (value[0][0] not in (np.float16, np.float128)) ] @@ -392,7 +390,6 @@ def test_FFT_small_real(par): pars_fft_random_real = [ dict(zip(par_lists_fft_random_real.keys(), value)) for value in itertools.product(*par_lists_fft_random_real.values()) - if (value[3] != "mkl_fft") or (value[1][0] not in (np.float16, np.float128)) ] @@ -453,7 +450,6 @@ def test_FFT_random_real(par): pars_fft_small_cpx = [ dict(zip(par_lists_fft_small_cpx.keys(), value)) for value in itertools.product(*par_lists_fft_small_cpx.values()) - if (value[4] != "mkl_fft") or (value[0][0] not in (np.complex256,)) ] @@ -526,8 +522,6 @@ def test_FFT_small_complex(par): pars_fft_random_cpx = [ dict(zip(par_lists_fft_random_cpx.keys(), value)) for value in itertools.product(*par_lists_fft_random_cpx.values()) - if (value[4] != "mkl_fft") - or (value[1][0] not in (np.float16, np.float128, np.complex256)) ] @@ -611,7 +605,6 @@ def test_FFT_random_complex(par): pars_fft2d_random_real = [ dict(zip(par_lists_fft2d_random_real.keys(), value)) for value in itertools.product(*par_lists_fft2d_random_real.values()) - if (value[3] != "mkl_fft") or (value[1][0] not in (np.float16, np.float128)) ] @@ -677,8 +670,6 @@ def test_FFT2D_random_real(par): pars_fft2d_random_cpx = [ dict(zip(par_lists_fft2d_random_cpx.keys(), value)) for value in itertools.product(*par_lists_fft2d_random_cpx.values()) - if (value[4] != "mkl_fft") - or (value[1][0] not in (np.float16, np.float128, np.complex256)) ] @@ -760,7 +751,6 @@ def test_FFT2D_random_complex(par): pars_fftnd_random_real = [ dict(zip(par_lists_fftnd_random_real.keys(), value)) for value in itertools.product(*par_lists_fftnd_random_real.values()) - if (value[2] != "mkl_fft") or (value[1][0] not in (np.float16, np.float128)) ] @@ -822,8 +812,6 @@ def test_FFTND_random_real(par): pars_fftnd_random_cpx = [ dict(zip(par_lists_fftnd_random_cpx.keys(), value)) for value in itertools.product(*par_lists_fftnd_random_cpx.values()) - if (value[2] != "mkl_fft") - or (value[1][0] not in (np.float16, np.float128, np.complex256)) ] @@ -900,7 +888,6 @@ def test_FFTND_random_complex(par): pars_fft2dnd_small_cpx = [ dict(zip(par_lists_fft2dnd_small_cpx.keys(), value)) for value in itertools.product(*par_lists_fft2dnd_small_cpx.values()) - if (value[2] != "mkl_fft") or (value[0][0] not in (np.complex256,)) ] From 129f485e30a6e21cc393021d518b83a3c3e3ca6a Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Sun, 31 Aug 2025 02:10:36 +0530 Subject: [PATCH 05/19] Remove mkl_fft from gpu reqs --- environment-dev-gpu.yml | 1 - requirements-dev-gpu.txt | 2 -- 2 files changed, 3 deletions(-) diff --git a/environment-dev-gpu.yml b/environment-dev-gpu.yml index 1ae98691..8436d228 100644 --- a/environment-dev-gpu.yml +++ b/environment-dev-gpu.yml @@ -21,7 +21,6 @@ dependencies: - autopep8 - isort - black - - mkl_fft - pip: - torch - pytest-runner diff --git a/requirements-dev-gpu.txt b/requirements-dev-gpu.txt index 94005a2c..22ed24a9 100644 --- a/requirements-dev-gpu.txt +++ b/requirements-dev-gpu.txt @@ -24,5 +24,3 @@ isort black flake8 mypy -mkl-fft; sys_platform != "darwin" - From 1fe6b34937179a77402504161bfa7a3345b97032 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Mon, 1 Sep 2025 14:13:12 +0530 Subject: [PATCH 06/19] Update index url and check for macos Add check to skip mkl_fft on darwin Minor update --- pylops/signalprocessing/fft.py | 29 +---------------------------- pylops/signalprocessing/fft2d.py | 29 +---------------------------- pylops/signalprocessing/fftnd.py | 29 +---------------------------- pytests/test_ffts.py | 31 +++++++++++++++++++++++++++++++ requirements-dev.txt | 4 +++- requirements-doc.txt | 4 +++- 6 files changed, 40 insertions(+), 86 deletions(-) diff --git a/pylops/signalprocessing/fft.py b/pylops/signalprocessing/fft.py index 57ebcfc5..7746e05a 100644 --- a/pylops/signalprocessing/fft.py +++ b/pylops/signalprocessing/fft.py @@ -22,7 +22,7 @@ import pyfftw if mkl_fft_message is None: - import mkl_fft.interfaces.scipy_fft as mkl_backend + import mkl_fft.interfaces.numpy_fft as mkl_backend logger = logging.getLogger(__name__) @@ -414,21 +414,6 @@ def __init__( dtype: DTypeLike = "complex128", **kwargs_fft, ) -> None: - if np.dtype(dtype) == np.float16: - warnings.warn( - "mkl_fft backend is unavailable with float16 dtype. Will use float32." - ) - dtype = np.float32 - elif np.dtype(dtype) == np.complex256: - warnings.warn( - "mkl_fft backend is unavailable with complex256 dtype. Will use complex128." - ) - dtype = np.complex128 - elif np.dtype(dtype) == np.float128: - warnings.warn( - "mkl_fft backend is unavailable with float128 dtype. Will use float64." - ) - dtype = np.float64 super().__init__( dims=dims, axis=axis, @@ -452,12 +437,6 @@ def __init__( @reshaped def _matvec(self, x: NDArray) -> NDArray: - if x.dtype == np.float16: - x = x.astype(np.float32) - elif x.dtype == np.float128: - x = x.astype(np.float64) - elif x.dtype == np.complex256: - x = x.astype(np.complex128) if self.ifftshift_before: x = mkl_backend.ifftshift(x, axes=self.axis) if not self.clinear: @@ -477,12 +456,6 @@ def _matvec(self, x: NDArray) -> NDArray: @reshaped def _rmatvec(self, x: NDArray) -> NDArray: - if x.dtype == np.float16: - x = x.astype(np.float32) - elif x.dtype == np.float128: - x = x.astype(np.float64) - elif x.dtype == np.complex256: - x = x.astype(np.complex128) if self.fftshift_after: x = mkl_backend.ifftshift(x, axes=self.axis) if self.real: diff --git a/pylops/signalprocessing/fft2d.py b/pylops/signalprocessing/fft2d.py index 915ceffc..f301b4c3 100644 --- a/pylops/signalprocessing/fft2d.py +++ b/pylops/signalprocessing/fft2d.py @@ -16,7 +16,7 @@ mkl_fft_message = deps.mkl_fft_import("the mkl fft module") if mkl_fft_message is None: - import mkl_fft.interfaces.scipy_fft as mkl_backend + import mkl_fft.interfaces.numpy_fft as mkl_backend class _FFT2D_numpy(_BaseFFTND): @@ -257,21 +257,6 @@ def __init__( dtype: DTypeLike = "complex128", **kwargs_fft, ) -> None: - if np.dtype(dtype) == np.float16: - warnings.warn( - "mkl_fft backend is unavailable with float16 dtype. Will use float32." - ) - dtype = np.float32 - elif np.dtype(dtype) == np.complex256: - warnings.warn( - "mkl_fft backend is unavailable with complex256 dtype. Will use complex128." - ) - dtype = np.complex128 - elif np.dtype(dtype) == np.float128: - warnings.warn( - "mkl_fft backend is unavailable with float128 dtype. Will use float64." - ) - dtype = np.float64 super().__init__( dims=dims, axes=axes, @@ -305,12 +290,6 @@ def __init__( @reshaped def _matvec(self, x): - if x.dtype == np.float16: - x = x.astype(np.float32) - elif x.dtype == np.float128: - x = x.astype(np.float64) - elif x.dtype == np.complex256: - x = x.astype(np.complex128) if self.ifftshift_before.any(): x = mkl_backend.ifftshift(x, axes=self.axes[self.ifftshift_before]) if not self.clinear: @@ -335,12 +314,6 @@ def _matvec(self, x): @reshaped def _rmatvec(self, x): - if x.dtype == np.float16: - x = x.astype(np.float32) - elif x.dtype == np.float128: - x = x.astype(np.float64) - elif x.dtype == np.complex256: - x = x.astype(np.complex128) if self.fftshift_after.any(): x = mkl_backend.ifftshift(x, axes=self.axes[self.fftshift_after]) if self.real: diff --git a/pylops/signalprocessing/fftnd.py b/pylops/signalprocessing/fftnd.py index 06f95275..4549a074 100644 --- a/pylops/signalprocessing/fftnd.py +++ b/pylops/signalprocessing/fftnd.py @@ -16,7 +16,7 @@ mkl_fft_message = deps.mkl_fft_import("the mkl fft module") if mkl_fft_message is None: - import mkl_fft.interfaces.scipy_fft as mkl_backend + import mkl_fft.interfaces.numpy_fft as mkl_backend class _FFTND_numpy(_BaseFFTND): @@ -238,21 +238,6 @@ def __init__( dtype: DTypeLike = "complex128", **kwargs_fft, ) -> None: - if np.dtype(dtype) == np.float16: - warnings.warn( - "mkl_fft backend is unavailable with float16 dtype. Will use float32." - ) - dtype = np.float32 - elif np.dtype(dtype) == np.complex256: - warnings.warn( - "mkl_fft backend is unavailable with complex256 dtype. Will use complex128." - ) - dtype = np.complex128 - elif np.dtype(dtype) == np.float128: - warnings.warn( - "mkl_fft backend is unavailable with float128 dtype. Will use float64." - ) - dtype = np.float64 super().__init__( dims=dims, axes=axes, @@ -275,12 +260,6 @@ def __init__( @reshaped def _matvec(self, x: NDArray) -> NDArray: - if x.dtype == np.float16: - x = x.astype(np.float32) - elif x.dtype == np.float128: - x = x.astype(np.float64) - elif x.dtype == np.complex256: - x = x.astype(np.complex128) if self.ifftshift_before.any(): x = mkl_backend.ifftshift(x, axes=self.axes[self.ifftshift_before]) if not self.clinear: @@ -305,12 +284,6 @@ def _matvec(self, x: NDArray) -> NDArray: @reshaped def _rmatvec(self, x: NDArray) -> NDArray: - if x.dtype == np.float16: - x = x.astype(np.float32) - elif x.dtype == np.float128: - x = x.astype(np.float64) - elif x.dtype == np.complex256: - x = x.astype(np.complex128) if self.fftshift_after.any(): x = mkl_backend.ifftshift(x, axes=self.axes[self.fftshift_after]) if self.real: diff --git a/pytests/test_ffts.py b/pytests/test_ffts.py index f6d7c5bf..5fdc6066 100644 --- a/pytests/test_ffts.py +++ b/pytests/test_ffts.py @@ -1,5 +1,6 @@ import itertools import os +import sys if int(os.environ.get("TEST_CUPY_PYLOPS", 0)): import cupy as np @@ -322,6 +323,8 @@ def test_unknown_engine(par): @pytest.mark.parametrize("par", pars_fft_small_real) def test_FFT_small_real(par): + if par["engine"] == "mkl_fft" and sys.platform == "darwin": + pytest.skip("mkl_fft not supported on macOS") np.random.seed(5) if backend == "numpy" or (backend == "cupy" and par["engine"] == "numpy"): @@ -399,6 +402,8 @@ def test_FFT_small_real(par): ) @pytest.mark.parametrize("par", pars_fft_random_real) def test_FFT_random_real(par): + if par["engine"] == "mkl_fft" and sys.platform == "darwin": + pytest.skip("mkl_fft not supported on macOS") np.random.seed(5) shape = par["shape"] @@ -455,6 +460,8 @@ def test_FFT_random_real(par): @pytest.mark.parametrize("par", pars_fft_small_cpx) def test_FFT_small_complex(par): + if par["engine"] == "mkl_fft" and sys.platform == "darwin": + pytest.skip("mkl_fft not supported on macOS") np.random.seed(5) dtype, decimal = par["dtype_precision"] norm = par["norm"] @@ -527,6 +534,8 @@ def test_FFT_small_complex(par): @pytest.mark.parametrize("par", pars_fft_random_cpx) def test_FFT_random_complex(par): + if par["engine"] == "mkl_fft" and sys.platform == "darwin": + pytest.skip("mkl_fft not supported on macOS") np.random.seed(5) if backend == "numpy" or (backend == "cupy" and par["engine"] == "numpy"): shape = par["shape"] @@ -614,6 +623,8 @@ def test_FFT_random_complex(par): ) @pytest.mark.parametrize("par", pars_fft2d_random_real) def test_FFT2D_random_real(par): + if par["engine"] == "mkl_fft" and sys.platform == "darwin": + pytest.skip("mkl_fft not supported on macOS") np.random.seed(5) if backend == "numpy" or (backend == "cupy" and par["engine"] == "numpy"): shape = par["shape"] @@ -675,6 +686,8 @@ def test_FFT2D_random_real(par): @pytest.mark.parametrize("par", pars_fft2d_random_cpx) def test_FFT2D_random_complex(par): + if par["engine"] == "mkl_fft" and sys.platform == "darwin": + pytest.skip("mkl_fft not supported on macOS") np.random.seed(5) if backend == "numpy" or (backend == "cupy" and par["engine"] == "numpy"): shape = par["shape"] @@ -756,6 +769,8 @@ def test_FFT2D_random_complex(par): @pytest.mark.parametrize("par", pars_fftnd_random_real) def test_FFTND_random_real(par): + if par["engine"] == "mkl_fft" and sys.platform == "darwin": + pytest.skip("mkl_fft not supported on macOS") np.random.seed(5) if backend == "numpy" or (backend == "cupy" and par["engine"] == "numpy"): shape = par["shape"] @@ -817,6 +832,8 @@ def test_FFTND_random_real(par): @pytest.mark.parametrize("par", pars_fftnd_random_cpx) def test_FFTND_random_complex(par): + if par["engine"] == "mkl_fft" and sys.platform == "darwin": + pytest.skip("mkl_fft not supported on macOS") np.random.seed(5) shape = par["shape"] dtype, decimal = par["dtype_precision"] @@ -893,6 +910,8 @@ def test_FFTND_random_complex(par): @pytest.mark.parametrize("par", pars_fft2dnd_small_cpx) def test_FFT2D_small_complex(par): + if par["engine"] == "mkl_fft" and sys.platform == "darwin": + pytest.skip("mkl_fft not supported on macOS") np.random.seed(5) dtype, decimal = par["dtype_precision"] norm = par["norm"] @@ -943,6 +962,8 @@ def test_FFT2D_small_complex(par): @pytest.mark.parametrize("par", pars_fft2dnd_small_cpx) def test_FFTND_small_complex(par): + if par["engine"] == "mkl_fft" and sys.platform == "darwin": + pytest.skip("mkl_fft not supported on macOS") np.random.seed(5) dtype, decimal = par["dtype_precision"] norm = par["norm"] @@ -1018,6 +1039,8 @@ def test_FFTND_small_complex(par): ], ) def test_FFT_1dsignal(par): + if par["engine"] == "mkl_fft" and sys.platform == "darwin": + pytest.skip("mkl_fft not supported on macOS") np.random.seed(5) """Dot-test and inversion for FFT operator for 1d signal""" decimal = 3 if np.real(np.ones(1, par["dtype"])).dtype == np.float32 else 8 @@ -1135,6 +1158,8 @@ def test_FFT_2dsignal(par): """Dot-test and inversion for fft operator for 2d signal (fft on single dimension) """ + if par["engine"] == "mkl_fft" and sys.platform == "darwin": + pytest.skip("mkl_fft not supported on macOS") np.random.seed(5) decimal = 3 if np.real(np.ones(1, par["dtype"])).dtype == np.float32 else 8 @@ -1350,6 +1375,8 @@ def test_FFT_3dsignal(par): """Dot-test and inversion for fft operator for 3d signal (fft on single dimension) """ + if par["engine"] == "mkl_fft" and sys.platform == "darwin": + pytest.skip("mkl_fft not supported on macOS") np.random.seed(5) decimal = 3 if np.real(np.ones(1, par["dtype"])).dtype == np.float32 else 8 @@ -1576,6 +1603,8 @@ def test_FFT_3dsignal(par): ) def test_FFT2D(par): """Dot-test and inversion for FFT2D operator for 2d signal""" + if par["engine"] == "mkl_fft" and sys.platform == "darwin": + pytest.skip("mkl_fft not supported on macOS") np.random.seed(5) decimal = 3 if np.real(np.ones(1, par["dtype"])).dtype == np.float32 else 8 @@ -1711,6 +1740,8 @@ def test_FFT2D(par): ) def test_FFT3D(par): """Dot-test and inversion for FFTND operator for 3d signal""" + if par["engine"] == "mkl_fft" and sys.platform == "darwin": + pytest.skip("mkl_fft not supported on macOS") np.random.seed(5) decimal = 3 if np.real(np.ones(1, par["dtype"])).dtype == np.float32 else 8 diff --git a/requirements-dev.txt b/requirements-dev.txt index c1b49522..9049bf4a 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -1,3 +1,6 @@ +--index-url https://software.repos.intel.com/python/pypi +--extra-index-url https://pypi.org/simple +mkl_fft; sys_platform != "darwin" numpy>=2.0.0 scipy>=1.13.0 jax @@ -30,4 +33,3 @@ black flake8 mypy pytensor>=2.28.0 -mkl-fft; sys_platform != "darwin" diff --git a/requirements-doc.txt b/requirements-doc.txt index 26a4403d..cf520658 100644 --- a/requirements-doc.txt +++ b/requirements-doc.txt @@ -1,3 +1,6 @@ +--index-url https://software.repos.intel.com/python/pypi +--extra-index-url https://pypi.org/simple +mkl_fft; sys_platform != "darwin" numpy>=2.0.0 scipy>=1.13.0 jax @@ -31,4 +34,3 @@ flake8 mypy pytensor>=2.28.0 pymc>=5.21.0 -mkl-fft; sys_platform != "darwin" From 3ad2b57d94f87f5a4b785c88a1b1946a86b83475 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Mon, 1 Sep 2025 15:12:32 +0530 Subject: [PATCH 07/19] Update installation docs --- docs/source/installation.rst | 25 ++++++++++++++++++++----- environment-dev-arm.yml | 1 - environment-dev.yml | 1 + 3 files changed, 21 insertions(+), 6 deletions(-) diff --git a/docs/source/installation.rst b/docs/source/installation.rst index 9f148a12..40c15048 100755 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -358,11 +358,11 @@ Install it via ``pip`` with FFTW ---- -Three different "engines" are provided by the :py:class:`pylops.signalprocessing.FFT` operator: -``engine="numpy"`` (default), ``engine="scipy"`` and ``engine="fftw"``. +Four different "engines" are provided by the :py:class:`pylops.signalprocessing.FFT` operator: +``engine="numpy"`` (default), ``engine="scipy"``, ``engine="fftw"`` and ``engine="mkl_fft"``. The first two engines are part of the required PyLops dependencies. -The latter implements the well-known `FFTW `_ +The third implements the well-known `FFTW `_ via the Python wrapper :py:class:`pyfftw.FFTW`. While this optimized FFT tends to outperform the other two in many cases, it is not included by default. To use this library, install it manually either via ``conda``: @@ -381,9 +381,24 @@ or via pip: FFTW is only available for :py:class:`pylops.signalprocessing.FFT`, not :py:class:`pylops.signalprocessing.FFT2D` or :py:class:`pylops.signalprocessing.FFTND`. -.. warning:: - Intel MKL FFT is not supported. +The fourth implements Intel MKL FFT via the Python interface `mkl_fft `_. +This provides access to Intel’s oneMKL Fourier Transform routines, enabling efficient FFT computations with performance +close to native C/Intel® oneMKL + +To use this library, you can install it using ``conda``: + +.. code-block:: bash + + >> conda install --channel https://software.repos.intel.com/python/conda --channel conda-forge mkl_fft +or via pip: + +.. code-block:: bash + + >> pip install --index-url https://software.repos.intel.com/python/pypi --extra-index-url https://pypi.org/simple mkl_fft + +.. note:: + `mkl_fft` is not supported on macOS Numba ----- diff --git a/environment-dev-arm.yml b/environment-dev-arm.yml index ae485474..49ae8824 100644 --- a/environment-dev-arm.yml +++ b/environment-dev-arm.yml @@ -23,7 +23,6 @@ dependencies: - autopep8 - isort - black - - mkl_fft - pip: - torch - devito diff --git a/environment-dev.yml b/environment-dev.yml index 8ca130c0..87fbffd0 100644 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -1,5 +1,6 @@ name: pylops channels: + - https://software.repos.intel.com/python/conda - defaults - conda-forge - numba From 7f38d4519dab745facc8244021358806aae7da57 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Sat, 6 Sep 2025 17:14:31 +0530 Subject: [PATCH 08/19] Add GA for build-mkl Add new requirements Update code and tests Minor change Fix-Pass login shell as default Add support for python3.10 --- .github/workflows/build-mkl.yaml | 41 +++++++++++++++++++++ Makefile | 11 +++++- docs/source/installation.rst | 4 +- environment-dev-intel-mkl.yml | 45 +++++++++++++++++++++++ environment-dev.yml | 2 - pylops/signalprocessing/fft.py | 15 +++++--- pylops/signalprocessing/fft2d.py | 16 +++++--- pylops/signalprocessing/fftnd.py | 16 +++++--- pytests/test_ffts.py | 63 ++++++++++++++++---------------- requirements-dev.txt | 3 -- requirements-intel-mkl.txt | 4 ++ 11 files changed, 164 insertions(+), 56 deletions(-) create mode 100644 .github/workflows/build-mkl.yaml create mode 100644 environment-dev-intel-mkl.yml create mode 100644 requirements-intel-mkl.txt diff --git a/.github/workflows/build-mkl.yaml b/.github/workflows/build-mkl.yaml new file mode 100644 index 00000000..c086d826 --- /dev/null +++ b/.github/workflows/build-mkl.yaml @@ -0,0 +1,41 @@ +name: PyLops Testing with Intel oneAPI Math Kernel Library(oneMKL) + +on: [push, pull_request] + +jobs: + build: + strategy: + matrix: + platform: [ubuntu-latest] + python-version: ["3.10", "3.11", "3.12"] + + runs-on: ${{ matrix.platform }} + defaults: + run: + shell: bash -l {0} + steps: + - uses: actions/checkout@v4 + - name: Get history and tags for SCM versioning to work + run: | + git fetch --prune --unshallow + git fetch --depth=1 origin +refs/tags/*:refs/tags/* + - uses: conda-incubator/setup-miniconda@v3.2.0 + with: + use-mamba: true + channels: https://software.repos.intel.com/python/conda, conda-forge + conda-remove-defaults: true + python-version: ${{ matrix.python-version }} + activate-environment: mkl-test-env + - name: Install dependencies + run: | + conda install -y pyfftw + pip install -r requirements-intel-mkl.txt + pip install -r requirements-dev.txt + pip install -r requirements-torch.txt + - name: Install pylops + run: | + python -m setuptools_scm + pip install . + - name: Tests with pytest + run: | + pytest diff --git a/Makefile b/Makefile index 06c36ca5..3df51181 100644 --- a/Makefile +++ b/Makefile @@ -1,7 +1,7 @@ PIP := $(shell command -v pip3 2> /dev/null || command which pip 2> /dev/null) PYTHON := $(shell command -v python3 2> /dev/null || command which python 2> /dev/null) -.PHONY: install dev-install dev-install_gpu install_conda dev-install_conda dev-install_conda_arm tests tests_cpu_ongpu tests_gpu doc docupdate servedoc lint typeannot coverage +.PHONY: install dev-install dev-install_intel_mkl dev-install_gpu install_conda dev-install_conda dev-install_conda_intel_mkl dev-install_conda_arm tests tests_cpu_ongpu tests_gpu doc docupdate servedoc lint typeannot coverage pipcheck: ifndef PIP @@ -24,6 +24,12 @@ dev-install: $(PIP) install -r requirements-dev.txt &&\ $(PIP) install -r requirements-torch.txt && $(PIP) install -e . +dev-install_intel_mkl: + make pipcheck + $(PIP) install -r requirements-intel-mkl.txt &&\ + $(PIP) install -r requirements-dev.txt &&\ + $(PIP) install -r requirements-torch.txt && $(PIP) install -e . + dev-install_gpu: make pipcheck $(PIP) install -r requirements-dev-gpu.txt &&\ @@ -35,6 +41,9 @@ install_conda: dev-install_conda: conda env create -f environment-dev.yml && source ${CONDA_PREFIX}/etc/profile.d/conda.sh && conda activate pylops && pip install -e . +dev-install_conda_intel_mkl: + conda env create -f environment-dev-intel-mkl.yml && source ${CONDA_PREFIX}/etc/profile.d/conda.sh && conda activate pylops && pip install -e . + dev-install_conda_arm: conda env create -f environment-dev-arm.yml && source ${CONDA_PREFIX}/etc/profile.d/conda.sh && conda activate pylops && pip install -e . diff --git a/docs/source/installation.rst b/docs/source/installation.rst index 40c15048..1776bc5d 100755 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -356,8 +356,8 @@ Install it via ``pip`` with >> pip install devito -FFTW ----- +FFTW and MKL-FFT +---------------- Four different "engines" are provided by the :py:class:`pylops.signalprocessing.FFT` operator: ``engine="numpy"`` (default), ``engine="scipy"``, ``engine="fftw"`` and ``engine="mkl_fft"``. diff --git a/environment-dev-intel-mkl.yml b/environment-dev-intel-mkl.yml new file mode 100644 index 00000000..e5bf97e2 --- /dev/null +++ b/environment-dev-intel-mkl.yml @@ -0,0 +1,45 @@ +name: pylops +channels: + - https://software.repos.intel.com/python/conda + - conda-forge + - defaults + - numba +dependencies: + - python==3.12.8 + - pip + - numpy>=2.0.0 + - scipy>=1.13.0 + - pyfftw + - pywavelets + - sympy + - pymc>=5 + - pytensor + - matplotlib + - ipython + - pytest + - Sphinx + - numpydoc + - numba + - icc_rt + - pre-commit + - autopep8 + - isort + - black + - mkl_fft + - pip: + - torch + - devito + # - dtcwt (until numpy>=2.0.0 is supported) + - scikit-fmm + - spgl1 + - jax + - pytest-runner + - setuptools_scm + - pydata-sphinx-theme + - pooch + - sphinx-gallery + - nbsphinx + - sphinxemoji + - image + - flake8 + - mypy diff --git a/environment-dev.yml b/environment-dev.yml index 87fbffd0..fccd5663 100644 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -1,6 +1,5 @@ name: pylops channels: - - https://software.repos.intel.com/python/conda - defaults - conda-forge - numba @@ -25,7 +24,6 @@ dependencies: - autopep8 - isort - black - - mkl_fft - pip: - torch - devito diff --git a/pylops/signalprocessing/fft.py b/pylops/signalprocessing/fft.py index 7746e05a..c5d4a15b 100644 --- a/pylops/signalprocessing/fft.py +++ b/pylops/signalprocessing/fft.py @@ -22,7 +22,8 @@ import pyfftw if mkl_fft_message is None: - import mkl_fft.interfaces.numpy_fft as mkl_backend + import mkl_fft.interfaces.scipy_fft as mkl_backend + from mkl_fft.interfaces import _float_utils logger = logging.getLogger(__name__) @@ -437,8 +438,10 @@ def __init__( @reshaped def _matvec(self, x: NDArray) -> NDArray: + x = _float_utils._downcast_float128_array(x) + x = _float_utils._upcast_float16_array(x) if self.ifftshift_before: - x = mkl_backend.ifftshift(x, axes=self.axis) + x = scipy.fft.ifftshift(x, axes=self.axis) if not self.clinear: x = np.real(x) if self.real: @@ -451,13 +454,15 @@ def _matvec(self, x: NDArray) -> NDArray: if self.norm is _FFTNorms.ONE_OVER_N: y *= self._scale if self.fftshift_after: - y = mkl_backend.fftshift(y, axes=self.axis) + y = scipy.fft.fftshift(y, axes=self.axis) return y @reshaped def _rmatvec(self, x: NDArray) -> NDArray: + x = _float_utils._downcast_float128_array(x) + x = _float_utils._upcast_float16_array(x) if self.fftshift_after: - x = mkl_backend.ifftshift(x, axes=self.axis) + x = scipy.fft.ifftshift(x, axes=self.axis) if self.real: x = x.copy() x = np.swapaxes(x, -1, self.axis) @@ -477,7 +482,7 @@ def _rmatvec(self, x: NDArray) -> NDArray: if not self.clinear: y = np.real(y) if self.ifftshift_before: - y = mkl_backend.fftshift(y, axes=self.axis) + y = scipy.fft.fftshift(y, axes=self.axis) return y def __truediv__(self, y): diff --git a/pylops/signalprocessing/fft2d.py b/pylops/signalprocessing/fft2d.py index f301b4c3..3f79364f 100644 --- a/pylops/signalprocessing/fft2d.py +++ b/pylops/signalprocessing/fft2d.py @@ -16,7 +16,8 @@ mkl_fft_message = deps.mkl_fft_import("the mkl fft module") if mkl_fft_message is None: - import mkl_fft.interfaces.numpy_fft as mkl_backend + import mkl_fft.interfaces.scipy_fft as mkl_backend + from mkl_fft.interfaces import _float_utils class _FFT2D_numpy(_BaseFFTND): @@ -290,8 +291,10 @@ def __init__( @reshaped def _matvec(self, x): + x = _float_utils._downcast_float128_array(x) + x = _float_utils._upcast_float16_array(x) if self.ifftshift_before.any(): - x = mkl_backend.ifftshift(x, axes=self.axes[self.ifftshift_before]) + x = scipy.fft.ifftshift(x, axes=self.axes[self.ifftshift_before]) if not self.clinear: x = np.real(x) if self.real: @@ -309,13 +312,15 @@ def _matvec(self, x): y *= self._scale y = y.astype(self.cdtype) if self.fftshift_after.any(): - y = mkl_backend.fftshift(y, axes=self.axes[self.fftshift_after]) + y = scipy.fft.fftshift(y, axes=self.axes[self.fftshift_after]) return y @reshaped def _rmatvec(self, x): + x = _float_utils._downcast_float128_array(x) + x = _float_utils._upcast_float16_array(x) if self.fftshift_after.any(): - x = mkl_backend.ifftshift(x, axes=self.axes[self.fftshift_after]) + x = scipy.fft.ifftshift(x, axes=self.axes[self.fftshift_after]) if self.real: x = x.copy() x = np.swapaxes(x, -1, self.axes[-1]) @@ -330,7 +335,6 @@ def _rmatvec(self, x): ) if self.norm is _FFTNorms.NONE: y *= self._scale - print(y.shape, self.dims[self.axes[0]]) if self.nffts[0] > self.dims[self.axes[0]]: y = np.take(y, np.arange(self.dims[self.axes[0]]), axis=self.axes[0]) if self.nffts[1] > self.dims[self.axes[1]]: @@ -341,7 +345,7 @@ def _rmatvec(self, x): y = np.real(y) y = y.astype(self.rdtype) if self.ifftshift_before.any(): - y = np.fft.fftshift(y, axes=self.axes[self.ifftshift_before]) + y = scipy.fft.fftshift(y, axes=self.axes[self.ifftshift_before]) return y def __truediv__(self, y): diff --git a/pylops/signalprocessing/fftnd.py b/pylops/signalprocessing/fftnd.py index 4549a074..2754cf8d 100644 --- a/pylops/signalprocessing/fftnd.py +++ b/pylops/signalprocessing/fftnd.py @@ -5,6 +5,7 @@ import numpy as np import numpy.typing as npt +import scipy.fft from pylops import LinearOperator from pylops.signalprocessing._baseffts import _BaseFFTND, _FFTNorms @@ -16,7 +17,8 @@ mkl_fft_message = deps.mkl_fft_import("the mkl fft module") if mkl_fft_message is None: - import mkl_fft.interfaces.numpy_fft as mkl_backend + import mkl_fft.interfaces.scipy_fft as mkl_backend + from mkl_fft.interfaces import _float_utils class _FFTND_numpy(_BaseFFTND): @@ -260,8 +262,10 @@ def __init__( @reshaped def _matvec(self, x: NDArray) -> NDArray: + x = _float_utils._downcast_float128_array(x) + x = _float_utils._upcast_float16_array(x) if self.ifftshift_before.any(): - x = mkl_backend.ifftshift(x, axes=self.axes[self.ifftshift_before]) + x = scipy.fft.ifftshift(x, axes=self.axes[self.ifftshift_before]) if not self.clinear: x = np.real(x) if self.real: @@ -279,13 +283,15 @@ def _matvec(self, x: NDArray) -> NDArray: if self.norm is _FFTNorms.ONE_OVER_N: y *= self._scale if self.fftshift_after.any(): - y = mkl_backend.fftshift(y, axes=self.axes[self.fftshift_after]) + y = scipy.fft.fftshift(y, axes=self.axes[self.fftshift_after]) return y @reshaped def _rmatvec(self, x: NDArray) -> NDArray: + x = _float_utils._downcast_float128_array(x) + x = _float_utils._upcast_float16_array(x) if self.fftshift_after.any(): - x = mkl_backend.ifftshift(x, axes=self.axes[self.fftshift_after]) + x = scipy.fft.ifftshift(x, axes=self.axes[self.fftshift_after]) if self.real: # Apply scaling to obtain a correct adjoint for this operator x = x.copy() @@ -309,7 +315,7 @@ def _rmatvec(self, x: NDArray) -> NDArray: if not self.clinear: y = np.real(y) if self.ifftshift_before.any(): - y = mkl_backend.fftshift(y, axes=self.axes[self.ifftshift_before]) + y = scipy.fft.fftshift(y, axes=self.axes[self.ifftshift_before]) return y def __truediv__(self, y: npt.ArrayLike) -> npt.ArrayLike: diff --git a/pytests/test_ffts.py b/pytests/test_ffts.py index 5fdc6066..c98e85a4 100644 --- a/pytests/test_ffts.py +++ b/pytests/test_ffts.py @@ -1,6 +1,5 @@ import itertools import os -import sys if int(os.environ.get("TEST_CUPY_PYLOPS", 0)): import cupy as np @@ -17,7 +16,7 @@ from pylops.optimization.basic import lsqr from pylops.signalprocessing import FFT, FFT2D, FFTND -from pylops.utils import dottest +from pylops.utils import dottest, mkl_fft_enabled # Utility function @@ -323,8 +322,8 @@ def test_unknown_engine(par): @pytest.mark.parametrize("par", pars_fft_small_real) def test_FFT_small_real(par): - if par["engine"] == "mkl_fft" and sys.platform == "darwin": - pytest.skip("mkl_fft not supported on macOS") + if par["engine"] == "mkl_fft" and not mkl_fft_enabled: + pytest.skip("mkl_fft is not installed") np.random.seed(5) if backend == "numpy" or (backend == "cupy" and par["engine"] == "numpy"): @@ -402,8 +401,8 @@ def test_FFT_small_real(par): ) @pytest.mark.parametrize("par", pars_fft_random_real) def test_FFT_random_real(par): - if par["engine"] == "mkl_fft" and sys.platform == "darwin": - pytest.skip("mkl_fft not supported on macOS") + if par["engine"] == "mkl_fft" and not mkl_fft_enabled: + pytest.skip("mkl_fft is not installed") np.random.seed(5) shape = par["shape"] @@ -460,8 +459,8 @@ def test_FFT_random_real(par): @pytest.mark.parametrize("par", pars_fft_small_cpx) def test_FFT_small_complex(par): - if par["engine"] == "mkl_fft" and sys.platform == "darwin": - pytest.skip("mkl_fft not supported on macOS") + if par["engine"] == "mkl_fft" and not mkl_fft_enabled: + pytest.skip("mkl_fft is not installed") np.random.seed(5) dtype, decimal = par["dtype_precision"] norm = par["norm"] @@ -534,8 +533,8 @@ def test_FFT_small_complex(par): @pytest.mark.parametrize("par", pars_fft_random_cpx) def test_FFT_random_complex(par): - if par["engine"] == "mkl_fft" and sys.platform == "darwin": - pytest.skip("mkl_fft not supported on macOS") + if par["engine"] == "mkl_fft" and not mkl_fft_enabled: + pytest.skip("mkl_fft is not installed") np.random.seed(5) if backend == "numpy" or (backend == "cupy" and par["engine"] == "numpy"): shape = par["shape"] @@ -623,8 +622,8 @@ def test_FFT_random_complex(par): ) @pytest.mark.parametrize("par", pars_fft2d_random_real) def test_FFT2D_random_real(par): - if par["engine"] == "mkl_fft" and sys.platform == "darwin": - pytest.skip("mkl_fft not supported on macOS") + if par["engine"] == "mkl_fft" and not mkl_fft_enabled: + pytest.skip("mkl_fft is not installed") np.random.seed(5) if backend == "numpy" or (backend == "cupy" and par["engine"] == "numpy"): shape = par["shape"] @@ -686,8 +685,8 @@ def test_FFT2D_random_real(par): @pytest.mark.parametrize("par", pars_fft2d_random_cpx) def test_FFT2D_random_complex(par): - if par["engine"] == "mkl_fft" and sys.platform == "darwin": - pytest.skip("mkl_fft not supported on macOS") + if par["engine"] == "mkl_fft" and not mkl_fft_enabled: + pytest.skip("mkl_fft is not installed") np.random.seed(5) if backend == "numpy" or (backend == "cupy" and par["engine"] == "numpy"): shape = par["shape"] @@ -769,8 +768,8 @@ def test_FFT2D_random_complex(par): @pytest.mark.parametrize("par", pars_fftnd_random_real) def test_FFTND_random_real(par): - if par["engine"] == "mkl_fft" and sys.platform == "darwin": - pytest.skip("mkl_fft not supported on macOS") + if par["engine"] == "mkl_fft" and not mkl_fft_enabled: + pytest.skip("mkl_fft is not installed") np.random.seed(5) if backend == "numpy" or (backend == "cupy" and par["engine"] == "numpy"): shape = par["shape"] @@ -832,8 +831,8 @@ def test_FFTND_random_real(par): @pytest.mark.parametrize("par", pars_fftnd_random_cpx) def test_FFTND_random_complex(par): - if par["engine"] == "mkl_fft" and sys.platform == "darwin": - pytest.skip("mkl_fft not supported on macOS") + if par["engine"] == "mkl_fft" and not mkl_fft_enabled: + pytest.skip("mkl_fft is not installed") np.random.seed(5) shape = par["shape"] dtype, decimal = par["dtype_precision"] @@ -910,8 +909,8 @@ def test_FFTND_random_complex(par): @pytest.mark.parametrize("par", pars_fft2dnd_small_cpx) def test_FFT2D_small_complex(par): - if par["engine"] == "mkl_fft" and sys.platform == "darwin": - pytest.skip("mkl_fft not supported on macOS") + if par["engine"] == "mkl_fft" and not mkl_fft_enabled: + pytest.skip("mkl_fft is not installed") np.random.seed(5) dtype, decimal = par["dtype_precision"] norm = par["norm"] @@ -962,8 +961,8 @@ def test_FFT2D_small_complex(par): @pytest.mark.parametrize("par", pars_fft2dnd_small_cpx) def test_FFTND_small_complex(par): - if par["engine"] == "mkl_fft" and sys.platform == "darwin": - pytest.skip("mkl_fft not supported on macOS") + if par["engine"] == "mkl_fft" and not mkl_fft_enabled: + pytest.skip("mkl_fft is not installed") np.random.seed(5) dtype, decimal = par["dtype_precision"] norm = par["norm"] @@ -1039,8 +1038,8 @@ def test_FFTND_small_complex(par): ], ) def test_FFT_1dsignal(par): - if par["engine"] == "mkl_fft" and sys.platform == "darwin": - pytest.skip("mkl_fft not supported on macOS") + if par["engine"] == "mkl_fft" and not mkl_fft_enabled: + pytest.skip("mkl_fft is not installed") np.random.seed(5) """Dot-test and inversion for FFT operator for 1d signal""" decimal = 3 if np.real(np.ones(1, par["dtype"])).dtype == np.float32 else 8 @@ -1158,8 +1157,8 @@ def test_FFT_2dsignal(par): """Dot-test and inversion for fft operator for 2d signal (fft on single dimension) """ - if par["engine"] == "mkl_fft" and sys.platform == "darwin": - pytest.skip("mkl_fft not supported on macOS") + if par["engine"] == "mkl_fft" and not mkl_fft_enabled: + pytest.skip("mkl_fft is not installed") np.random.seed(5) decimal = 3 if np.real(np.ones(1, par["dtype"])).dtype == np.float32 else 8 @@ -1375,8 +1374,8 @@ def test_FFT_3dsignal(par): """Dot-test and inversion for fft operator for 3d signal (fft on single dimension) """ - if par["engine"] == "mkl_fft" and sys.platform == "darwin": - pytest.skip("mkl_fft not supported on macOS") + if par["engine"] == "mkl_fft" and not mkl_fft_enabled: + pytest.skip("mkl_fft is not installed") np.random.seed(5) decimal = 3 if np.real(np.ones(1, par["dtype"])).dtype == np.float32 else 8 @@ -1603,8 +1602,8 @@ def test_FFT_3dsignal(par): ) def test_FFT2D(par): """Dot-test and inversion for FFT2D operator for 2d signal""" - if par["engine"] == "mkl_fft" and sys.platform == "darwin": - pytest.skip("mkl_fft not supported on macOS") + if par["engine"] == "mkl_fft" and not mkl_fft_enabled: + pytest.skip("mkl_fft is not installed") np.random.seed(5) decimal = 3 if np.real(np.ones(1, par["dtype"])).dtype == np.float32 else 8 @@ -1740,8 +1739,8 @@ def test_FFT2D(par): ) def test_FFT3D(par): """Dot-test and inversion for FFTND operator for 3d signal""" - if par["engine"] == "mkl_fft" and sys.platform == "darwin": - pytest.skip("mkl_fft not supported on macOS") + if par["engine"] == "mkl_fft" and not mkl_fft_enabled: + pytest.skip("mkl_fft is not installed") np.random.seed(5) decimal = 3 if np.real(np.ones(1, par["dtype"])).dtype == np.float32 else 8 diff --git a/requirements-dev.txt b/requirements-dev.txt index 9049bf4a..717584cd 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -1,6 +1,3 @@ ---index-url https://software.repos.intel.com/python/pypi ---extra-index-url https://pypi.org/simple -mkl_fft; sys_platform != "darwin" numpy>=2.0.0 scipy>=1.13.0 jax diff --git a/requirements-intel-mkl.txt b/requirements-intel-mkl.txt new file mode 100644 index 00000000..d407245e --- /dev/null +++ b/requirements-intel-mkl.txt @@ -0,0 +1,4 @@ +--index-url https://software.repos.intel.com/python/pypi +numpy>=2.0.0 +scipy>=1.13.0 +mkl_fft; sys_platform != "darwin" From da56788f50027e90434c7ef786cb373fc0f97f86 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Tue, 9 Sep 2025 00:37:15 +0530 Subject: [PATCH 09/19] Add warning for pyfftw --- docs/source/installation.rst | 24 ++++++++++++++++++++++-- requirements-doc.txt | 4 +--- 2 files changed, 23 insertions(+), 5 deletions(-) diff --git a/docs/source/installation.rst b/docs/source/installation.rst index 1776bc5d..a42d7fbb 100755 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -230,8 +230,7 @@ These are an option for an environment without ``conda`` that needs Intel MKL wi .. warning:: ``intel-numpy`` and ``intel-scipy`` not only link against Intel MKL, but also substitute NumPy and - SciPy FFTs for `Intel MKL FFT `_. **MKL FFT is not supported - and may break PyLops**. + SciPy FFTs for `Intel MKL FFT `_. Multithreading @@ -400,6 +399,27 @@ or via pip: .. note:: `mkl_fft` is not supported on macOS +While the library will work without NumPy with MKL and Intel Python, for maximum performance it is recommended to use +NumPy built with Intel’s Math Kernel Library (MKL) together with Intel Python. + +.. warning:: + + ``pyFFTW`` may not work correctly with NumPy + MKL. To avoid issues, it is recommended to build ``pyFFTW`` from + source after setting the ``STATIC_FFTW_DIR`` environment variable to the absolute path of the static FFTW + libraries. + + If the following environment variables are set before installing ``pyFFTW``, compatibility problems with MKL + should not occur: + + 1. ``export STATIC_FFTW_DIR=${PREFIX}/lib`` + (where ``${PREFIX}`` is the base of the current Anaconda environment with + the ``fftw`` package installed) + + 2. ``export CFLAGS="$CFLAGS -Wl,-Bsymbolic"`` + + Alternatively, you can install ``pyFFTW`` directly with ``conda``, since the updated recipe is already available + and works without any manual adjustments. + Numba ----- Although we always strive to write code for forward and adjoint operators that takes advantage of diff --git a/requirements-doc.txt b/requirements-doc.txt index cf520658..1d564081 100644 --- a/requirements-doc.txt +++ b/requirements-doc.txt @@ -1,6 +1,3 @@ ---index-url https://software.repos.intel.com/python/pypi ---extra-index-url https://pypi.org/simple -mkl_fft; sys_platform != "darwin" numpy>=2.0.0 scipy>=1.13.0 jax @@ -34,3 +31,4 @@ flake8 mypy pytensor>=2.28.0 pymc>=5.21.0 +mkl_fft; sys_platform != "darwin" From a5645947e211051ae9d573fadc542ddb68cd723b Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sat, 27 Sep 2025 22:05:25 +0100 Subject: [PATCH 10/19] minor: added mkl_fft warning and small docstring improvements --- docs/source/installation.rst | 7 +++---- examples/plot_fft.py | 2 +- pylops/signalprocessing/fft.py | 6 ++++-- pylops/signalprocessing/fftnd.py | 4 ++-- 4 files changed, 10 insertions(+), 9 deletions(-) diff --git a/docs/source/installation.rst b/docs/source/installation.rst index a42d7fbb..eeb313d6 100755 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -359,6 +359,9 @@ FFTW and MKL-FFT ---------------- Four different "engines" are provided by the :py:class:`pylops.signalprocessing.FFT` operator: ``engine="numpy"`` (default), ``engine="scipy"``, ``engine="fftw"`` and ``engine="mkl_fft"``. +Similarly, the :py:class:`pylops.signalprocessing.FFT2D` and +the :py:class:`pylops.signalprocessing.FFTND` operators come with three "engines", namely +``engine="numpy"`` (default), ``engine="scipy"``, and ``engine="mkl_fft"``. The first two engines are part of the required PyLops dependencies. The third implements the well-known `FFTW `_ @@ -376,10 +379,6 @@ or via pip: >> pip install pyfftw -.. note:: - FFTW is only available for :py:class:`pylops.signalprocessing.FFT`, - not :py:class:`pylops.signalprocessing.FFT2D` or :py:class:`pylops.signalprocessing.FFTND`. - The fourth implements Intel MKL FFT via the Python interface `mkl_fft `_. This provides access to Intel’s oneMKL Fourier Transform routines, enabling efficient FFT computations with performance close to native C/Intel® oneMKL diff --git a/examples/plot_fft.py b/examples/plot_fft.py index 5718e252..cd50ac59 100644 --- a/examples/plot_fft.py +++ b/examples/plot_fft.py @@ -67,7 +67,7 @@ plt.tight_layout() ############################################################################### -# PyLops also has a third FFT engine (engine='mkl_fft') that uses the well-known +# PyLops also has a third FFT engine (``engine='mkl_fft'``) that uses the well-known # `Intel MKL FFT `_. This is a Python wrapper around # the `Intel® oneAPI Math Kernel Library (oneMKL) `_ # Fourier Transform functions. It lets PyLops run discrete Fourier transforms faster diff --git a/pylops/signalprocessing/fft.py b/pylops/signalprocessing/fft.py index c5d4a15b..2b87838e 100644 --- a/pylops/signalprocessing/fft.py +++ b/pylops/signalprocessing/fft.py @@ -16,7 +16,7 @@ from pylops.utils.typing import DTypeLike, InputDimsLike, NDArray pyfftw_message = deps.pyfftw_import("the fft module") -mkl_fft_message = deps.mkl_fft_import("the mkl fft module") +mkl_fft_message = deps.mkl_fft_import("the fft module") if pyfftw_message is None: import pyfftw @@ -696,6 +696,8 @@ def FFT( ): if engine == "fftw" and pyfftw_message is not None: logger.warning(pyfftw_message) + if engine == "mkl_fft" and mkl_fft_message is not None: + logger.warning(mkl_fft_message) f = _FFT_numpy( dims, axis=axis, @@ -722,6 +724,6 @@ def FFT( **kwargs_fft, ) else: - raise NotImplementedError("engine must be numpy, fftw, scipy or mkl_fft") + raise NotImplementedError("engine must be numpy, scipy, fftw, or mkl_fft") f.name = name return f diff --git a/pylops/signalprocessing/fftnd.py b/pylops/signalprocessing/fftnd.py index 2754cf8d..cb46a9f3 100644 --- a/pylops/signalprocessing/fftnd.py +++ b/pylops/signalprocessing/fftnd.py @@ -425,7 +425,7 @@ def FFTND( engine : :obj:`str`, optional .. versionadded:: 1.17.0 - Engine used for fft computation (``numpy`` or ``scipy``). Choose + Engine used for fft computation (``numpy`` or ``scipy`` or ``mkl_fft``). Choose ``numpy`` when working with cupy and jax arrays. dtype : :obj:`str`, optional Type of elements in input array. Note that the ``dtype`` of the operator @@ -483,7 +483,7 @@ def FFTND( the same dimension ``axes``. - If ``norm`` is not one of "ortho", "none", or "1/n". NotImplementedError - If ``engine`` is neither ``numpy``, nor ``scipy``. + If ``engine`` is neither ``numpy``, ``scipy`` nor ``mkl_fft``. Notes ----- From c463da2de88247858a7cf9a6c9f6a540c1c87d00 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sat, 27 Sep 2025 22:10:48 +0100 Subject: [PATCH 11/19] minor: improved mkl_fft_message --- pylops/utils/deps.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/pylops/utils/deps.py b/pylops/utils/deps.py index 4232cc5c..dad41381 100644 --- a/pylops/utils/deps.py +++ b/pylops/utils/deps.py @@ -254,8 +254,11 @@ def mkl_fft_import(message): "mkl_fft not available, reverting to numpy. " "In order to be able to use " f"{message} run " - f'"pip install mkl_fft" or ' - f'"conda install -c conda-forge mkl_fft".' + '"pip install --index-url ' + "https://software.repos.intel.com/python/pypi " + '--extra-index-url https://pypi.org/simple mkl_fft" ' + 'or "conda install -c https://software.repos.intel.com/python/conda ' + '-c conda-forge mkl_fft".' ) return mkl_fft_message From 734ed762f300cc05a570fb6c5c74a2786212eac6 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sat, 27 Sep 2025 22:18:12 +0100 Subject: [PATCH 12/19] minor: rearrange steps in _FFT2D_mklfft init to align with _FFT2D_numpy --- pylops/signalprocessing/fft.py | 4 ++-- pylops/signalprocessing/fft2d.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/pylops/signalprocessing/fft.py b/pylops/signalprocessing/fft.py index 2b87838e..b459519f 100644 --- a/pylops/signalprocessing/fft.py +++ b/pylops/signalprocessing/fft.py @@ -340,7 +340,7 @@ def _matvec(self, x: NDArray) -> NDArray: elif self.doifftpad: x = np.take(x, range(0, self.nfft), axis=self.axis) - # self.fftplan() always uses byte-alligned self.x as input array and + # self.fftplan() always uses byte-aligned self.x as input array and # returns self.y as output array. As such, self.x must be copied so as # not to be overwritten on a subsequent call to _matvec. np.copyto(self.x, x) @@ -362,7 +362,7 @@ def _rmatvec(self, x: NDArray) -> NDArray: if self.fftshift_after: x = np.fft.ifftshift(x, axes=self.axis) - # self.ifftplan() always uses byte-alligned self.y as input array. + # self.ifftplan() always uses byte-aligned self.y as input array. # We copy here so we don't need to copy again in the case of `real=True`, # which only performs operations that preserve byte-allignment. np.copyto(self.y, x) diff --git a/pylops/signalprocessing/fft2d.py b/pylops/signalprocessing/fft2d.py index 3f79364f..34d01d9a 100644 --- a/pylops/signalprocessing/fft2d.py +++ b/pylops/signalprocessing/fft2d.py @@ -276,10 +276,10 @@ def __init__( if self.naxes != 2: raise ValueError("FFT2D must be applied along exactly two dimensions") - self._kwargs_fft = kwargs_fft self.f1, self.f2 = self.fs del self.fs + self._kwargs_fft = kwargs_fft self._norm_kwargs: Dict[str, Union[None, str]] = {"norm": None} if self.norm is _FFTNorms.ORTHO: self._norm_kwargs["norm"] = "ortho" From 253ec9a72216f58778e4e21e654239a60e2279f1 Mon Sep 17 00:00:00 2001 From: Carlos da Costa Date: Sun, 28 Sep 2025 13:35:07 -0700 Subject: [PATCH 13/19] fix: doc improvements --- examples/plot_fft.py | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/examples/plot_fft.py b/examples/plot_fft.py index cd50ac59..a0b3f3e0 100644 --- a/examples/plot_fft.py +++ b/examples/plot_fft.py @@ -6,6 +6,7 @@ and :py:class:`pylops.signalprocessing.FFTND` operators to apply the Fourier Transform to the model and the inverse Fourier Transform to the data. """ + import matplotlib.pyplot as plt import numpy as np @@ -15,7 +16,7 @@ ############################################################################### # Let's start by applying the one dimensional FFT to a one dimensional -# sinusoidal signal :math:`d(t)=sin(2 \pi f_0t)` using a time axis of +# sinusoidal signal :math:`d(t)=\sin(2 \pi f_0t)` using a time axis of # lenght :math:`nt` and sampling :math:`dt` dt = 0.005 nt = 100 @@ -91,7 +92,7 @@ plt.tight_layout() ############################################################################### -# We can also apply the one dimensional FFT to to a two-dimensional +# We can also apply the one dimensional FFT to a two-dimensional # signal (along one of the first axis) dt = 0.005 nt, nx = 100, 20 @@ -124,7 +125,7 @@ fig.tight_layout() ############################################################################### -# We can also apply the two dimensional FFT to to a two-dimensional signal +# We can also apply the two dimensional FFT to a two-dimensional signal dt, dx = 0.005, 5 nt, nx = 100, 201 t = np.arange(nt) * dt @@ -161,7 +162,7 @@ ############################################################################### -# Finally can apply the three dimensional FFT to to a three-dimensional signal +# Finally can apply the three dimensional FFT to a three-dimensional signal dt, dx, dy = 0.005, 5, 3 nt, nx, ny = 30, 21, 11 t = np.arange(nt) * dt @@ -200,3 +201,15 @@ axs[1][1].set_title("Error") axs[1][1].axis("tight") fig.tight_layout() + +############################################################################### +# Supported Backends +# ~~~~~~~~~~~~~~~~~~ +# ======= ==================== ========= +# Backend Supported Dimensions Dependecy +# ======= ==================== ========= +# Numpy 1D, 2D, ND ``numpy`` (included) +# Scipy 1D, 2D, ND ``scipy`` (included) +# FFTW 1D ``pyfftw`` +# MKL 1D, 2D, ND ``mkl_fft``, or ``intel-numpy``/``intel-scipy`` via standard "numpy"/"scipy" engines +# ======= ==================== ========= From fb56ac03bd10d897066bbab4d5f17d148f1a8ec4 Mon Sep 17 00:00:00 2001 From: Carlos da Costa Date: Sun, 28 Sep 2025 19:39:47 -0700 Subject: [PATCH 14/19] fix: annotations --- pylops/utils/deps.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/pylops/utils/deps.py b/pylops/utils/deps.py index dad41381..67a55464 100644 --- a/pylops/utils/deps.py +++ b/pylops/utils/deps.py @@ -20,7 +20,7 @@ # error message at import of available package -def cupy_import(message: Optional[str] = None) -> str: +def cupy_import(message: Optional[str] = None) -> str | None: # detect if cupy is available and the user is expecting to be used cupy_test = ( util.find_spec("cupy") is not None and int(os.getenv("CUPY_PYLOPS", 1)) == 1 @@ -54,7 +54,7 @@ def cupy_import(message: Optional[str] = None) -> str: return cupy_message -def jax_import(message: Optional[str] = None) -> str: +def jax_import(message: Optional[str] = None) -> str | None: jax_test = ( util.find_spec("jax") is not None and int(os.getenv("JAX_PYLOPS", 1)) == 1 ) @@ -82,7 +82,7 @@ def jax_import(message: Optional[str] = None) -> str: return jax_message -def devito_import(message: Optional[str] = None) -> str: +def devito_import(message: Optional[str] = None) -> str | None: if devito_enabled: try: import_module("devito") # noqa: F401 @@ -99,7 +99,7 @@ def devito_import(message: Optional[str] = None) -> str: return devito_message -def dtcwt_import(message: Optional[str] = None) -> str: +def dtcwt_import(message: Optional[str] = None) -> str | None: if dtcwt_enabled: try: import dtcwt # noqa: F401 @@ -116,7 +116,7 @@ def dtcwt_import(message: Optional[str] = None) -> str: return dtcwt_message -def numba_import(message: Optional[str] = None) -> str: +def numba_import(message: Optional[str] = None) -> str | None: if numba_enabled: try: import_module("numba") # noqa: F401 @@ -135,7 +135,7 @@ def numba_import(message: Optional[str] = None) -> str: return numba_message -def pyfftw_import(message: Optional[str] = None) -> str: +def pyfftw_import(message: Optional[str] = None) -> str | None: if pyfftw_enabled: try: import_module("pyfftw") # noqa: F401 @@ -154,7 +154,7 @@ def pyfftw_import(message: Optional[str] = None) -> str: return pyfftw_message -def pywt_import(message: Optional[str] = None) -> str: +def pywt_import(message: Optional[str] = None) -> str | None: if pywt_enabled: try: import_module("pywt") # noqa: F401 @@ -173,7 +173,7 @@ def pywt_import(message: Optional[str] = None) -> str: return pywt_message -def skfmm_import(message: Optional[str] = None) -> str: +def skfmm_import(message: Optional[str] = None) -> str | None: if skfmm_enabled: try: import_module("skfmm") # noqa: F401 @@ -191,7 +191,7 @@ def skfmm_import(message: Optional[str] = None) -> str: return skfmm_message -def spgl1_import(message: Optional[str] = None) -> str: +def spgl1_import(message: Optional[str] = None) -> str | None: if spgl1_enabled: try: import_module("spgl1") # noqa: F401 @@ -208,7 +208,7 @@ def spgl1_import(message: Optional[str] = None) -> str: return spgl1_message -def sympy_import(message: Optional[str] = None) -> str: +def sympy_import(message: Optional[str] = None) -> str | None: if sympy_enabled: try: import_module("sympy") # noqa: F401 @@ -225,7 +225,7 @@ def sympy_import(message: Optional[str] = None) -> str: return sympy_message -def pytensor_import(message: Optional[str] = None) -> str: +def pytensor_import(message: Optional[str] = None) -> str | None: if pytensor_enabled: try: import_module("pytensor") # noqa: F401 @@ -242,7 +242,7 @@ def pytensor_import(message: Optional[str] = None) -> str: return pytensor_message -def mkl_fft_import(message): +def mkl_fft_import(message: Optional[str]) -> str | None: if mkl_fft_enabled: try: import_module("mkl_fft") # noqa: F401 From b4db490514450bb48ea1f332789874ca6cfc4c0f Mon Sep 17 00:00:00 2001 From: Carlos da Costa Date: Sun, 28 Sep 2025 20:29:23 -0700 Subject: [PATCH 15/19] fix: add cupy --- examples/plot_fft.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/examples/plot_fft.py b/examples/plot_fft.py index a0b3f3e0..47cf98a1 100644 --- a/examples/plot_fft.py +++ b/examples/plot_fft.py @@ -205,11 +205,11 @@ ############################################################################### # Supported Backends # ~~~~~~~~~~~~~~~~~~ -# ======= ==================== ========= -# Backend Supported Dimensions Dependecy -# ======= ==================== ========= -# Numpy 1D, 2D, ND ``numpy`` (included) -# Scipy 1D, 2D, ND ``scipy`` (included) -# FFTW 1D ``pyfftw`` -# MKL 1D, 2D, ND ``mkl_fft``, or ``intel-numpy``/``intel-scipy`` via standard "numpy"/"scipy" engines -# ======= ==================== ========= +# ========== ==================== ========= +# Backend Supported Dimensions Dependecy +# ========== ==================== ========= +# Numpy/CuPy 1D, 2D, ND ``numpy`` (included) +# Scipy 1D, 2D, ND ``scipy`` (included) +# FFTW 1D ``pyfftw`` +# MKL 1D, 2D, ND ``mkl_fft``, or ``intel-numpy``/``intel-scipy`` via standard "numpy"/"scipy" engines +# ========== ==================== ========= From 0b8e92d608d307ecb422da9558d3899b1930d1de Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Fri, 3 Oct 2025 15:35:09 +0530 Subject: [PATCH 16/19] Update GA and add content for mkl-fft --- .github/workflows/build-mkl.yaml | 9 +-------- docs/source/installation.rst | 10 +++++++--- 2 files changed, 8 insertions(+), 11 deletions(-) diff --git a/.github/workflows/build-mkl.yaml b/.github/workflows/build-mkl.yaml index c086d826..4e11dcd7 100644 --- a/.github/workflows/build-mkl.yaml +++ b/.github/workflows/build-mkl.yaml @@ -22,16 +22,9 @@ jobs: - uses: conda-incubator/setup-miniconda@v3.2.0 with: use-mamba: true - channels: https://software.repos.intel.com/python/conda, conda-forge - conda-remove-defaults: true python-version: ${{ matrix.python-version }} activate-environment: mkl-test-env - - name: Install dependencies - run: | - conda install -y pyfftw - pip install -r requirements-intel-mkl.txt - pip install -r requirements-dev.txt - pip install -r requirements-torch.txt + environment-file: environment-dev-intel-mkl.yml - name: Install pylops run: | python -m setuptools_scm diff --git a/docs/source/installation.rst b/docs/source/installation.rst index eeb313d6..f7b80aea 100755 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -379,7 +379,7 @@ or via pip: >> pip install pyfftw -The fourth implements Intel MKL FFT via the Python interface `mkl_fft `_. +The fourth implements **Intel MKL FFT** via the Python interface `mkl_fft `_. This provides access to Intel’s oneMKL Fourier Transform routines, enabling efficient FFT computations with performance close to native C/Intel® oneMKL @@ -398,8 +398,12 @@ or via pip: .. note:: `mkl_fft` is not supported on macOS -While the library will work without NumPy with MKL and Intel Python, for maximum performance it is recommended to use -NumPy built with Intel’s Math Kernel Library (MKL) together with Intel Python. +Installing ``mkl-fft`` triggers the installation of Intel-optimized versions of `NumPy `__ and +`SciPy `__, which redirects ``numpy.fft`` and ``scipy.fft`` to use MKL FFT routines. +As a result, all FFT operations and computational backends leverage Intel MKL for optimal performance. + +Although the library can run without Intel-optimized NumPy and SciPy, maximum performance is achieved when using NumPy and +SciPy built with Intel’s Math Kernel Library (MKL) alongside Intel Python. .. warning:: From 1d407fb8446efe6af86bd6f1a4b7cf2a13fd347a Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Fri, 3 Oct 2025 15:51:06 +0530 Subject: [PATCH 17/19] Minor change in py version --- environment-dev-intel-mkl.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/environment-dev-intel-mkl.yml b/environment-dev-intel-mkl.yml index e5bf97e2..7ce56bd6 100644 --- a/environment-dev-intel-mkl.yml +++ b/environment-dev-intel-mkl.yml @@ -5,7 +5,7 @@ channels: - defaults - numba dependencies: - - python==3.12.8 + - python>=3.10.0 - pip - numpy>=2.0.0 - scipy>=1.13.0 From bdb62bd40fcff10ba43af7d6ce988e66499a9bb3 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Sat, 4 Oct 2025 00:31:42 +0530 Subject: [PATCH 18/19] Add pyfftw in requirements-pyfftw.txt --- .github/workflows/build-mkl.yaml | 9 ++++++++- .github/workflows/build.yaml | 1 + .github/workflows/codacy-coverage-reporter.yaml | 1 + Makefile | 1 + azure-pipelines.yml | 2 ++ requirements-dev.txt | 1 - requirements-pyfftw.txt | 1 + 7 files changed, 14 insertions(+), 2 deletions(-) create mode 100644 requirements-pyfftw.txt diff --git a/.github/workflows/build-mkl.yaml b/.github/workflows/build-mkl.yaml index 4e11dcd7..c086d826 100644 --- a/.github/workflows/build-mkl.yaml +++ b/.github/workflows/build-mkl.yaml @@ -22,9 +22,16 @@ jobs: - uses: conda-incubator/setup-miniconda@v3.2.0 with: use-mamba: true + channels: https://software.repos.intel.com/python/conda, conda-forge + conda-remove-defaults: true python-version: ${{ matrix.python-version }} activate-environment: mkl-test-env - environment-file: environment-dev-intel-mkl.yml + - name: Install dependencies + run: | + conda install -y pyfftw + pip install -r requirements-intel-mkl.txt + pip install -r requirements-dev.txt + pip install -r requirements-torch.txt - name: Install pylops run: | python -m setuptools_scm diff --git a/.github/workflows/build.yaml b/.github/workflows/build.yaml index ac5788ad..c430ae7b 100644 --- a/.github/workflows/build.yaml +++ b/.github/workflows/build.yaml @@ -25,6 +25,7 @@ jobs: python -m pip install --upgrade pip setuptools pip install flake8 pytest pip install -r requirements-dev.txt + pip install -r requirements-pyfftw.txt pip install -r requirements-torch.txt - name: Install pylops run: | diff --git a/.github/workflows/codacy-coverage-reporter.yaml b/.github/workflows/codacy-coverage-reporter.yaml index 113357ad..239af9f0 100644 --- a/.github/workflows/codacy-coverage-reporter.yaml +++ b/.github/workflows/codacy-coverage-reporter.yaml @@ -27,6 +27,7 @@ jobs: python -m pip install --upgrade pip pip install flake8 pytest pip install -r requirements-dev.txt + pip install -r requirements-pyfftw.txt pip install -r requirements-torch.txt - name: Install pylops run: | diff --git a/Makefile b/Makefile index 3df51181..64940b37 100644 --- a/Makefile +++ b/Makefile @@ -22,6 +22,7 @@ install: dev-install: make pipcheck $(PIP) install -r requirements-dev.txt &&\ + $(PIP) install -r requirements-pyfftw.txt &&\ $(PIP) install -r requirements-torch.txt && $(PIP) install -e . dev-install_intel_mkl: diff --git a/azure-pipelines.yml b/azure-pipelines.yml index 049f19bf..132eedff 100644 --- a/azure-pipelines.yml +++ b/azure-pipelines.yml @@ -62,6 +62,7 @@ jobs: - script: | python -m pip install --upgrade pip setuptools wheel django pip install -r requirements-dev.txt + pip install -r requirements-pyfftw.txt pip install -r requirements-torch.txt pip install . displayName: 'Install prerequisites and library' @@ -92,6 +93,7 @@ jobs: - script: | python -m pip install --upgrade pip setuptools wheel django pip install -r requirements-dev.txt + pip install -r requirements-pyfftw.txt pip install -r requirements-torch.txt pip install . displayName: 'Install prerequisites and library' diff --git a/requirements-dev.txt b/requirements-dev.txt index 717584cd..6d18813e 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -2,7 +2,6 @@ numpy>=2.0.0 scipy>=1.13.0 jax numba -pyfftw PyWavelets spgl1 scikit-fmm diff --git a/requirements-pyfftw.txt b/requirements-pyfftw.txt new file mode 100644 index 00000000..0ee0a829 --- /dev/null +++ b/requirements-pyfftw.txt @@ -0,0 +1 @@ +pyfftw From 059ab717bc5c1cc7a3309a8468f06918f514bc08 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sat, 4 Oct 2025 23:04:01 +0100 Subject: [PATCH 19/19] minor: small updates to doc --- docs/source/installation.rst | 53 +++++++++++++++++++----------------- examples/plot_fft.py | 12 ++++++-- 2 files changed, 37 insertions(+), 28 deletions(-) diff --git a/docs/source/installation.rst b/docs/source/installation.rst index f7b80aea..f4ff87fa 100755 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -209,7 +209,7 @@ library will ensure optimal performance of PyLops when using only *required depe We strongly encourage using the Anaconda Python distribution as NumPy and SciPy will, when available, be automatically linked to `Intel MKL `_, the most performant library for basic linear algebra -operations to date (see `Markus Beuckelmann's benchmarks `_). +operations to date. The PyPI version installed with ``pip``, however, will default to `OpenBLAS `_. For more information, see `NumPy's section on BLAS `_. @@ -224,13 +224,13 @@ run the following commands in a Python interpreter: print(sp.__config__.show()) -Intel also provides `NumPy `__ and `SciPy `__ replacement packages in PyPI ``intel-numpy`` and ``intel-scipy``, respectively, which link to Intel MKL. +Intel also provides `NumPy `__ and `SciPy `__ replacement packages in PyPI, namely ``intel-numpy`` and ``intel-scipy``, which link to Intel MKL. These are an option for an environment without ``conda`` that needs Intel MKL without requiring manual compilation. .. warning:: ``intel-numpy`` and ``intel-scipy`` not only link against Intel MKL, but also substitute NumPy and - SciPy FFTs for `Intel MKL FFT `_. + SciPy FFTs with `Intel MKL FFT `_. Multithreading @@ -296,7 +296,7 @@ of PyLops in such a way that if an *optional* dependency is not present in your a safe fallback to one of the required dependencies will be enforced. When available in your system, we recommend using the Conda package manager and install all the -required and optional dependencies of PyLops at once using the command: +required and some of the optional dependencies of PyLops at once using the command: .. code-block:: bash @@ -304,17 +304,19 @@ required and optional dependencies of PyLops at once using the command: in this case all dependencies will be installed from their Conda distributions. -Alternatively, from version ``1.4.0`` optional dependencies can also be installed as -part of the pip installation via: +Alternatively, from version ``1.4.0`` some of the optional dependencies can also be +installed as part of the pip installation via: .. code-block:: bash >> pip install pylops[advanced] Dependencies are however installed from their PyPI wheels. -An exception is however represented by CuPy. This library is **not** installed + +Finally, note that CuPy and JAX are not **not** installed automatically. Users interested to accelerate their computations with the aid -of GPUs should install it prior to installing PyLops as described in :ref:`OptionalGPU`. +of GPUs should install either or both of them prior to installing PyLops as +described in :ref:`OptionalGPU`. .. note:: @@ -323,11 +325,13 @@ of GPUs should install it prior to installing PyLops as described in :ref:`Optio PyLops via ``make dev-install_conda`` (``conda``) or ``make dev-install`` (``pip``). -In alphabetic order: +More details about the installation process for the different optional dependencies are described +in the following (an asterisc is used to indicate those dependencies that are automatically installed +when installing PyLops from conda-forge or via ``pip install pylops[advanced]``): -dtcwt ------ +dtcwt* +------ `dtcwt `_ is a library used to implement the DT-CWT operators. @@ -355,8 +359,8 @@ Install it via ``pip`` with >> pip install devito -FFTW and MKL-FFT ----------------- +FFTW* and MKL-FFT +----------------- Four different "engines" are provided by the :py:class:`pylops.signalprocessing.FFT` operator: ``engine="numpy"`` (default), ``engine="scipy"``, ``engine="fftw"`` and ``engine="mkl_fft"``. Similarly, the :py:class:`pylops.signalprocessing.FFT2D` and @@ -395,9 +399,6 @@ or via pip: >> pip install --index-url https://software.repos.intel.com/python/pypi --extra-index-url https://pypi.org/simple mkl_fft -.. note:: - `mkl_fft` is not supported on macOS - Installing ``mkl-fft`` triggers the installation of Intel-optimized versions of `NumPy `__ and `SciPy `__, which redirects ``numpy.fft`` and ``scipy.fft`` to use MKL FFT routines. As a result, all FFT operations and computational backends leverage Intel MKL for optimal performance. @@ -405,6 +406,9 @@ As a result, all FFT operations and computational backends leverage Intel MKL fo Although the library can run without Intel-optimized NumPy and SciPy, maximum performance is achieved when using NumPy and SciPy built with Intel’s Math Kernel Library (MKL) alongside Intel Python. +.. note:: + `mkl_fft` is not supported on macOS + .. warning:: ``pyFFTW`` may not work correctly with NumPy + MKL. To avoid issues, it is recommended to build ``pyFFTW`` from @@ -423,8 +427,8 @@ SciPy built with Intel’s Math Kernel Library (MKL) alongside Intel Python. Alternatively, you can install ``pyFFTW`` directly with ``conda``, since the updated recipe is already available and works without any manual adjustments. -Numba ------ +Numba* +------ Although we always strive to write code for forward and adjoint operators that takes advantage of the perks of NumPy and SciPy (e.g., broadcasting, ufunc), in some case we may end up using for loops that may lead to poor performance. In those cases we may decide to implement alternative (optional) @@ -483,7 +487,6 @@ It can also be checked dynamically with ``numba.config.NUMBA_DEFAULT_NUM_THREADS PyMC and PyTensor ----------------- - `PyTensor `_ is used to allow seamless integration between PyLops and `PyMC `_ operators. Install both of them via ``conda`` with: @@ -502,8 +505,8 @@ or via ``pip`` with OSX users may experience a ``CompileError`` error when using PyTensor. This can be solved by adding ``pytensor.config.gcc__cxxflags = "-Wno-c++11-narrowing"`` after ``import pytensor``. -PyWavelets ----------- +PyWavelets* +----------- `PyWavelets `_ is used to implement the wavelet operators. Install it via ``conda`` with: @@ -518,8 +521,8 @@ or via ``pip`` with >> pip install PyWavelets -scikit-fmm ----------- +scikit-fmm* +----------- `scikit-fmm `_ is a library which implements the fast marching method. It is used in PyLops to compute traveltime tables in the initialization of :py:class:`pylops.waveeqprocessing.Kirchhoff` @@ -537,8 +540,8 @@ or with ``pip`` via >> pip install scikit-fmm -SPGL1 ------ +SPGL1* +------ `SPGL1 `_ is used to solve sparsity-promoting basis pursuit, basis pursuit denoise, and Lasso problems in :py:func:`pylops.optimization.sparsity.SPGL1` solver. diff --git a/examples/plot_fft.py b/examples/plot_fft.py index 47cf98a1..a41e8ffe 100644 --- a/examples/plot_fft.py +++ b/examples/plot_fft.py @@ -203,11 +203,17 @@ fig.tight_layout() ############################################################################### +# To conclude, we provide a summary table of the different backends +# supported by :py:class:`pylops.signalprocessing.FFT`, +# :py:class:`pylops.signalprocessing.FFT2D` +# and :py:class:`pylops.signalprocessing.FFTND` operators and +# third-party dependencies are required to be able to use them. +# # Supported Backends # ~~~~~~~~~~~~~~~~~~ -# ========== ==================== ========= -# Backend Supported Dimensions Dependecy -# ========== ==================== ========= +# ========== ==================== ========== +# Backend Supported Dimensions Dependency +# ========== ==================== ========== # Numpy/CuPy 1D, 2D, ND ``numpy`` (included) # Scipy 1D, 2D, ND ``scipy`` (included) # FFTW 1D ``pyfftw``