From d30cc4e6c4e8de1b0282a7f87ca99082357dbd53 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Tue, 7 Feb 2023 20:13:18 +0530 Subject: [PATCH 01/48] Added mkl-fft along with examples and tests --- environment-dev.yml | 3 + examples/plot_fft.py | 130 +++++++++++++++++++++++++++ pylops/signalprocessing/fft.py | 112 ++++++++++++++++++++++- pylops/signalprocessing/fft2d.py | 120 ++++++++++++++++++++++++- pylops/signalprocessing/fftnd.py | 149 ++++++++++++++++++++++++++++++- pytests/test_ffts.py | 119 +++++++++++++++++++----- requirements-dev.txt | 3 + 7 files changed, 604 insertions(+), 32 deletions(-) diff --git a/environment-dev.yml b/environment-dev.yml index f00830c2..99687625 100755 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -7,6 +7,9 @@ channels: dependencies: - python>=3.6.4 - pip + - mkl + - mkl-fft + - mkl-service - numpy>=1.21.0 - scipy>=1.4.0 - pytorch>=1.2.0 diff --git a/examples/plot_fft.py b/examples/plot_fft.py index a8af9da6..6bb6b01a 100644 --- a/examples/plot_fft.py +++ b/examples/plot_fft.py @@ -66,6 +66,26 @@ axs[1].set_xlim([0, 3 * f0]) plt.tight_layout() +############################################################################### +# PyLops implements a third engine (``engine='mkl_fft'``) which uses the +# well-known `IntelPython mkl_fft ` +FFTop = pylops.signalprocessing.FFT(dims=nt, nfft=nfft, sampling=dt, engine="mkl_fft") +D = FFTop * d + +# Adjoint = 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) @@ -99,6 +119,39 @@ axs[1][1].axis("tight") fig.tight_layout() +############################################################################### +# We can also apply the one dimensional FFT to a two-dimensional +# signal (along one of the first axis) using intel mkl_fft +dt = 0.005 +nt, nx = 100, 20 +t = np.arange(nt) * dt +f0 = 10 +nfft = 2**10 +d = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1) + +FFTop = pylops.signalprocessing.FFT(dims=(nt, nx), axis=0, nfft=nfft, sampling=dt, engine="mkl_fft") +D = FFTop * d.ravel() + +# Adjoint = inverse for FFT +dinv = FFTop.H * D +dinv = FFTop / D +dinv = np.real(dinv).reshape(nt, nx) + +fig, axs = plt.subplots(2, 2, figsize=(10, 6)) +axs[0][0].imshow(d, vmin=-20, vmax=20, cmap="bwr") +axs[0][0].set_title("Signal") +axs[0][0].axis("tight") +axs[0][1].imshow(np.abs(D.reshape(nfft, nx)[:200, :]), cmap="bwr") +axs[0][1].set_title("Fourier Transform using mkl_fft") +axs[0][1].axis("tight") +axs[1][0].imshow(dinv, vmin=-20, vmax=20, cmap="bwr") +axs[1][0].set_title("Inverted") +axs[1][0].axis("tight") +axs[1][1].imshow(d - dinv, vmin=-20, vmax=20, cmap="bwr") +axs[1][1].set_title("Error") +axs[1][1].axis("tight") +fig.tight_layout() + ############################################################################### # We can also apply the two dimensional FFT to to a two-dimensional signal dt, dx = 0.005, 5 @@ -135,6 +188,42 @@ axs[1][1].axis("tight") fig.tight_layout() +############################################################################### +# We can also apply the two-dimensional FFT to a two-dimensional signal using intel mkl_fft +dt, dx = 0.005, 5 +nt, nx = 100, 201 +t = np.arange(nt) * dt +x = np.arange(nx) * dx +f0 = 10 +nfft = 2**10 +d = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1) + +FFTop = pylops.signalprocessing.FFT2D( + dims=(nt, nx), nffts=(nfft, nfft), sampling=(dt, dx), engine='mkl_fft' +) +D = FFTop * d.ravel() + +dinv = FFTop.H * D +dinv = FFTop / D +dinv = np.real(dinv).reshape(nt, nx) + +fig, axs = plt.subplots(2, 2, figsize=(10, 6)) +axs[0][0].imshow(d, vmin=-100, vmax=100, cmap="bwr") +axs[0][0].set_title("Signal") +axs[0][0].axis("tight") +axs[0][1].imshow( + np.abs(np.fft.fftshift(D.reshape(nfft, nfft), axes=1)[:200, :]), cmap="bwr" +) +axs[0][1].set_title("Fourier Transform 2D mkl_fft") +axs[0][1].axis("tight") +axs[1][0].imshow(dinv, vmin=-100, vmax=100, cmap="bwr") +axs[1][0].set_title("Inverted") +axs[1][0].axis("tight") +axs[1][1].imshow(d - dinv, vmin=-100, vmax=100, cmap="bwr") +axs[1][1].set_title("Error") +axs[1][1].axis("tight") +fig.tight_layout() + ############################################################################### # Finally can apply the three dimensional FFT to to a three-dimensional signal @@ -176,3 +265,44 @@ axs[1][1].set_title("Error") axs[1][1].axis("tight") fig.tight_layout() + +############################################################################### +# Finally can apply the three-dimensional FFT to a three-dimensional signal using mkl_fft +dt, dx, dy = 0.005, 5, 3 +nt, nx, ny = 30, 21, 11 +t = np.arange(nt) * dt +x = np.arange(nx) * dx +y = np.arange(nx) * dy +f0 = 10 +nfft = 2**6 +nfftk = 2**5 + +d = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1) +d = np.tile(d[:, :, np.newaxis], [1, 1, ny]) + +FFTop = pylops.signalprocessing.FFTND( + dims=(nt, nx, ny), nffts=(nfft, nfftk, nfftk), sampling=(dt, dx, dy), engine="mkl_fft" +) +D = FFTop * d.ravel() + +dinv = FFTop.H * D +dinv = FFTop / D +dinv = np.real(dinv).reshape(nt, nx, ny) + +fig, axs = plt.subplots(2, 2, figsize=(10, 6)) +axs[0][0].imshow(d[:, :, ny // 2], vmin=-20, vmax=20, cmap="bwr") +axs[0][0].set_title("Signal") +axs[0][0].axis("tight") +axs[0][1].imshow( + np.abs(np.fft.fftshift(D.reshape(nfft, nfftk, nfftk), axes=1)[:20, :, nfftk // 2]), + cmap="bwr", +) +axs[0][1].set_title("Fourier Transform 3D mkl fft") +axs[0][1].axis("tight") +axs[1][0].imshow(dinv[:, :, ny // 2], vmin=-20, vmax=20, cmap="bwr") +axs[1][0].set_title("Inverted") +axs[1][0].axis("tight") +axs[1][1].imshow(d[:, :, ny // 2] - dinv[:, :, ny // 2], vmin=-20, vmax=20, cmap="bwr") +axs[1][1].set_title("Error") +axs[1][1].axis("tight") +fig.tight_layout() diff --git a/pylops/signalprocessing/fft.py b/pylops/signalprocessing/fft.py index 64444bcd..50d6dbc2 100644 --- a/pylops/signalprocessing/fft.py +++ b/pylops/signalprocessing/fft.py @@ -7,6 +7,8 @@ import numpy as np import numpy.typing as npt import scipy.fft +from mkl_fft import _numpy_fft +from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift from pylops import LinearOperator from pylops.signalprocessing._baseffts import _BaseFFT, _FFTNorms @@ -365,6 +367,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", + ) -> 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._norm_kwargs = {"norm": None} # equivalent to "backward" in Numpy/Scipy + 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_iffshift(x, axes=self.axis) + if not self.clinear: + x = np.real(x) + if self.real: + y = _numpy_fft.rfft(x, n=self.nfft, axis=self.axis, **self._norm_kwargs) + # Apply scaling to obtain a correct adjoint for this operator + 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 = _numpy_fft.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_fftshift(y, axes=self.axis) + return y + + @reshaped + def _rmatvec(self, x: NDArray) -> NDArray: + if self.fftshift_after: + x = mkl_iffshift(x, axes=self.axis) + if self.real: + # Apply scaling to obtain a correct adjoint for this operator + 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 = _numpy_fft.irfft(x, n=self.nfft, axis=self.axis, **self._norm_kwargs) + else: + y = _numpy_fft.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_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, @@ -395,6 +485,10 @@ def FFT( forward mode, and to :py:func:`scipy.fft.ifft` (or :py:func:`scipy.fft.irfft` for real models) in adjoint mode. + When the mkl_fft engine is chosen, the overloads are of :py:func: 'mkl_fft._numpy_fft.fft' + (or :py:func:`mkl_fft._numpy_fft.rfft` for real models) in forward mode and to :py:func:`mkl_fft._numpy_fft.ifft` + (or :py:func:`mkl_fft._numpy_fft.irfft`for real models) in adjoint mode. + When using ``real=True``, the result of the forward is also multiplied by :math:`\sqrt{2}` for all frequency bins except zero and Nyquist, and the input of the adjoint is multiplied by :math:`1 / \sqrt{2}` for the same frequencies. @@ -452,7 +546,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 arrays. .. note:: Since version 1.17.0, accepts "scipy". @@ -505,7 +599,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 -------- @@ -576,7 +670,19 @@ def FFT( fftshift_after=fftshift_after, dtype=dtype, ) + elif engine == "mkl_fft": + f = _FFT_mklfft( + dims, + axis=axis, + nfft=nfft, + sampling=sampling, + norm=norm, + real=real, + ifftshift_before=ifftshift_before, + fftshift_after=fftshift_after, + dtype=dtype, + ) 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 f54e2972..74e7ed9e 100644 --- a/pylops/signalprocessing/fft2d.py +++ b/pylops/signalprocessing/fft2d.py @@ -6,6 +6,8 @@ import numpy as np import scipy.fft +from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift +from .fftnd import _FFTND_mklfft from pylops import LinearOperator from pylops.signalprocessing._baseffts import _BaseFFTND, _FFTNorms @@ -215,6 +217,102 @@ 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", + ) -> 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.f1, self.f2 = self.fs + del self.fs + + self._norm_kwargs: Dict[str, Union[None, str]] = { + "norm": None + } # equivalent to "backward" in Numpy/Scipy + 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_iffshift(x, axes=self.axes[self.ifftshift_before]) + if not self.clinear: + x = np.real(x) + if self.real: + y = _FFTND_mklfft.rfftn(x, s=self.nffts, axes=self.axes, **self._norm_kwargs) + # 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 = _FFTND_mklfft.fftn(x, s=self.nffts, axes=self.axes, **self._norm_kwargs) + if self.norm is _FFTNorms.ONE_OVER_N: + y *= self._scale + if self.fftshift_after.any(): + y = mkl_fftshift(y, axes=self.axes[self.fftshift_after]) + return y + + @reshaped + def _rmatvec(self, x): + if self.fftshift_after.any(): + x = mkl_iffshift(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 = _FFTND_mklfft.irfftn(x, s=self.nffts, axes=self.axes, **self._norm_kwargs) + else: + y = _FFTND_mklfft.ifftn(x, s=self.nffts, axes=self.axes, **self._norm_kwargs) + if self.norm is _FFTNorms.NONE: + y *= self._scale + y = np.take(y, range(self.dims[self.axes[0]]), axis=self.axes[0]) + y = np.take(y, range(self.dims[self.axes[1]]), axis=self.axes[1]) + if not self.clinear: + y = np.real(y) + if self.ifftshift_before.any(): + y = mkl_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), @@ -242,6 +340,10 @@ def FFT2D( forward mode, and to :py:func:`scipy.fft.ifft2` (or :py:func:`scipy.fft.irfft2` for real models) in adjoint mode. + When the mkl_fft engine is chosen, the overloads are of :py:func: 'mkl_fft._numpy_fft.fft2' + (or :py:func:`mkl_fft._numpy_fft.fft2` for real models) in forward mode and to :py:func:`mkl_fft._numpy_fft.ifft2` + (or :py:func:`mkl_fft._numpy_fft.irfft2`for real models) in adjoint mode. + When using ``real=True``, the result of the forward is also multiplied by :math:`\sqrt{2}` for all frequency bins except zero and Nyquist, and the input of the adjoint is multiplied by :math:`1 / \sqrt{2}` for the same frequencies. @@ -310,7 +412,7 @@ def FFT2D( engine : :obj:`str`, optional .. versionadded:: 1.17.0 - Engine used for fft computation (``numpy`` or ``scipy``). + Engine used for fft computation (``numpy`` or ``scipy`` or ``mkl_fft``). dtype : :obj:`str`, optional Type of elements in input array. Note that the ``dtype`` of the operator is the corresponding complex type even when a real type is provided. @@ -361,7 +463,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 -------- @@ -418,7 +520,19 @@ def FFT2D( fftshift_after=fftshift_after, dtype=dtype, ) + elif engine == "mkl_fft": + 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, + ) 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 a33f4918..6b5182d3 100644 --- a/pylops/signalprocessing/fftnd.py +++ b/pylops/signalprocessing/fftnd.py @@ -6,6 +6,9 @@ import numpy as np import numpy.typing as npt +from mkl_fft import _numpy_fft as pymkl_fft +from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift +from mkl_fft._numpy_fft import _cook_nd_args, asarray from pylops.signalprocessing._baseffts import _BaseFFTND, _FFTNorms from pylops.utils.backend import get_sp_fft @@ -197,6 +200,132 @@ 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", + ) -> 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._norm_kwargs = {"norm": None} # equivalent to "backward" in Numpy/Scipy + if self.norm is _FFTNorms.ORTHO: + self._norm_kwargs["norm"] = "ortho" + self._scale = np.sqrt(1 / np.prod(self.nffts)) + 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_iffshift(x, axes=self.axes[self.ifftshift_before]) + if not self.clinear: + x = np.real(x) + if self.real: + y = self.rfftn(x, s=self.nffts, axes=tuple(self.axes), **self._norm_kwargs) + # 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 = self.fftn(x, s=self.nffts, axes=tuple(self.axes), **self._norm_kwargs) + if self.norm is _FFTNorms.ONE_OVER_N: + y *= self._scale + if self.fftshift_after.any(): + y = mkl_fftshift(y, axes=self.axes[self.fftshift_after]) + return y + + @reshaped + def _rmatvec(self, x: NDArray) -> NDArray: + if self.fftshift_after.any(): + x = mkl_iffshift(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 = self.irfftn(x, s=self.nffts, axes=tuple(self.axes), **self._norm_kwargs) + else: + y = self.ifftn(x, s=self.nffts, axes=tuple(self.axes), **self._norm_kwargs) + 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_fftshift(y, axes=self.axes[self.ifftshift_before]) + return y + + @staticmethod + def rfftn(a, s=None, axes=None, norm=None): + a = asarray(a) + s, axes = _cook_nd_args(a, s, axes) + a = pymkl_fft.rfft(a, s[-1], axes[-1], norm) + for ii in range(len(axes) - 1): + a = pymkl_fft.fft(a, s[ii], axes[ii], norm) + return a + + @staticmethod + def irfftn(a, s=None, axes=None, norm=None): + a = asarray(a) + s, axes = _cook_nd_args(a, s, axes, invreal=1) + for ii in range(len(axes) - 1): + a = pymkl_fft.ifft(a, s[ii], axes[ii], norm) + a = pymkl_fft.irfft(a, s[-1], axes[-1], norm) + return a + + @staticmethod + def fftn(a, s=None, axes=None, norm=None): + a = asarray(a) + s, axes = _cook_nd_args(a, s, axes) + itl = list(range(len(axes))) + itl.reverse() + for ii in itl: + a = pymkl_fft.fft(a, n=s[ii], axis=axes[ii], norm=norm) + return a + + @staticmethod + def ifftn(a, s=None, axes=None, norm=None): + a = asarray(a) + s, axes = _cook_nd_args(a, s, axes) + itl = list(range(len(axes))) + itl.reverse() + for ii in itl: + a = pymkl_fft.ifft(a, n=s[ii], axis=axes[ii], norm=norm) + return a + + 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) + + def FFTND( dims: Union[int, InputDimsLike], axes: Union[int, InputDimsLike] = (-3, -2, -1), @@ -224,6 +353,10 @@ def FFTND( forward mode, and to :py:func:`scipy.fft.ifftn` (or :py:func:`scipy.fft.irfftn` for real models) in adjoint mode. + When the mkl_fft engine is chosen, the overloads are of :py:func: `mkl_fft._numpy_fft.fftn` + (or :py:func:`mkl_fft._numpy_fft.rfftn` for real models) in forward mode and to :py:func:`mkl_fft._numpy_fft.ifftn` + (or :py:func:`mkl_fft._numpy_fft.irfftn`for real models) in adjoint mode. + When using ``real=True``, the result of the forward is also multiplied by :math:`\sqrt{2}` for all frequency bins except zero and Nyquist along the last ``axes``, and the input of the adjoint is multiplied by @@ -297,7 +430,7 @@ def FFTND( engine : :obj:`str`, optional .. versionadded:: 1.17.0 - Engine used for fft computation (``numpy`` or ``scipy``). + Engine used for fft computation (``numpy`` or ``scipy`` or ``mkl_fft``). dtype : :obj:`str`, optional Type of elements in input array. Note that the ``dtype`` of the operator is the corresponding complex type even when a real type is provided. @@ -350,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 ----- @@ -409,6 +542,18 @@ def FFTND( fftshift_after=fftshift_after, dtype=dtype, ) + elif engine == "mkl_fft": + 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, + ) else: raise NotImplementedError("engine must be numpy or scipy") f.name = name diff --git a/pytests/test_ffts.py b/pytests/test_ffts.py index 152f7573..6adfbd4d 100755 --- a/pytests/test_ffts.py +++ b/pytests/test_ffts.py @@ -141,6 +141,59 @@ def _choose_random_axes(ndim, n_choices=2): "dtype": np.complex128, } # nfftnt, complex input, fftw engine +par3t = { + "nt": 41, + "nx": 31, + "ny": 10, + "nfft": None, + "real": True, + "engine": "mkl_fft", + "ifftshift_before": False, + "dtype": np.float64, +} # nfft=nt, real input, fftw engine +par4t = { + "nt": 41, + "nx": 31, + "ny": 10, + "nfft": 64, + "real": True, + "engine": "mkl_fft", + "ifftshift_before": False, + "dtype": np.float32, +} # nfft>nt, real input, fftw engine +par5t = { + "nt": 41, + "nx": 31, + "ny": 10, + "nfft": 16, + "real": False, + "engine": "mkl_fft", + "ifftshift_before": False, + "dtype": np.complex128, +} # nfft=1.21.0 scipy>=1.4.0 torch>=1.2.0 From 300b9c51ee126c70fdd23908ec1b273719f05c79 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Tue, 7 Feb 2023 20:17:46 +0530 Subject: [PATCH 02/48] Added mkl-fft along with examples and tests --- pytests/test_ffts.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/pytests/test_ffts.py b/pytests/test_ffts.py index 6adfbd4d..f9083a2c 100755 --- a/pytests/test_ffts.py +++ b/pytests/test_ffts.py @@ -215,7 +215,7 @@ def test_unknown_engine(par): (np.float16, 1), (np.float32, 4), (np.float64, 11), - # (np.float128, 11), + (np.float128, 11), ], norm=["ortho", "none", "1/n"], ifftshift_before=[False, True], @@ -287,7 +287,7 @@ def test_FFT_small_real(par): (np.float16, 1), (np.float32, 3), (np.float64, 11), - # (np.float128, 11), + (np.float128, 11), ], ifftshift_before=[False, True], engine=["numpy", "fftw", "scipy", "mkl_fft"], @@ -333,7 +333,7 @@ def test_FFT_random_real(par): par_lists_fft_small_cpx = dict( - dtype_precision=[(np.complex64, 4), (np.complex128, 11)], + dtype_precision=[(np.complex64, 4), (np.complex128, 11), (np.complex256, 11)], norm=["ortho", "none", "1/n"], ifftshift_before=[False, True], fftshift_after=[False, True], @@ -397,10 +397,10 @@ def test_FFT_small_complex(par): (np.float16, 1), (np.float32, 3), (np.float64, 11), - # (np.float128, 11), + (np.float128, 11), (np.complex64, 3), (np.complex128, 11), - # (np.complex256, 11), + (np.complex256, 11), ], ifftshift_before=[False, True], fftshift_after=[False, True], @@ -479,7 +479,7 @@ def test_FFT_random_complex(par): (np.float16, 1), (np.float32, 3), (np.float64, 11), - # (np.float128, 11), + (np.float128, 11), ], ifftshift_before=[False, True], engine=["numpy", "scipy", "mkl_fft"], @@ -537,10 +537,10 @@ def test_FFT2D_random_real(par): (np.float16, 1), (np.float32, 3), (np.float64, 11), - # (np.float128, 11), + (np.float128, 11), (np.complex64, 3), (np.complex128, 11), - # (np.complex256, 11), + (np.complex256, 11), ], ifftshift_before=itertools.product([False, True], [False, True]), fftshift_after=itertools.product([False, True], [False, True]), @@ -619,7 +619,7 @@ def test_FFT2D_random_complex(par): (np.float16, 1), (np.float32, 3), (np.float64, 11), - # (np.float128, 11), + (np.float128, 11), ], engine=["numpy", "scipy", "mkl_fft"], ) @@ -678,10 +678,10 @@ def test_FFTND_random_real(par): (np.float16, 1), (np.float32, 3), (np.float64, 11), - # (np.float128, 11), + (np.float128, 11), (np.complex64, 3), (np.complex128, 11), - # (np.complex256, 11), + (np.complex256, 11), ], engine=["numpy", "scipy", "mkl_fft"], ) @@ -753,7 +753,7 @@ def test_FFTND_random_complex(par): par_lists_fft2dnd_small_cpx = dict( - dtype_precision=[(np.complex64, 5), (np.complex128, 11)], + dtype_precision=[(np.complex64, 5), (np.complex128, 11), (np.complex256, 11)], norm=["ortho", "none", "1/n"], engine=["numpy", "scipy", "mkl_fft"], ) From 1ca769ed4e9ad5023b5dd5ef30ab1772579e16ba Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Tue, 7 Feb 2023 20:45:50 +0530 Subject: [PATCH 03/48] Renamed _numpy_fft to pymkl_fft --- pylops/signalprocessing/fft.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/pylops/signalprocessing/fft.py b/pylops/signalprocessing/fft.py index 50d6dbc2..a98172d3 100644 --- a/pylops/signalprocessing/fft.py +++ b/pylops/signalprocessing/fft.py @@ -7,7 +7,7 @@ import numpy as np import numpy.typing as npt import scipy.fft -from mkl_fft import _numpy_fft +from mkl_fft import _numpy_fft as pymkl_fft from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift from pylops import LinearOperator @@ -409,13 +409,13 @@ def _matvec(self, x: NDArray) -> NDArray: if not self.clinear: x = np.real(x) if self.real: - y = _numpy_fft.rfft(x, n=self.nfft, axis=self.axis, **self._norm_kwargs) + y = pymkl_fft.rfft(x, n=self.nfft, axis=self.axis, **self._norm_kwargs) # Apply scaling to obtain a correct adjoint for this operator 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 = _numpy_fft.fft(x, n=self.nfft, axis=self.axis, **self._norm_kwargs) + y = pymkl_fft.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: @@ -432,9 +432,9 @@ def _rmatvec(self, x: NDArray) -> NDArray: 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 = _numpy_fft.irfft(x, n=self.nfft, axis=self.axis, **self._norm_kwargs) + y = pymkl_fft.irfft(x, n=self.nfft, axis=self.axis, **self._norm_kwargs) else: - y = _numpy_fft.ifft(x, n=self.nfft, axis=self.axis, **self._norm_kwargs) + y = pymkl_fft.ifft(x, n=self.nfft, axis=self.axis, **self._norm_kwargs) if self.norm is _FFTNorms.NONE: y *= self._scale From b002cb03edc625362a049984ca3ebf9251e34033 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Tue, 7 Feb 2023 21:00:05 +0530 Subject: [PATCH 04/48] Updated requirements-doc.txt --- requirements-doc.txt | 3 +++ 1 file changed, 3 insertions(+) diff --git a/requirements-doc.txt b/requirements-doc.txt index d137b621..fd3f56f6 100644 --- a/requirements-doc.txt +++ b/requirements-doc.txt @@ -1,3 +1,6 @@ +mkl +mkl-fft +mkl-service numpy>=1.21.0 scipy>=1.4.0 torch From 1c557bfa5ab494baa3a36efe8e98cfc16cf14496 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Wed, 8 Feb 2023 16:58:12 +0530 Subject: [PATCH 05/48] Quick Fix --- pylops/signalprocessing/fft.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/pylops/signalprocessing/fft.py b/pylops/signalprocessing/fft.py index a98172d3..09acd1df 100644 --- a/pylops/signalprocessing/fft.py +++ b/pylops/signalprocessing/fft.py @@ -1,3 +1,7 @@ +import os + +os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE" + __all__ = ["FFT"] import logging From cf62f9a7a33bd6e75c85be43ade2e6e658734ad7 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Wed, 8 Feb 2023 17:13:54 +0530 Subject: [PATCH 06/48] Test fix --- pylops/signalprocessing/fft.py | 6 ++++-- pylops/signalprocessing/fftnd.py | 23 ++++++++++++++--------- 2 files changed, 18 insertions(+), 11 deletions(-) diff --git a/pylops/signalprocessing/fft.py b/pylops/signalprocessing/fft.py index 09acd1df..f2e66379 100644 --- a/pylops/signalprocessing/fft.py +++ b/pylops/signalprocessing/fft.py @@ -11,8 +11,6 @@ import numpy as np import numpy.typing as npt import scipy.fft -from mkl_fft import _numpy_fft as pymkl_fft -from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift from pylops import LinearOperator from pylops.signalprocessing._baseffts import _BaseFFT, _FFTNorms @@ -408,6 +406,8 @@ def __init__( @reshaped def _matvec(self, x: NDArray) -> NDArray: + from mkl_fft import _numpy_fft as pymkl_fft + from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift if self.ifftshift_before: x = mkl_iffshift(x, axes=self.axis) if not self.clinear: @@ -428,6 +428,8 @@ def _matvec(self, x: NDArray) -> NDArray: @reshaped def _rmatvec(self, x: NDArray) -> NDArray: + from mkl_fft import _numpy_fft as pymkl_fft + from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift if self.fftshift_after: x = mkl_iffshift(x, axes=self.axis) if self.real: diff --git a/pylops/signalprocessing/fftnd.py b/pylops/signalprocessing/fftnd.py index 6b5182d3..ce29eed1 100644 --- a/pylops/signalprocessing/fftnd.py +++ b/pylops/signalprocessing/fftnd.py @@ -6,9 +6,6 @@ import numpy as np import numpy.typing as npt -from mkl_fft import _numpy_fft as pymkl_fft -from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift -from mkl_fft._numpy_fft import _cook_nd_args, asarray from pylops.signalprocessing._baseffts import _BaseFFTND, _FFTNorms from pylops.utils.backend import get_sp_fft @@ -238,6 +235,8 @@ def __init__( @reshaped def _matvec(self, x: NDArray) -> NDArray: + from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift + if self.ifftshift_before.any(): x = mkl_iffshift(x, axes=self.axes[self.ifftshift_before]) if not self.clinear: @@ -258,6 +257,8 @@ def _matvec(self, x: NDArray) -> NDArray: @reshaped def _rmatvec(self, x: NDArray) -> NDArray: + from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift + if self.fftshift_after.any(): x = mkl_iffshift(x, axes=self.axes[self.fftshift_after]) if self.real: @@ -284,40 +285,44 @@ def _rmatvec(self, x: NDArray) -> NDArray: @staticmethod def rfftn(a, s=None, axes=None, norm=None): + from mkl_fft._numpy_fft import rfft, fft, _cook_nd_args, asarray a = asarray(a) s, axes = _cook_nd_args(a, s, axes) - a = pymkl_fft.rfft(a, s[-1], axes[-1], norm) + a = rfft(a, s[-1], axes[-1], norm) for ii in range(len(axes) - 1): - a = pymkl_fft.fft(a, s[ii], axes[ii], norm) + a = fft(a, s[ii], axes[ii], norm) return a @staticmethod def irfftn(a, s=None, axes=None, norm=None): + from mkl_fft._numpy_fft import irfft, ifft, asarray, _cook_nd_args a = asarray(a) s, axes = _cook_nd_args(a, s, axes, invreal=1) for ii in range(len(axes) - 1): - a = pymkl_fft.ifft(a, s[ii], axes[ii], norm) - a = pymkl_fft.irfft(a, s[-1], axes[-1], norm) + a = ifft(a, s[ii], axes[ii], norm) + a = irfft(a, s[-1], axes[-1], norm) return a @staticmethod def fftn(a, s=None, axes=None, norm=None): + from mkl_fft._numpy_fft import fft, asarray, _cook_nd_args a = asarray(a) s, axes = _cook_nd_args(a, s, axes) itl = list(range(len(axes))) itl.reverse() for ii in itl: - a = pymkl_fft.fft(a, n=s[ii], axis=axes[ii], norm=norm) + a = fft(a, n=s[ii], axis=axes[ii], norm=norm) return a @staticmethod def ifftn(a, s=None, axes=None, norm=None): + from mkl_fft._numpy_fft import ifft, asarray, _cook_nd_args a = asarray(a) s, axes = _cook_nd_args(a, s, axes) itl = list(range(len(axes))) itl.reverse() for ii in itl: - a = pymkl_fft.ifft(a, n=s[ii], axis=axes[ii], norm=norm) + a = ifft(a, n=s[ii], axis=axes[ii], norm=norm) return a def __truediv__(self, y: npt.ArrayLike) -> npt.ArrayLike: From da2a621b0dcec1b17370ba2b27b392374dd8aba1 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Wed, 8 Feb 2023 17:24:21 +0530 Subject: [PATCH 07/48] Removed KMP duplication --- pylops/signalprocessing/fft.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/pylops/signalprocessing/fft.py b/pylops/signalprocessing/fft.py index f2e66379..d36a5b9d 100644 --- a/pylops/signalprocessing/fft.py +++ b/pylops/signalprocessing/fft.py @@ -1,7 +1,3 @@ -import os - -os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE" - __all__ = ["FFT"] import logging From 4868d3736fb78c0e3e1d23d7509ce80742cdb4ce Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Wed, 8 Feb 2023 22:19:49 +0530 Subject: [PATCH 08/48] Added nomkl to env-dev --- environment-dev.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/environment-dev.yml b/environment-dev.yml index 99687625..bc4ed0d7 100755 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -27,6 +27,7 @@ dependencies: - autopep8 - isort - black + - nomkl - pip: - devito - scikit-fmm From 2fa67327a6f3f16d6e352b50f91a0a1e1b87d5ec Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Thu, 9 Feb 2023 17:43:08 +0530 Subject: [PATCH 09/48] FFT2D Fix --- pylops/signalprocessing/fft2d.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pylops/signalprocessing/fft2d.py b/pylops/signalprocessing/fft2d.py index 74e7ed9e..57595428 100644 --- a/pylops/signalprocessing/fft2d.py +++ b/pylops/signalprocessing/fft2d.py @@ -6,7 +6,6 @@ import numpy as np import scipy.fft -from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift from .fftnd import _FFTND_mklfft from pylops import LinearOperator @@ -266,6 +265,7 @@ def __init__( @reshaped def _matvec(self, x): + from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift if self.ifftshift_before.any(): x = mkl_iffshift(x, axes=self.axes[self.ifftshift_before]) if not self.clinear: @@ -286,6 +286,7 @@ def _matvec(self, x): @reshaped def _rmatvec(self, x): + from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift if self.fftshift_after.any(): x = mkl_iffshift(x, axes=self.axes[self.fftshift_after]) if self.real: From 266ca7ab0f6561ff6730e9d1cd8a54900971163f Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Thu, 9 Feb 2023 18:20:44 +0530 Subject: [PATCH 10/48] FFT2D updated --- pylops/signalprocessing/fft2d.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pylops/signalprocessing/fft2d.py b/pylops/signalprocessing/fft2d.py index 57595428..134f7778 100644 --- a/pylops/signalprocessing/fft2d.py +++ b/pylops/signalprocessing/fft2d.py @@ -6,7 +6,7 @@ import numpy as np import scipy.fft -from .fftnd import _FFTND_mklfft +from pylops.signalprocessing.fftnd import _FFTND_mklfft from pylops import LinearOperator from pylops.signalprocessing._baseffts import _BaseFFTND, _FFTNorms From 73c48a2e96f16c65b67113409e821b4218164bed Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Thu, 9 Feb 2023 19:11:41 +0530 Subject: [PATCH 11/48] Removed mkl scipy backend --- pylops/signalprocessing/fft.py | 12 ++++++------ pylops/signalprocessing/fft2d.py | 12 ++++++------ pylops/signalprocessing/fftnd.py | 14 ++++++-------- 3 files changed, 18 insertions(+), 20 deletions(-) diff --git a/pylops/signalprocessing/fft.py b/pylops/signalprocessing/fft.py index d36a5b9d..006631ca 100644 --- a/pylops/signalprocessing/fft.py +++ b/pylops/signalprocessing/fft.py @@ -403,9 +403,9 @@ def __init__( @reshaped def _matvec(self, x: NDArray) -> NDArray: from mkl_fft import _numpy_fft as pymkl_fft - from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift + # from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift if self.ifftshift_before: - x = mkl_iffshift(x, axes=self.axis) + x = np.fft.ifftshift(x, axes=self.axis) if not self.clinear: x = np.real(x) if self.real: @@ -419,15 +419,15 @@ def _matvec(self, x: NDArray) -> NDArray: if self.norm is _FFTNorms.ONE_OVER_N: y *= self._scale if self.fftshift_after: - y = mkl_fftshift(y, axes=self.axis) + y = np.fft.fftshift(y, axes=self.axis) return y @reshaped def _rmatvec(self, x: NDArray) -> NDArray: from mkl_fft import _numpy_fft as pymkl_fft - from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift + # from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift if self.fftshift_after: - x = mkl_iffshift(x, axes=self.axis) + x = np.fft.ifftshift(x, axes=self.axis) if self.real: # Apply scaling to obtain a correct adjoint for this operator x = x.copy() @@ -448,7 +448,7 @@ def _rmatvec(self, x: NDArray) -> NDArray: if not self.clinear: y = np.real(y) if self.ifftshift_before: - y = mkl_fftshift(y, axes=self.axis) + y = np.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 134f7778..cf150994 100644 --- a/pylops/signalprocessing/fft2d.py +++ b/pylops/signalprocessing/fft2d.py @@ -265,9 +265,9 @@ def __init__( @reshaped def _matvec(self, x): - from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift + # from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift if self.ifftshift_before.any(): - x = mkl_iffshift(x, axes=self.axes[self.ifftshift_before]) + x = np.fft.ifftshift(x, axes=self.axes[self.ifftshift_before]) if not self.clinear: x = np.real(x) if self.real: @@ -281,14 +281,14 @@ def _matvec(self, x): if self.norm is _FFTNorms.ONE_OVER_N: y *= self._scale if self.fftshift_after.any(): - y = mkl_fftshift(y, axes=self.axes[self.fftshift_after]) + y = np.fft.fftshift(y, axes=self.axes[self.fftshift_after]) return y @reshaped def _rmatvec(self, x): - from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift + # from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift if self.fftshift_after.any(): - x = mkl_iffshift(x, axes=self.axes[self.fftshift_after]) + x = np.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() @@ -305,7 +305,7 @@ def _rmatvec(self, x): if not self.clinear: y = np.real(y) if self.ifftshift_before.any(): - y = mkl_fftshift(y, axes=self.axes[self.ifftshift_before]) + y = np.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 ce29eed1..9ce04af7 100644 --- a/pylops/signalprocessing/fftnd.py +++ b/pylops/signalprocessing/fftnd.py @@ -235,10 +235,9 @@ def __init__( @reshaped def _matvec(self, x: NDArray) -> NDArray: - from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift - + # from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift if self.ifftshift_before.any(): - x = mkl_iffshift(x, axes=self.axes[self.ifftshift_before]) + x = np.fft.ifftshift(x, axes=self.axes[self.ifftshift_before]) if not self.clinear: x = np.real(x) if self.real: @@ -252,15 +251,14 @@ def _matvec(self, x: NDArray) -> NDArray: if self.norm is _FFTNorms.ONE_OVER_N: y *= self._scale if self.fftshift_after.any(): - y = mkl_fftshift(y, axes=self.axes[self.fftshift_after]) + y = np.fft.fftshift(y, axes=self.axes[self.fftshift_after]) return y @reshaped def _rmatvec(self, x: NDArray) -> NDArray: - from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift - + # from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift if self.fftshift_after.any(): - x = mkl_iffshift(x, axes=self.axes[self.fftshift_after]) + x = np.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() @@ -280,7 +278,7 @@ def _rmatvec(self, x: NDArray) -> NDArray: if not self.clinear: y = np.real(y) if self.ifftshift_before.any(): - y = mkl_fftshift(y, axes=self.axes[self.ifftshift_before]) + y = np.fft.fftshift(y, axes=self.axes[self.ifftshift_before]) return y @staticmethod From 8e6cb9a695f1c54c08868a1b357f289a13190282 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Thu, 9 Feb 2023 19:24:15 +0530 Subject: [PATCH 12/48] Removed nomkl --- environment-dev.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/environment-dev.yml b/environment-dev.yml index bc4ed0d7..99687625 100755 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -27,7 +27,6 @@ dependencies: - autopep8 - isort - black - - nomkl - pip: - devito - scikit-fmm From 5c05501d44487505f487bb3dbf69ea77944a0b76 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Sun, 12 Feb 2023 21:56:22 +0530 Subject: [PATCH 13/48] Added imports above numpy --- pylops/signalprocessing/fft.py | 15 +++++++-------- pylops/signalprocessing/fft2d.py | 13 ++++++------- pylops/signalprocessing/fftnd.py | 24 +++++++++++++++--------- 3 files changed, 28 insertions(+), 24 deletions(-) diff --git a/pylops/signalprocessing/fft.py b/pylops/signalprocessing/fft.py index 006631ca..7a3051b3 100644 --- a/pylops/signalprocessing/fft.py +++ b/pylops/signalprocessing/fft.py @@ -3,6 +3,9 @@ import logging import warnings from typing import Optional, Union +from mkl_fft import _numpy_fft as pymkl_fft +from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_ifftshift + import numpy as np import numpy.typing as npt @@ -402,10 +405,8 @@ def __init__( @reshaped def _matvec(self, x: NDArray) -> NDArray: - from mkl_fft import _numpy_fft as pymkl_fft - # from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift if self.ifftshift_before: - x = np.fft.ifftshift(x, axes=self.axis) + x = mkl_ifftshift(x, axes=self.axis) if not self.clinear: x = np.real(x) if self.real: @@ -419,15 +420,13 @@ def _matvec(self, x: NDArray) -> NDArray: if self.norm is _FFTNorms.ONE_OVER_N: y *= self._scale if self.fftshift_after: - y = np.fft.fftshift(y, axes=self.axis) + y = mkl_fftshift(y, axes=self.axis) return y @reshaped def _rmatvec(self, x: NDArray) -> NDArray: - from mkl_fft import _numpy_fft as pymkl_fft - # from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift if self.fftshift_after: - x = np.fft.ifftshift(x, axes=self.axis) + x = mkl_ifftshift(x, axes=self.axis) if self.real: # Apply scaling to obtain a correct adjoint for this operator x = x.copy() @@ -448,7 +447,7 @@ def _rmatvec(self, x: NDArray) -> NDArray: if not self.clinear: y = np.real(y) if self.ifftshift_before: - y = np.fft.fftshift(y, axes=self.axis) + y = mkl_fftshift(y, axes=self.axis) return y def __truediv__(self, y): diff --git a/pylops/signalprocessing/fft2d.py b/pylops/signalprocessing/fft2d.py index cf150994..e1957b25 100644 --- a/pylops/signalprocessing/fft2d.py +++ b/pylops/signalprocessing/fft2d.py @@ -4,9 +4,10 @@ import warnings from typing import Dict, Optional, Sequence, Union +from pylops.signalprocessing.fftnd import _FFTND_mklfft +from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_ifftshift import numpy as np import scipy.fft -from pylops.signalprocessing.fftnd import _FFTND_mklfft from pylops import LinearOperator from pylops.signalprocessing._baseffts import _BaseFFTND, _FFTNorms @@ -265,9 +266,8 @@ def __init__( @reshaped def _matvec(self, x): - # from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift if self.ifftshift_before.any(): - x = np.fft.ifftshift(x, axes=self.axes[self.ifftshift_before]) + x = mkl_ifftshift(x, axes=self.axes[self.ifftshift_before]) if not self.clinear: x = np.real(x) if self.real: @@ -281,14 +281,13 @@ def _matvec(self, x): if self.norm is _FFTNorms.ONE_OVER_N: y *= self._scale if self.fftshift_after.any(): - y = np.fft.fftshift(y, axes=self.axes[self.fftshift_after]) + y = mkl_fftshift(y, axes=self.axes[self.fftshift_after]) return y @reshaped def _rmatvec(self, x): - # from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift if self.fftshift_after.any(): - x = np.fft.ifftshift(x, axes=self.axes[self.fftshift_after]) + x = mkl_ifftshift(x, axes=self.axes[self.fftshift_after]) if self.real: # Apply scaling to obtain a correct adjoint for this operator x = x.copy() @@ -305,7 +304,7 @@ def _rmatvec(self, x): if not self.clinear: y = np.real(y) if self.ifftshift_before.any(): - y = np.fft.fftshift(y, axes=self.axes[self.ifftshift_before]) + y = mkl_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 9ce04af7..ae00b6e1 100644 --- a/pylops/signalprocessing/fftnd.py +++ b/pylops/signalprocessing/fftnd.py @@ -235,9 +235,10 @@ def __init__( @reshaped def _matvec(self, x: NDArray) -> NDArray: - # from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift + from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift + if self.ifftshift_before.any(): - x = np.fft.ifftshift(x, axes=self.axes[self.ifftshift_before]) + x = mkl_fftshift(x, axes=self.axes[self.ifftshift_before]) if not self.clinear: x = np.real(x) if self.real: @@ -251,14 +252,15 @@ def _matvec(self, x: NDArray) -> NDArray: if self.norm is _FFTNorms.ONE_OVER_N: y *= self._scale if self.fftshift_after.any(): - y = np.fft.fftshift(y, axes=self.axes[self.fftshift_after]) + y = mkl_iffshift(y, axes=self.axes[self.fftshift_after]) return y @reshaped def _rmatvec(self, x: NDArray) -> NDArray: - # from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift + from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift + if self.fftshift_after.any(): - x = np.fft.ifftshift(x, axes=self.axes[self.fftshift_after]) + x = mkl_iffshift(x, axes=self.axes[self.fftshift_after]) if self.real: # Apply scaling to obtain a correct adjoint for this operator x = x.copy() @@ -278,12 +280,13 @@ def _rmatvec(self, x: NDArray) -> NDArray: if not self.clinear: y = np.real(y) if self.ifftshift_before.any(): - y = np.fft.fftshift(y, axes=self.axes[self.ifftshift_before]) + y = mkl_fftshift(y, axes=self.axes[self.ifftshift_before]) return y @staticmethod def rfftn(a, s=None, axes=None, norm=None): - from mkl_fft._numpy_fft import rfft, fft, _cook_nd_args, asarray + from mkl_fft._numpy_fft import asarray, _cook_nd_args, fft, rfft + a = asarray(a) s, axes = _cook_nd_args(a, s, axes) a = rfft(a, s[-1], axes[-1], norm) @@ -293,7 +296,8 @@ def rfftn(a, s=None, axes=None, norm=None): @staticmethod def irfftn(a, s=None, axes=None, norm=None): - from mkl_fft._numpy_fft import irfft, ifft, asarray, _cook_nd_args + from mkl_fft._numpy_fft import ifft, asarray, _cook_nd_args, irfft + a = asarray(a) s, axes = _cook_nd_args(a, s, axes, invreal=1) for ii in range(len(axes) - 1): @@ -303,7 +307,8 @@ def irfftn(a, s=None, axes=None, norm=None): @staticmethod def fftn(a, s=None, axes=None, norm=None): - from mkl_fft._numpy_fft import fft, asarray, _cook_nd_args + from mkl_fft._numpy_fft import asarray, _cook_nd_args, fft + a = asarray(a) s, axes = _cook_nd_args(a, s, axes) itl = list(range(len(axes))) @@ -315,6 +320,7 @@ def fftn(a, s=None, axes=None, norm=None): @staticmethod def ifftn(a, s=None, axes=None, norm=None): from mkl_fft._numpy_fft import ifft, asarray, _cook_nd_args + a = asarray(a) s, axes = _cook_nd_args(a, s, axes) itl = list(range(len(axes))) From 3f221fc7d41c0589e5d51ee3031e4074a721911f Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Mon, 13 Feb 2023 00:03:11 +0530 Subject: [PATCH 14/48] Updated ffts --- pylops/signalprocessing/fft.py | 9 ++++++--- pylops/signalprocessing/fft2d.py | 5 ++++- pylops/signalprocessing/fftnd.py | 4 ++-- 3 files changed, 12 insertions(+), 6 deletions(-) diff --git a/pylops/signalprocessing/fft.py b/pylops/signalprocessing/fft.py index 7a3051b3..9cf2cd72 100644 --- a/pylops/signalprocessing/fft.py +++ b/pylops/signalprocessing/fft.py @@ -3,9 +3,6 @@ import logging import warnings from typing import Optional, Union -from mkl_fft import _numpy_fft as pymkl_fft -from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_ifftshift - import numpy as np import numpy.typing as npt @@ -405,6 +402,9 @@ def __init__( @reshaped def _matvec(self, x: NDArray) -> NDArray: + from mkl_fft import _numpy_fft as pymkl_fft + from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_ifftshift + if self.ifftshift_before: x = mkl_ifftshift(x, axes=self.axis) if not self.clinear: @@ -425,6 +425,9 @@ def _matvec(self, x: NDArray) -> NDArray: @reshaped def _rmatvec(self, x: NDArray) -> NDArray: + from mkl_fft import _numpy_fft as pymkl_fft + from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_ifftshift + if self.fftshift_after: x = mkl_ifftshift(x, axes=self.axis) if self.real: diff --git a/pylops/signalprocessing/fft2d.py b/pylops/signalprocessing/fft2d.py index e1957b25..8cf39470 100644 --- a/pylops/signalprocessing/fft2d.py +++ b/pylops/signalprocessing/fft2d.py @@ -5,7 +5,6 @@ from typing import Dict, Optional, Sequence, Union from pylops.signalprocessing.fftnd import _FFTND_mklfft -from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_ifftshift import numpy as np import scipy.fft @@ -266,6 +265,8 @@ def __init__( @reshaped def _matvec(self, x): + from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_ifftshift + if self.ifftshift_before.any(): x = mkl_ifftshift(x, axes=self.axes[self.ifftshift_before]) if not self.clinear: @@ -286,6 +287,8 @@ def _matvec(self, x): @reshaped def _rmatvec(self, x): + from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_ifftshift + if self.fftshift_after.any(): x = mkl_ifftshift(x, axes=self.axes[self.fftshift_after]) if self.real: diff --git a/pylops/signalprocessing/fftnd.py b/pylops/signalprocessing/fftnd.py index ae00b6e1..095ce672 100644 --- a/pylops/signalprocessing/fftnd.py +++ b/pylops/signalprocessing/fftnd.py @@ -238,7 +238,7 @@ def _matvec(self, x: NDArray) -> NDArray: from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift if self.ifftshift_before.any(): - x = mkl_fftshift(x, axes=self.axes[self.ifftshift_before]) + x = mkl_iffshift(x, axes=self.axes[self.ifftshift_before]) if not self.clinear: x = np.real(x) if self.real: @@ -252,7 +252,7 @@ def _matvec(self, x: NDArray) -> NDArray: if self.norm is _FFTNorms.ONE_OVER_N: y *= self._scale if self.fftshift_after.any(): - y = mkl_iffshift(y, axes=self.axes[self.fftshift_after]) + y = mkl_fftshift(y, axes=self.axes[self.fftshift_after]) return y @reshaped From 4ee0f71d5dbb10aaee91711a591dfb866e820c6a Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Wed, 15 Feb 2023 21:50:51 +0530 Subject: [PATCH 15/48] Added KMP_DUPLICATE_LIB_OK --- azure-pipelines.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/azure-pipelines.yml b/azure-pipelines.yml index d778e70a..16d56018 100644 --- a/azure-pipelines.yml +++ b/azure-pipelines.yml @@ -51,6 +51,7 @@ jobs: variables: NUMBA_NUM_THREADS: 1 + KMP_DUPLICATE_LIB_OK: TRUE steps: - task: UsePythonVersion@0 From 6101ba4275d5da21f90f572802786dbb52a63310 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Wed, 15 Feb 2023 21:51:26 +0530 Subject: [PATCH 16/48] Added KMP_DUPLICATE_LIB_OK --- azure-pipelines.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/azure-pipelines.yml b/azure-pipelines.yml index 16d56018..dae33c1d 100644 --- a/azure-pipelines.yml +++ b/azure-pipelines.yml @@ -81,6 +81,8 @@ jobs: variables: NUMBA_NUM_THREADS: 1 + KMP_DUPLICATE_LIB_OK: TRUE + steps: - task: UsePythonVersion@0 From fb535f54ead8c6802018a45df41c6e53c4fb3222 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Wed, 15 Feb 2023 22:58:53 +0530 Subject: [PATCH 17/48] Add MKL_THREADING --- azure-pipelines.yml | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/azure-pipelines.yml b/azure-pipelines.yml index dae33c1d..58eac648 100644 --- a/azure-pipelines.yml +++ b/azure-pipelines.yml @@ -52,6 +52,9 @@ jobs: variables: NUMBA_NUM_THREADS: 1 KMP_DUPLICATE_LIB_OK: TRUE + MKL_THREADING_LAYER: GNU + MKL_INTERFACE_LAYER: GNU + OMP_NESTED: FALSE steps: - task: UsePythonVersion@0 @@ -81,8 +84,9 @@ jobs: variables: NUMBA_NUM_THREADS: 1 - KMP_DUPLICATE_LIB_OK: TRUE - + MKL_THREADING_LAYER: GNU + MKL_INTERFACE_LAYER: GNU + OMP_NESTED: FALSE steps: - task: UsePythonVersion@0 From c27728adbcdb01a90c160b73382479fd39071497 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Thu, 16 Feb 2023 11:59:51 +0530 Subject: [PATCH 18/48] Updated plot_fft.py --- examples/plot_fft.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/plot_fft.py b/examples/plot_fft.py index 6bb6b01a..e804a4f6 100644 --- a/examples/plot_fft.py +++ b/examples/plot_fft.py @@ -6,11 +6,11 @@ and :py:class:`pylops.signalprocessing.FFTND` operators to apply the Fourier Transform to the model and the inverse Fourier Transform to the data. """ +import pylops + import matplotlib.pyplot as plt import numpy as np -import pylops - plt.close("all") ############################################################################### From 4911d0ea0e528bcfdddcd529bbf6e888f4406bc5 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Thu, 16 Feb 2023 13:04:12 +0530 Subject: [PATCH 19/48] Added torch fix --- pylops/signalprocessing/fft.py | 1 + pylops/signalprocessing/fft2d.py | 1 + pylops/signalprocessing/fftnd.py | 1 + 3 files changed, 3 insertions(+) diff --git a/pylops/signalprocessing/fft.py b/pylops/signalprocessing/fft.py index 9cf2cd72..d14bb3f5 100644 --- a/pylops/signalprocessing/fft.py +++ b/pylops/signalprocessing/fft.py @@ -1,5 +1,6 @@ __all__ = ["FFT"] +import torch import logging import warnings from typing import Optional, Union diff --git a/pylops/signalprocessing/fft2d.py b/pylops/signalprocessing/fft2d.py index 8cf39470..4641601e 100644 --- a/pylops/signalprocessing/fft2d.py +++ b/pylops/signalprocessing/fft2d.py @@ -1,5 +1,6 @@ __all__ = ["FFT2D"] +import torch import logging import warnings from typing import Dict, Optional, Sequence, Union diff --git a/pylops/signalprocessing/fftnd.py b/pylops/signalprocessing/fftnd.py index 095ce672..e407f7e6 100644 --- a/pylops/signalprocessing/fftnd.py +++ b/pylops/signalprocessing/fftnd.py @@ -1,5 +1,6 @@ __all__ = ["FFTND"] +import torch import logging import warnings from typing import Optional, Sequence, Union From acfe2ce7fea857ac8387a01d1de8f32e868957c8 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Thu, 16 Feb 2023 15:27:16 +0530 Subject: [PATCH 20/48] Added deps and fft, fft2d, fftnd --- azure-pipelines.yml | 7 ------- examples/plot_fft.py | 4 ++-- pylops/signalprocessing/fft.py | 28 ++++++++++++------------- pylops/signalprocessing/fft2d.py | 24 ++++++++++----------- pylops/signalprocessing/fftnd.py | 36 +++++++++++++------------------- pylops/utils/deps.py | 19 +++++++++++++++++ 6 files changed, 61 insertions(+), 57 deletions(-) diff --git a/azure-pipelines.yml b/azure-pipelines.yml index 58eac648..d778e70a 100644 --- a/azure-pipelines.yml +++ b/azure-pipelines.yml @@ -51,10 +51,6 @@ jobs: variables: NUMBA_NUM_THREADS: 1 - KMP_DUPLICATE_LIB_OK: TRUE - MKL_THREADING_LAYER: GNU - MKL_INTERFACE_LAYER: GNU - OMP_NESTED: FALSE steps: - task: UsePythonVersion@0 @@ -84,9 +80,6 @@ jobs: variables: NUMBA_NUM_THREADS: 1 - MKL_THREADING_LAYER: GNU - MKL_INTERFACE_LAYER: GNU - OMP_NESTED: FALSE steps: - task: UsePythonVersion@0 diff --git a/examples/plot_fft.py b/examples/plot_fft.py index e804a4f6..6bb6b01a 100644 --- a/examples/plot_fft.py +++ b/examples/plot_fft.py @@ -6,11 +6,11 @@ and :py:class:`pylops.signalprocessing.FFTND` operators to apply the Fourier Transform to the model and the inverse Fourier Transform to the data. """ -import pylops - import matplotlib.pyplot as plt import numpy as np +import pylops + plt.close("all") ############################################################################### diff --git a/pylops/signalprocessing/fft.py b/pylops/signalprocessing/fft.py index d14bb3f5..5210c856 100644 --- a/pylops/signalprocessing/fft.py +++ b/pylops/signalprocessing/fft.py @@ -1,6 +1,5 @@ __all__ = ["FFT"] -import torch import logging import warnings from typing import Optional, Union @@ -16,10 +15,15 @@ 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._numpy_fft as pymkl_fft + from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_ifftshift + logging.basicConfig(format="%(levelname)s: %(message)s", level=logging.WARNING) @@ -403,9 +407,6 @@ def __init__( @reshaped def _matvec(self, x: NDArray) -> NDArray: - from mkl_fft import _numpy_fft as pymkl_fft - from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_ifftshift - if self.ifftshift_before: x = mkl_ifftshift(x, axes=self.axis) if not self.clinear: @@ -426,9 +427,6 @@ def _matvec(self, x: NDArray) -> NDArray: @reshaped def _rmatvec(self, x: NDArray) -> NDArray: - from mkl_fft import _numpy_fft as pymkl_fft - from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_ifftshift - if self.fftshift_after: x = mkl_ifftshift(x, axes=self.axis) if self.real: @@ -649,10 +647,8 @@ def FFT( dtype=dtype, **kwargs_fftw, ) - elif engine == "numpy" or (engine == "fftw" and pyfftw_message is not None): - if engine == "fftw" and pyfftw_message is not None: - logging.warning(pyfftw_message) - f = _FFT_numpy( + elif engine == "mkl_fft" and mkl_fft_message is None: + f = _FFT_mklfft( dims, axis=axis, nfft=nfft, @@ -663,8 +659,10 @@ def FFT( fftshift_after=fftshift_after, dtype=dtype, ) - elif engine == "scipy": - f = _FFT_scipy( + 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: + logging.warning(pyfftw_message) + f = _FFT_numpy( dims, axis=axis, nfft=nfft, @@ -675,8 +673,8 @@ def FFT( fftshift_after=fftshift_after, dtype=dtype, ) - elif engine == "mkl_fft": - f = _FFT_mklfft( + elif engine == "scipy": + f = _FFT_scipy( dims, axis=axis, nfft=nfft, diff --git a/pylops/signalprocessing/fft2d.py b/pylops/signalprocessing/fft2d.py index 4641601e..b7163e2e 100644 --- a/pylops/signalprocessing/fft2d.py +++ b/pylops/signalprocessing/fft2d.py @@ -1,11 +1,9 @@ __all__ = ["FFT2D"] -import torch import logging import warnings from typing import Dict, Optional, Sequence, Union -from pylops.signalprocessing.fftnd import _FFTND_mklfft import numpy as np import scipy.fft @@ -13,6 +11,12 @@ from pylops.signalprocessing._baseffts import _BaseFFTND, _FFTNorms from pylops.utils.decorators import reshaped from pylops.utils.typing import DTypeLike, InputDimsLike +from pylops.utils import deps +mkl_fft_message = deps.mkl_fft_import("the mkl fft module") + +if mkl_fft_message is None: + from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_ifftshift + from pylops.signalprocessing.fftnd import _FFTND_mklfft logging.basicConfig(format="%(levelname)s: %(message)s", level=logging.WARNING) @@ -266,8 +270,6 @@ def __init__( @reshaped def _matvec(self, x): - from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_ifftshift - if self.ifftshift_before.any(): x = mkl_ifftshift(x, axes=self.axes[self.ifftshift_before]) if not self.clinear: @@ -288,8 +290,6 @@ def _matvec(self, x): @reshaped def _rmatvec(self, x): - from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_ifftshift - if self.fftshift_after.any(): x = mkl_ifftshift(x, axes=self.axes[self.fftshift_after]) if self.real: @@ -500,8 +500,8 @@ def FFT2D( signals. """ - if engine == "numpy": - f = _FFT2D_numpy( + if engine == "mkl_fft" and mkl_fft_message is None: + f = _FFT2D_mklfft( dims=dims, axes=axes, nffts=nffts, @@ -512,8 +512,8 @@ def FFT2D( fftshift_after=fftshift_after, dtype=dtype, ) - elif engine == "scipy": - f = _FFT2D_scipy( + elif engine == "numpy" or (engine == "mkl_fft" and mkl_fft_message is not None): + f = _FFT2D_numpy( dims=dims, axes=axes, nffts=nffts, @@ -524,8 +524,8 @@ def FFT2D( fftshift_after=fftshift_after, dtype=dtype, ) - elif engine == "mkl_fft": - f = _FFT2D_mklfft( + elif engine == "scipy": + f = _FFT2D_scipy( dims=dims, axes=axes, nffts=nffts, diff --git a/pylops/signalprocessing/fftnd.py b/pylops/signalprocessing/fftnd.py index e407f7e6..c5a77d65 100644 --- a/pylops/signalprocessing/fftnd.py +++ b/pylops/signalprocessing/fftnd.py @@ -1,6 +1,5 @@ __all__ = ["FFTND"] -import torch import logging import warnings from typing import Optional, Sequence, Union @@ -12,6 +11,13 @@ from pylops.utils.backend import get_sp_fft from pylops.utils.decorators import reshaped from pylops.utils.typing import DTypeLike, InputDimsLike, NDArray +from pylops.utils import deps + +mkl_fft_message = deps.mkl_fft_import("the mkl fft module") + +if mkl_fft_message is None: + from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_ifftshift + from mkl_fft._numpy_fft import asarray, _cook_nd_args, fft, rfft, irfft, ifft logging.basicConfig(format="%(levelname)s: %(message)s", level=logging.WARNING) @@ -236,10 +242,8 @@ def __init__( @reshaped def _matvec(self, x: NDArray) -> NDArray: - from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift - if self.ifftshift_before.any(): - x = mkl_iffshift(x, axes=self.axes[self.ifftshift_before]) + x = mkl_ifftshift(x, axes=self.axes[self.ifftshift_before]) if not self.clinear: x = np.real(x) if self.real: @@ -258,10 +262,8 @@ def _matvec(self, x: NDArray) -> NDArray: @reshaped def _rmatvec(self, x: NDArray) -> NDArray: - from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_iffshift - if self.fftshift_after.any(): - x = mkl_iffshift(x, axes=self.axes[self.fftshift_after]) + x = mkl_ifftshift(x, axes=self.axes[self.fftshift_after]) if self.real: # Apply scaling to obtain a correct adjoint for this operator x = x.copy() @@ -286,8 +288,6 @@ def _rmatvec(self, x: NDArray) -> NDArray: @staticmethod def rfftn(a, s=None, axes=None, norm=None): - from mkl_fft._numpy_fft import asarray, _cook_nd_args, fft, rfft - a = asarray(a) s, axes = _cook_nd_args(a, s, axes) a = rfft(a, s[-1], axes[-1], norm) @@ -297,8 +297,6 @@ def rfftn(a, s=None, axes=None, norm=None): @staticmethod def irfftn(a, s=None, axes=None, norm=None): - from mkl_fft._numpy_fft import ifft, asarray, _cook_nd_args, irfft - a = asarray(a) s, axes = _cook_nd_args(a, s, axes, invreal=1) for ii in range(len(axes) - 1): @@ -308,8 +306,6 @@ def irfftn(a, s=None, axes=None, norm=None): @staticmethod def fftn(a, s=None, axes=None, norm=None): - from mkl_fft._numpy_fft import asarray, _cook_nd_args, fft - a = asarray(a) s, axes = _cook_nd_args(a, s, axes) itl = list(range(len(axes))) @@ -320,8 +316,6 @@ def fftn(a, s=None, axes=None, norm=None): @staticmethod def ifftn(a, s=None, axes=None, norm=None): - from mkl_fft._numpy_fft import ifft, asarray, _cook_nd_args - a = asarray(a) s, axes = _cook_nd_args(a, s, axes) itl = list(range(len(axes))) @@ -528,8 +522,8 @@ def FFTND( for real input signals. """ - if engine == "numpy": - f = _FFTND_numpy( + if engine == "mkl_fft" and mkl_fft_message is None: + f = _FFTND_mklfft( dims=dims, axes=axes, nffts=nffts, @@ -540,8 +534,8 @@ def FFTND( fftshift_after=fftshift_after, dtype=dtype, ) - elif engine == "scipy": - f = _FFTND_scipy( + elif engine == "numpy" or (engine == "mkl_fft" and mkl_fft_message is not None): + f = _FFTND_numpy( dims=dims, axes=axes, nffts=nffts, @@ -552,8 +546,8 @@ def FFTND( fftshift_after=fftshift_after, dtype=dtype, ) - elif engine == "mkl_fft": - f = _FFTND_mklfft( + elif engine == "scipy": + f = _FFTND_scipy( dims=dims, axes=axes, nffts=nffts, diff --git a/pylops/utils/deps.py b/pylops/utils/deps.py index 7fad2838..65401d50 100644 --- a/pylops/utils/deps.py +++ b/pylops/utils/deps.py @@ -29,6 +29,7 @@ spgl1_enabled = util.find_spec("spgl1") is not None sympy_enabled = util.find_spec("sympy") is not None torch_enabled = util.find_spec("torch") is not None +mkl_fft_enabled = util.find_spec("mkl-fft") is not None # error message at import of available package @@ -87,6 +88,24 @@ def pyfftw_import(message): return pyfftw_message +def mkl_fft_import(message): + if pyfftw_enabled: + try: + import mkl_fft # noqa: F401 + mkl_fft_message = None + except Exception as e: + mkl_fft_message = f"Failed to import pyfftw (error:{e}), use numpy." + else: + mkl_fft_message = ( + "Pyfftw not available, reverting to numpy. " + "In order to be able to use " + f"{message} run " + f'"pip install pyFFTW" or ' + f'"conda install -c conda-forge pyfftw".' + ) + return mkl_fft_message + + def pywt_import(message): if pywt_enabled: try: From 37a48becfb257e15711de9f02b1b2474753a968c Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Thu, 16 Feb 2023 20:21:27 +0530 Subject: [PATCH 21/48] Updated FFTS and added noqa --- pylops/signalprocessing/fft.py | 20 ++++++++-------- pylops/signalprocessing/fft2d.py | 13 ++++++----- pylops/signalprocessing/fftnd.py | 40 ++++++++++++++++---------------- pylops/utils/deps.py | 5 ++-- 4 files changed, 40 insertions(+), 38 deletions(-) diff --git a/pylops/signalprocessing/fft.py b/pylops/signalprocessing/fft.py index 5210c856..4be01be3 100644 --- a/pylops/signalprocessing/fft.py +++ b/pylops/signalprocessing/fft.py @@ -21,8 +21,8 @@ import pyfftw if mkl_fft_message is None: - import mkl_fft._numpy_fft as pymkl_fft - from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_ifftshift + from mkl_fft import _numpy_fft + from mkl_fft import _scipy_fft_backend logging.basicConfig(format="%(levelname)s: %(message)s", level=logging.WARNING) @@ -408,36 +408,36 @@ def __init__( @reshaped def _matvec(self, x: NDArray) -> NDArray: if self.ifftshift_before: - x = mkl_ifftshift(x, axes=self.axis) + x = _scipy_fft_backend.ifftshift(x, axes=self.axis) if not self.clinear: x = np.real(x) if self.real: - y = pymkl_fft.rfft(x, n=self.nfft, axis=self.axis, **self._norm_kwargs) + y = _numpy_fft.rfft(x, n=self.nfft, axis=self.axis, **self._norm_kwargs) # Apply scaling to obtain a correct adjoint for this operator 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 = pymkl_fft.fft(x, n=self.nfft, axis=self.axis, **self._norm_kwargs) + y = _numpy_fft.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_fftshift(y, axes=self.axis) + y = _scipy_fft_backend.fftshift(y, axes=self.axis) return y @reshaped def _rmatvec(self, x: NDArray) -> NDArray: if self.fftshift_after: - x = mkl_ifftshift(x, axes=self.axis) + x = _scipy_fft_backend.ifftshift(x, axes=self.axis) if self.real: # Apply scaling to obtain a correct adjoint for this operator 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 = pymkl_fft.irfft(x, n=self.nfft, axis=self.axis, **self._norm_kwargs) + y = _numpy_fft.irfft(x, n=self.nfft, axis=self.axis, **self._norm_kwargs) else: - y = pymkl_fft.ifft(x, n=self.nfft, axis=self.axis, **self._norm_kwargs) + y = _numpy_fft.ifft(x, n=self.nfft, axis=self.axis, **self._norm_kwargs) if self.norm is _FFTNorms.NONE: y *= self._scale @@ -449,7 +449,7 @@ def _rmatvec(self, x: NDArray) -> NDArray: if not self.clinear: y = np.real(y) if self.ifftshift_before: - y = mkl_fftshift(y, axes=self.axis) + y = _scipy_fft_backend.fftshift(y, axes=self.axis) return y def __truediv__(self, y): diff --git a/pylops/signalprocessing/fft2d.py b/pylops/signalprocessing/fft2d.py index b7163e2e..788d5f40 100644 --- a/pylops/signalprocessing/fft2d.py +++ b/pylops/signalprocessing/fft2d.py @@ -12,11 +12,12 @@ from pylops.utils.decorators import reshaped from pylops.utils.typing import DTypeLike, InputDimsLike from pylops.utils import deps +from pylops.signalprocessing.fftnd import _FFTND_mklfft + mkl_fft_message = deps.mkl_fft_import("the mkl fft module") if mkl_fft_message is None: - from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_ifftshift - from pylops.signalprocessing.fftnd import _FFTND_mklfft + from mkl_fft import _scipy_fft_backend logging.basicConfig(format="%(levelname)s: %(message)s", level=logging.WARNING) @@ -271,7 +272,7 @@ def __init__( @reshaped def _matvec(self, x): if self.ifftshift_before.any(): - x = mkl_ifftshift(x, axes=self.axes[self.ifftshift_before]) + x = _scipy_fft_backend.ifftshift(x, axes=self.axes[self.ifftshift_before]) if not self.clinear: x = np.real(x) if self.real: @@ -285,13 +286,13 @@ def _matvec(self, x): if self.norm is _FFTNorms.ONE_OVER_N: y *= self._scale if self.fftshift_after.any(): - y = mkl_fftshift(y, axes=self.axes[self.fftshift_after]) + y = _scipy_fft_backend.fftshift(y, axes=self.axes[self.fftshift_after]) return y @reshaped def _rmatvec(self, x): if self.fftshift_after.any(): - x = mkl_ifftshift(x, axes=self.axes[self.fftshift_after]) + x = _scipy_fft_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() @@ -308,7 +309,7 @@ def _rmatvec(self, x): if not self.clinear: y = np.real(y) if self.ifftshift_before.any(): - y = mkl_fftshift(y, axes=self.axes[self.ifftshift_before]) + y = _scipy_fft_backend.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 c5a77d65..a016a2c0 100644 --- a/pylops/signalprocessing/fftnd.py +++ b/pylops/signalprocessing/fftnd.py @@ -16,8 +16,8 @@ mkl_fft_message = deps.mkl_fft_import("the mkl fft module") if mkl_fft_message is None: - from mkl_fft._scipy_fft_backend import fftshift as mkl_fftshift, ifftshift as mkl_ifftshift - from mkl_fft._numpy_fft import asarray, _cook_nd_args, fft, rfft, irfft, ifft + from mkl_fft import _numpy_fft + from mkl_fft import _scipy_fft_backend logging.basicConfig(format="%(levelname)s: %(message)s", level=logging.WARNING) @@ -243,7 +243,7 @@ def __init__( @reshaped def _matvec(self, x: NDArray) -> NDArray: if self.ifftshift_before.any(): - x = mkl_ifftshift(x, axes=self.axes[self.ifftshift_before]) + x = _scipy_fft_backend.ifftshift(x, axes=self.axes[self.ifftshift_before]) if not self.clinear: x = np.real(x) if self.real: @@ -257,13 +257,13 @@ def _matvec(self, x: NDArray) -> NDArray: if self.norm is _FFTNorms.ONE_OVER_N: y *= self._scale if self.fftshift_after.any(): - y = mkl_fftshift(y, axes=self.axes[self.fftshift_after]) + y = _scipy_fft_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_ifftshift(x, axes=self.axes[self.fftshift_after]) + x = _scipy_fft_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() @@ -283,45 +283,45 @@ def _rmatvec(self, x: NDArray) -> NDArray: if not self.clinear: y = np.real(y) if self.ifftshift_before.any(): - y = mkl_fftshift(y, axes=self.axes[self.ifftshift_before]) + y = _scipy_fft_backend.fftshift(y, axes=self.axes[self.ifftshift_before]) return y @staticmethod def rfftn(a, s=None, axes=None, norm=None): - a = asarray(a) - s, axes = _cook_nd_args(a, s, axes) - a = rfft(a, s[-1], axes[-1], norm) + a = np.asarray(a) + s, axes = _numpy_fft._cook_nd_args(a, s, axes) + a = _numpy_fft.rfft(a, s[-1], axes[-1], norm) for ii in range(len(axes) - 1): - a = fft(a, s[ii], axes[ii], norm) + a = _numpy_fft.fft(a, s[ii], axes[ii], norm) return a @staticmethod def irfftn(a, s=None, axes=None, norm=None): - a = asarray(a) - s, axes = _cook_nd_args(a, s, axes, invreal=1) + a = np.asarray(a) + s, axes = _numpy_fft._cook_nd_args(a, s, axes, invreal=1) for ii in range(len(axes) - 1): - a = ifft(a, s[ii], axes[ii], norm) - a = irfft(a, s[-1], axes[-1], norm) + a = _numpy_fft.ifft(a, s[ii], axes[ii], norm) + a = _numpy_fft.irfft(a, s[-1], axes[-1], norm) return a @staticmethod def fftn(a, s=None, axes=None, norm=None): - a = asarray(a) - s, axes = _cook_nd_args(a, s, axes) + a = np.asarray(a) + s, axes = _numpy_fft._cook_nd_args(a, s, axes) itl = list(range(len(axes))) itl.reverse() for ii in itl: - a = fft(a, n=s[ii], axis=axes[ii], norm=norm) + a = _numpy_fft.fft(a, n=s[ii], axis=axes[ii], norm=norm) return a @staticmethod def ifftn(a, s=None, axes=None, norm=None): - a = asarray(a) - s, axes = _cook_nd_args(a, s, axes) + a = np.asarray(a) + s, axes = _numpy_fft._cook_nd_args(a, s, axes) itl = list(range(len(axes))) itl.reverse() for ii in itl: - a = ifft(a, n=s[ii], axis=axes[ii], norm=norm) + a = _numpy_fft.ifft(a, n=s[ii], axis=axes[ii], norm=norm) return a def __truediv__(self, y: npt.ArrayLike) -> npt.ArrayLike: diff --git a/pylops/utils/deps.py b/pylops/utils/deps.py index 65401d50..fca95855 100644 --- a/pylops/utils/deps.py +++ b/pylops/utils/deps.py @@ -89,9 +89,10 @@ def pyfftw_import(message): def mkl_fft_import(message): - if pyfftw_enabled: + if mkl_fft_enabled: try: - import mkl_fft # noqa: F401 + from mkl_fft import _scipy_fft_backend # noqa: F401 + from mkl_fft import _numpy_fft # noqa: F401 mkl_fft_message = None except Exception as e: mkl_fft_message = f"Failed to import pyfftw (error:{e}), use numpy." From f838c1e0748a094cecf8a99a6e12348fc3abde1a Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Thu, 16 Feb 2023 21:14:40 +0530 Subject: [PATCH 22/48] Updated deps --- pylops/utils/deps.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pylops/utils/deps.py b/pylops/utils/deps.py index fca95855..17543226 100644 --- a/pylops/utils/deps.py +++ b/pylops/utils/deps.py @@ -9,6 +9,7 @@ "spgl1_enabled", "sympy_enabled", "torch_enabled", + "mkl_fft_enabled" ] import os @@ -29,7 +30,7 @@ spgl1_enabled = util.find_spec("spgl1") is not None sympy_enabled = util.find_spec("sympy") is not None torch_enabled = util.find_spec("torch") is not None -mkl_fft_enabled = util.find_spec("mkl-fft") is not None +mkl_fft_enabled = util.find_spec("mkl_fft") is not None # error message at import of available package From b4656a97ea0eff16ea4b9a3a7c1ecea770b05059 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Fri, 17 Feb 2023 09:57:14 +0530 Subject: [PATCH 23/48] Updated deps --- pylops/utils/deps.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pylops/utils/deps.py b/pylops/utils/deps.py index 17543226..ebb5c330 100644 --- a/pylops/utils/deps.py +++ b/pylops/utils/deps.py @@ -99,11 +99,11 @@ def mkl_fft_import(message): mkl_fft_message = f"Failed to import pyfftw (error:{e}), use numpy." else: mkl_fft_message = ( - "Pyfftw not available, reverting to numpy. " + "mkl_fft not available, reverting to numpy. " "In order to be able to use " f"{message} run " - f'"pip install pyFFTW" or ' - f'"conda install -c conda-forge pyfftw".' + f'"pip install mkl_fft" or ' + f'"conda install -c conda-forge mkl_fft".' ) return mkl_fft_message From 8da2c3ed06e7a8925b594a6258bc4cd7877eb5bd Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Fri, 17 Feb 2023 12:28:46 +0530 Subject: [PATCH 24/48] Updated azure-pipelines.yml --- azure-pipelines.yml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/azure-pipelines.yml b/azure-pipelines.yml index d778e70a..0e4d89d5 100644 --- a/azure-pipelines.yml +++ b/azure-pipelines.yml @@ -51,6 +51,9 @@ jobs: variables: NUMBA_NUM_THREADS: 1 + MKL_THREADING_LAYER: 'sequential' + MKL_ENABLE_INSTRUCTIONS: 'SSE4_2' + steps: - task: UsePythonVersion@0 @@ -80,6 +83,8 @@ jobs: variables: NUMBA_NUM_THREADS: 1 + MKL_THREADING_LAYER: 'sequential' + MKL_ENABLE_INSTRUCTIONS: 'SSE4_2' steps: - task: UsePythonVersion@0 From fafe25e3f6f11f94879ddf3817afaa177675dbb4 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Fri, 17 Feb 2023 22:25:23 +0530 Subject: [PATCH 25/48] Added mkl_fft_message --- pylops/utils/deps.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pylops/utils/deps.py b/pylops/utils/deps.py index ebb5c330..12a9cc23 100644 --- a/pylops/utils/deps.py +++ b/pylops/utils/deps.py @@ -24,13 +24,13 @@ ) devito_enabled = util.find_spec("devito") is not None numba_enabled = util.find_spec("numba") is not None +mkl_fft_enabled = util.find_spec("mkl_fft") is not None pyfftw_enabled = util.find_spec("pyfftw") is not None pywt_enabled = util.find_spec("pywt") is not None skfmm_enabled = util.find_spec("skfmm") is not None spgl1_enabled = util.find_spec("spgl1") is not None sympy_enabled = util.find_spec("sympy") is not None torch_enabled = util.find_spec("torch") is not None -mkl_fft_enabled = util.find_spec("mkl_fft") is not None # error message at import of available package @@ -96,7 +96,7 @@ def mkl_fft_import(message): from mkl_fft import _numpy_fft # noqa: F401 mkl_fft_message = None except Exception as e: - mkl_fft_message = f"Failed to import pyfftw (error:{e}), use numpy." + mkl_fft_message = f"Failed to import mkl_fft (error:{e}), use numpy." else: mkl_fft_message = ( "mkl_fft not available, reverting to numpy. " From a03d1e4ce03282af619a8c3d4f65f4b55af7a89d Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Fri, 17 Feb 2023 22:40:30 +0530 Subject: [PATCH 26/48] Updated Torch commands for mkl_fft --- pylops/torchoperator.py | 4 +++- pylops/utils/typing.py | 11 ----------- 2 files changed, 3 insertions(+), 12 deletions(-) diff --git a/pylops/torchoperator.py b/pylops/torchoperator.py index b5a968a5..dfc3f4da 100644 --- a/pylops/torchoperator.py +++ b/pylops/torchoperator.py @@ -10,6 +10,8 @@ from pylops.utils import deps if deps.torch_enabled: + import torch + TensorTypeLike = torch.Tensor from pylops._torchoperator import _TorchOperator else: torch_message = ( @@ -17,7 +19,7 @@ 'the twoway module run "pip install torch" or' '"conda install -c pytorch torch".' ) -from pylops.utils.typing import TensorTypeLike + TensorTypeLike = None class TorchOperator(LinearOperator): diff --git a/pylops/utils/typing.py b/pylops/utils/typing.py index 191efbfa..8945cc5a 100644 --- a/pylops/utils/typing.py +++ b/pylops/utils/typing.py @@ -5,7 +5,6 @@ "SamplingLike", "ShapeLike", "DTypeLike", - "TensorTypeLike", ] from typing import Sequence, Tuple, Union @@ -13,11 +12,6 @@ import numpy as np import numpy.typing as npt -from pylops.utils.deps import torch_enabled - -if torch_enabled: - import torch - IntNDArray = npt.NDArray[np.int_] NDArray = npt.NDArray @@ -25,8 +19,3 @@ SamplingLike = Union[Sequence[float], NDArray] ShapeLike = Tuple[int, ...] DTypeLike = npt.DTypeLike - -if torch_enabled: - TensorTypeLike = torch.Tensor -else: - TensorTypeLike = None From bee7e9feba253bbec3c6ac1e78cec2607aaa51b6 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Fri, 17 Feb 2023 22:40:43 +0530 Subject: [PATCH 27/48] Rearranged imports --- pylops/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pylops/__init__.py b/pylops/__init__.py index 55d4ce3d..a4623056 100755 --- a/pylops/__init__.py +++ b/pylops/__init__.py @@ -47,7 +47,6 @@ from .config import * from .linearoperator import * -from .torchoperator import * from .basicoperators import * from . import ( avo, @@ -57,6 +56,7 @@ utils, waveeqprocessing, ) +from .torchoperator import * from .avo.poststack import * from .avo.prestack import * from .optimization.basic import * From 899ddd41c806fa366c7c4deb81a8159e7a9dc775 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Fri, 17 Feb 2023 22:42:02 +0530 Subject: [PATCH 28/48] Reverted azure-pipelines.yml --- azure-pipelines.yml | 5 ----- 1 file changed, 5 deletions(-) diff --git a/azure-pipelines.yml b/azure-pipelines.yml index 0e4d89d5..d778e70a 100644 --- a/azure-pipelines.yml +++ b/azure-pipelines.yml @@ -51,9 +51,6 @@ jobs: variables: NUMBA_NUM_THREADS: 1 - MKL_THREADING_LAYER: 'sequential' - MKL_ENABLE_INSTRUCTIONS: 'SSE4_2' - steps: - task: UsePythonVersion@0 @@ -83,8 +80,6 @@ jobs: variables: NUMBA_NUM_THREADS: 1 - MKL_THREADING_LAYER: 'sequential' - MKL_ENABLE_INSTRUCTIONS: 'SSE4_2' steps: - task: UsePythonVersion@0 From 04c966ac7959f5f83a88ea25b4c2dab3b3a3b280 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Fri, 17 Feb 2023 22:59:43 +0530 Subject: [PATCH 29/48] Test Updated plot_fft.py --- examples/plot_fft.py | 130 +++++++++++++++++++++---------------------- 1 file changed, 65 insertions(+), 65 deletions(-) diff --git a/examples/plot_fft.py b/examples/plot_fft.py index 6bb6b01a..9a1c2d97 100644 --- a/examples/plot_fft.py +++ b/examples/plot_fft.py @@ -120,37 +120,37 @@ fig.tight_layout() ############################################################################### -# We can also apply the one dimensional FFT to a two-dimensional -# signal (along one of the first axis) using intel mkl_fft -dt = 0.005 -nt, nx = 100, 20 -t = np.arange(nt) * dt -f0 = 10 -nfft = 2**10 -d = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1) - -FFTop = pylops.signalprocessing.FFT(dims=(nt, nx), axis=0, nfft=nfft, sampling=dt, engine="mkl_fft") -D = FFTop * d.ravel() - -# Adjoint = inverse for FFT -dinv = FFTop.H * D -dinv = FFTop / D -dinv = np.real(dinv).reshape(nt, nx) - -fig, axs = plt.subplots(2, 2, figsize=(10, 6)) -axs[0][0].imshow(d, vmin=-20, vmax=20, cmap="bwr") -axs[0][0].set_title("Signal") -axs[0][0].axis("tight") -axs[0][1].imshow(np.abs(D.reshape(nfft, nx)[:200, :]), cmap="bwr") -axs[0][1].set_title("Fourier Transform using mkl_fft") -axs[0][1].axis("tight") -axs[1][0].imshow(dinv, vmin=-20, vmax=20, cmap="bwr") -axs[1][0].set_title("Inverted") -axs[1][0].axis("tight") -axs[1][1].imshow(d - dinv, vmin=-20, vmax=20, cmap="bwr") -axs[1][1].set_title("Error") -axs[1][1].axis("tight") -fig.tight_layout() +# # We can also apply the one dimensional FFT to a two-dimensional +# # signal (along one of the first axis) using intel mkl_fft +# dt = 0.005 +# nt, nx = 100, 20 +# t = np.arange(nt) * dt +# f0 = 10 +# nfft = 2**10 +# d = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1) +# +# FFTop = pylops.signalprocessing.FFT(dims=(nt, nx), axis=0, nfft=nfft, sampling=dt, engine="mkl_fft") +# D = FFTop * d.ravel() +# +# # Adjoint = inverse for FFT +# dinv = FFTop.H * D +# dinv = FFTop / D +# dinv = np.real(dinv).reshape(nt, nx) +# +# fig, axs = plt.subplots(2, 2, figsize=(10, 6)) +# axs[0][0].imshow(d, vmin=-20, vmax=20, cmap="bwr") +# axs[0][0].set_title("Signal") +# axs[0][0].axis("tight") +# axs[0][1].imshow(np.abs(D.reshape(nfft, nx)[:200, :]), cmap="bwr") +# axs[0][1].set_title("Fourier Transform using mkl_fft") +# axs[0][1].axis("tight") +# axs[1][0].imshow(dinv, vmin=-20, vmax=20, cmap="bwr") +# axs[1][0].set_title("Inverted") +# axs[1][0].axis("tight") +# axs[1][1].imshow(d - dinv, vmin=-20, vmax=20, cmap="bwr") +# axs[1][1].set_title("Error") +# axs[1][1].axis("tight") +# fig.tight_layout() ############################################################################### # We can also apply the two dimensional FFT to to a two-dimensional signal @@ -189,40 +189,40 @@ fig.tight_layout() ############################################################################### -# We can also apply the two-dimensional FFT to a two-dimensional signal using intel mkl_fft -dt, dx = 0.005, 5 -nt, nx = 100, 201 -t = np.arange(nt) * dt -x = np.arange(nx) * dx -f0 = 10 -nfft = 2**10 -d = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1) - -FFTop = pylops.signalprocessing.FFT2D( - dims=(nt, nx), nffts=(nfft, nfft), sampling=(dt, dx), engine='mkl_fft' -) -D = FFTop * d.ravel() - -dinv = FFTop.H * D -dinv = FFTop / D -dinv = np.real(dinv).reshape(nt, nx) - -fig, axs = plt.subplots(2, 2, figsize=(10, 6)) -axs[0][0].imshow(d, vmin=-100, vmax=100, cmap="bwr") -axs[0][0].set_title("Signal") -axs[0][0].axis("tight") -axs[0][1].imshow( - np.abs(np.fft.fftshift(D.reshape(nfft, nfft), axes=1)[:200, :]), cmap="bwr" -) -axs[0][1].set_title("Fourier Transform 2D mkl_fft") -axs[0][1].axis("tight") -axs[1][0].imshow(dinv, vmin=-100, vmax=100, cmap="bwr") -axs[1][0].set_title("Inverted") -axs[1][0].axis("tight") -axs[1][1].imshow(d - dinv, vmin=-100, vmax=100, cmap="bwr") -axs[1][1].set_title("Error") -axs[1][1].axis("tight") -fig.tight_layout() +# # We can also apply the two-dimensional FFT to a two-dimensional signal using intel mkl_fft +# dt, dx = 0.005, 5 +# nt, nx = 100, 201 +# t = np.arange(nt) * dt +# x = np.arange(nx) * dx +# f0 = 10 +# nfft = 2**10 +# d = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1) +# +# FFTop = pylops.signalprocessing.FFT2D( +# dims=(nt, nx), nffts=(nfft, nfft), sampling=(dt, dx), engine='mkl_fft' +# ) +# D = FFTop * d.ravel() +# +# dinv = FFTop.H * D +# dinv = FFTop / D +# dinv = np.real(dinv).reshape(nt, nx) +# +# fig, axs = plt.subplots(2, 2, figsize=(10, 6)) +# axs[0][0].imshow(d, vmin=-100, vmax=100, cmap="bwr") +# axs[0][0].set_title("Signal") +# axs[0][0].axis("tight") +# axs[0][1].imshow( +# np.abs(np.fft.fftshift(D.reshape(nfft, nfft), axes=1)[:200, :]), cmap="bwr" +# ) +# axs[0][1].set_title("Fourier Transform 2D mkl_fft") +# axs[0][1].axis("tight") +# axs[1][0].imshow(dinv, vmin=-100, vmax=100, cmap="bwr") +# axs[1][0].set_title("Inverted") +# axs[1][0].axis("tight") +# axs[1][1].imshow(d - dinv, vmin=-100, vmax=100, cmap="bwr") +# axs[1][1].set_title("Error") +# axs[1][1].axis("tight") +# fig.tight_layout() ############################################################################### From abf0253f012de503d6fb2c6b004e862afee3097e Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Fri, 17 Feb 2023 23:05:21 +0530 Subject: [PATCH 30/48] Test Updated plot_fft.py --- examples/plot_fft.py | 76 ++++++++++++++++++++++---------------------- 1 file changed, 38 insertions(+), 38 deletions(-) diff --git a/examples/plot_fft.py b/examples/plot_fft.py index 9a1c2d97..2addd4df 100644 --- a/examples/plot_fft.py +++ b/examples/plot_fft.py @@ -268,41 +268,41 @@ ############################################################################### # Finally can apply the three-dimensional FFT to a three-dimensional signal using mkl_fft -dt, dx, dy = 0.005, 5, 3 -nt, nx, ny = 30, 21, 11 -t = np.arange(nt) * dt -x = np.arange(nx) * dx -y = np.arange(nx) * dy -f0 = 10 -nfft = 2**6 -nfftk = 2**5 - -d = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1) -d = np.tile(d[:, :, np.newaxis], [1, 1, ny]) - -FFTop = pylops.signalprocessing.FFTND( - dims=(nt, nx, ny), nffts=(nfft, nfftk, nfftk), sampling=(dt, dx, dy), engine="mkl_fft" -) -D = FFTop * d.ravel() - -dinv = FFTop.H * D -dinv = FFTop / D -dinv = np.real(dinv).reshape(nt, nx, ny) - -fig, axs = plt.subplots(2, 2, figsize=(10, 6)) -axs[0][0].imshow(d[:, :, ny // 2], vmin=-20, vmax=20, cmap="bwr") -axs[0][0].set_title("Signal") -axs[0][0].axis("tight") -axs[0][1].imshow( - np.abs(np.fft.fftshift(D.reshape(nfft, nfftk, nfftk), axes=1)[:20, :, nfftk // 2]), - cmap="bwr", -) -axs[0][1].set_title("Fourier Transform 3D mkl fft") -axs[0][1].axis("tight") -axs[1][0].imshow(dinv[:, :, ny // 2], vmin=-20, vmax=20, cmap="bwr") -axs[1][0].set_title("Inverted") -axs[1][0].axis("tight") -axs[1][1].imshow(d[:, :, ny // 2] - dinv[:, :, ny // 2], vmin=-20, vmax=20, cmap="bwr") -axs[1][1].set_title("Error") -axs[1][1].axis("tight") -fig.tight_layout() +# dt, dx, dy = 0.005, 5, 3 +# nt, nx, ny = 30, 21, 11 +# t = np.arange(nt) * dt +# x = np.arange(nx) * dx +# y = np.arange(nx) * dy +# f0 = 10 +# nfft = 2**6 +# nfftk = 2**5 +# +# d = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1) +# d = np.tile(d[:, :, np.newaxis], [1, 1, ny]) +# +# FFTop = pylops.signalprocessing.FFTND( +# dims=(nt, nx, ny), nffts=(nfft, nfftk, nfftk), sampling=(dt, dx, dy), engine="mkl_fft" +# ) +# D = FFTop * d.ravel() +# +# dinv = FFTop.H * D +# dinv = FFTop / D +# dinv = np.real(dinv).reshape(nt, nx, ny) +# +# fig, axs = plt.subplots(2, 2, figsize=(10, 6)) +# axs[0][0].imshow(d[:, :, ny // 2], vmin=-20, vmax=20, cmap="bwr") +# axs[0][0].set_title("Signal") +# axs[0][0].axis("tight") +# axs[0][1].imshow( +# np.abs(np.fft.fftshift(D.reshape(nfft, nfftk, nfftk), axes=1)[:20, :, nfftk // 2]), +# cmap="bwr", +# ) +# axs[0][1].set_title("Fourier Transform 3D mkl fft") +# axs[0][1].axis("tight") +# axs[1][0].imshow(dinv[:, :, ny // 2], vmin=-20, vmax=20, cmap="bwr") +# axs[1][0].set_title("Inverted") +# axs[1][0].axis("tight") +# axs[1][1].imshow(d[:, :, ny // 2] - dinv[:, :, ny // 2], vmin=-20, vmax=20, cmap="bwr") +# axs[1][1].set_title("Error") +# axs[1][1].axis("tight") +# fig.tight_layout() From 63aa63f3a2b952cefa66a7427aca2356ffc914c5 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Fri, 17 Feb 2023 23:12:50 +0530 Subject: [PATCH 31/48] Test Updated plot_fft.py --- examples/plot_fft.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/plot_fft.py b/examples/plot_fft.py index 2addd4df..5c24c5bf 100644 --- a/examples/plot_fft.py +++ b/examples/plot_fft.py @@ -67,8 +67,6 @@ plt.tight_layout() ############################################################################### -# PyLops implements a third engine (``engine='mkl_fft'``) which uses the -# well-known `IntelPython mkl_fft ` FFTop = pylops.signalprocessing.FFT(dims=nt, nfft=nfft, sampling=dt, engine="mkl_fft") D = FFTop * d @@ -306,3 +304,4 @@ # axs[1][1].set_title("Error") # axs[1][1].axis("tight") # fig.tight_layout() +plt.show() From 9819b218384ae995130915ae1934936411d7e453 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Fri, 17 Feb 2023 23:13:23 +0530 Subject: [PATCH 32/48] Test Updated plot_fft.py --- examples/plot_fft.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/plot_fft.py b/examples/plot_fft.py index 5c24c5bf..872a4660 100644 --- a/examples/plot_fft.py +++ b/examples/plot_fft.py @@ -82,7 +82,7 @@ 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() +fig.tight_layout() ############################################################################### # We can also apply the one dimensional FFT to to a two-dimensional From d319332bc3ef812f2507c1df4ef3d2534b417289 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Fri, 17 Feb 2023 23:55:31 +0530 Subject: [PATCH 33/48] Updated plot_fft.py --- examples/plot_fft.py | 116 ++----------------------------------------- 1 file changed, 3 insertions(+), 113 deletions(-) diff --git a/examples/plot_fft.py b/examples/plot_fft.py index 872a4660..3017fc46 100644 --- a/examples/plot_fft.py +++ b/examples/plot_fft.py @@ -67,6 +67,8 @@ plt.tight_layout() ############################################################################### +# PyLops implements a third engine (``engine='mkl_fft'``) which uses the +# well-known `mkl_fft `_ . FFTop = pylops.signalprocessing.FFT(dims=nt, nfft=nfft, sampling=dt, engine="mkl_fft") D = FFTop * d @@ -82,7 +84,7 @@ 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]) -fig.tight_layout() +plt.tight_layout() ############################################################################### # We can also apply the one dimensional FFT to to a two-dimensional @@ -117,39 +119,6 @@ axs[1][1].axis("tight") fig.tight_layout() -############################################################################### -# # We can also apply the one dimensional FFT to a two-dimensional -# # signal (along one of the first axis) using intel mkl_fft -# dt = 0.005 -# nt, nx = 100, 20 -# t = np.arange(nt) * dt -# f0 = 10 -# nfft = 2**10 -# d = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1) -# -# FFTop = pylops.signalprocessing.FFT(dims=(nt, nx), axis=0, nfft=nfft, sampling=dt, engine="mkl_fft") -# D = FFTop * d.ravel() -# -# # Adjoint = inverse for FFT -# dinv = FFTop.H * D -# dinv = FFTop / D -# dinv = np.real(dinv).reshape(nt, nx) -# -# fig, axs = plt.subplots(2, 2, figsize=(10, 6)) -# axs[0][0].imshow(d, vmin=-20, vmax=20, cmap="bwr") -# axs[0][0].set_title("Signal") -# axs[0][0].axis("tight") -# axs[0][1].imshow(np.abs(D.reshape(nfft, nx)[:200, :]), cmap="bwr") -# axs[0][1].set_title("Fourier Transform using mkl_fft") -# axs[0][1].axis("tight") -# axs[1][0].imshow(dinv, vmin=-20, vmax=20, cmap="bwr") -# axs[1][0].set_title("Inverted") -# axs[1][0].axis("tight") -# axs[1][1].imshow(d - dinv, vmin=-20, vmax=20, cmap="bwr") -# axs[1][1].set_title("Error") -# axs[1][1].axis("tight") -# fig.tight_layout() - ############################################################################### # We can also apply the two dimensional FFT to to a two-dimensional signal dt, dx = 0.005, 5 @@ -186,43 +155,6 @@ axs[1][1].axis("tight") fig.tight_layout() -############################################################################### -# # We can also apply the two-dimensional FFT to a two-dimensional signal using intel mkl_fft -# dt, dx = 0.005, 5 -# nt, nx = 100, 201 -# t = np.arange(nt) * dt -# x = np.arange(nx) * dx -# f0 = 10 -# nfft = 2**10 -# d = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1) -# -# FFTop = pylops.signalprocessing.FFT2D( -# dims=(nt, nx), nffts=(nfft, nfft), sampling=(dt, dx), engine='mkl_fft' -# ) -# D = FFTop * d.ravel() -# -# dinv = FFTop.H * D -# dinv = FFTop / D -# dinv = np.real(dinv).reshape(nt, nx) -# -# fig, axs = plt.subplots(2, 2, figsize=(10, 6)) -# axs[0][0].imshow(d, vmin=-100, vmax=100, cmap="bwr") -# axs[0][0].set_title("Signal") -# axs[0][0].axis("tight") -# axs[0][1].imshow( -# np.abs(np.fft.fftshift(D.reshape(nfft, nfft), axes=1)[:200, :]), cmap="bwr" -# ) -# axs[0][1].set_title("Fourier Transform 2D mkl_fft") -# axs[0][1].axis("tight") -# axs[1][0].imshow(dinv, vmin=-100, vmax=100, cmap="bwr") -# axs[1][0].set_title("Inverted") -# axs[1][0].axis("tight") -# axs[1][1].imshow(d - dinv, vmin=-100, vmax=100, cmap="bwr") -# axs[1][1].set_title("Error") -# axs[1][1].axis("tight") -# fig.tight_layout() - - ############################################################################### # Finally can apply the three dimensional FFT to to a three-dimensional signal dt, dx, dy = 0.005, 5, 3 @@ -263,45 +195,3 @@ axs[1][1].set_title("Error") axs[1][1].axis("tight") fig.tight_layout() - -############################################################################### -# Finally can apply the three-dimensional FFT to a three-dimensional signal using mkl_fft -# dt, dx, dy = 0.005, 5, 3 -# nt, nx, ny = 30, 21, 11 -# t = np.arange(nt) * dt -# x = np.arange(nx) * dx -# y = np.arange(nx) * dy -# f0 = 10 -# nfft = 2**6 -# nfftk = 2**5 -# -# d = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1) -# d = np.tile(d[:, :, np.newaxis], [1, 1, ny]) -# -# FFTop = pylops.signalprocessing.FFTND( -# dims=(nt, nx, ny), nffts=(nfft, nfftk, nfftk), sampling=(dt, dx, dy), engine="mkl_fft" -# ) -# D = FFTop * d.ravel() -# -# dinv = FFTop.H * D -# dinv = FFTop / D -# dinv = np.real(dinv).reshape(nt, nx, ny) -# -# fig, axs = plt.subplots(2, 2, figsize=(10, 6)) -# axs[0][0].imshow(d[:, :, ny // 2], vmin=-20, vmax=20, cmap="bwr") -# axs[0][0].set_title("Signal") -# axs[0][0].axis("tight") -# axs[0][1].imshow( -# np.abs(np.fft.fftshift(D.reshape(nfft, nfftk, nfftk), axes=1)[:20, :, nfftk // 2]), -# cmap="bwr", -# ) -# axs[0][1].set_title("Fourier Transform 3D mkl fft") -# axs[0][1].axis("tight") -# axs[1][0].imshow(dinv[:, :, ny // 2], vmin=-20, vmax=20, cmap="bwr") -# axs[1][0].set_title("Inverted") -# axs[1][0].axis("tight") -# axs[1][1].imshow(d[:, :, ny // 2] - dinv[:, :, ny // 2], vmin=-20, vmax=20, cmap="bwr") -# axs[1][1].set_title("Error") -# axs[1][1].axis("tight") -# fig.tight_layout() -plt.show() From e473d84890bdafd58d64fc39781313942045ddca Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Sat, 18 Feb 2023 01:07:59 +0530 Subject: [PATCH 34/48] Updated plot_fft.py --- examples/plot_fft.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/examples/plot_fft.py b/examples/plot_fft.py index 3017fc46..fa59d2d9 100644 --- a/examples/plot_fft.py +++ b/examples/plot_fft.py @@ -6,10 +6,11 @@ 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 pylops import numpy as np -import pylops +import matplotlib.pyplot as plt + plt.close("all") @@ -195,3 +196,4 @@ axs[1][1].set_title("Error") axs[1][1].axis("tight") fig.tight_layout() +plt.show() From a2743d05797d285ab3377bb1156fed4ec272f2f1 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Sat, 18 Feb 2023 11:12:42 +0530 Subject: [PATCH 35/48] Updated environment-dev.yml and setup.py --- environment-dev.yml | 2 +- setup.py | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/environment-dev.yml b/environment-dev.yml index 99687625..b2a18d23 100755 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -8,7 +8,7 @@ dependencies: - python>=3.6.4 - pip - mkl - - mkl-fft + - mkl_fft - mkl-service - numpy>=1.21.0 - scipy>=1.4.0 diff --git a/setup.py b/setup.py index 8b82afa2..0e0e8eed 100755 --- a/setup.py +++ b/setup.py @@ -36,6 +36,9 @@ def src(pth): extras_require={ "advanced": [ "llvmlite", + "mkl", + "mkl_fft", + "mkl-service" "numba", "pyfftw", "PyWavelets", From ae9d48ac08b9186c46e11adb450df64107e038bd Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Sat, 18 Feb 2023 11:16:33 +0530 Subject: [PATCH 36/48] Updated plot_fft.py --- examples/plot_fft.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/examples/plot_fft.py b/examples/plot_fft.py index fa59d2d9..3017fc46 100644 --- a/examples/plot_fft.py +++ b/examples/plot_fft.py @@ -6,11 +6,10 @@ and :py:class:`pylops.signalprocessing.FFTND` operators to apply the Fourier Transform to the model and the inverse Fourier Transform to the data. """ -import pylops -import numpy as np - import matplotlib.pyplot as plt +import numpy as np +import pylops plt.close("all") @@ -196,4 +195,3 @@ axs[1][1].set_title("Error") axs[1][1].axis("tight") fig.tight_layout() -plt.show() From 258352039566e7d063acef8013b6e4446143a35a Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Sat, 18 Feb 2023 14:20:18 +0530 Subject: [PATCH 37/48] Test for python 3.8 docs --- readthedocs.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/readthedocs.yml b/readthedocs.yml index 8967c8a2..2c040e59 100644 --- a/readthedocs.yml +++ b/readthedocs.yml @@ -5,10 +5,10 @@ version: 2 build: os: ubuntu-20.04 tools: - python: "3.9" + python: "3.8" python: install: - requirements: requirements-doc.txt - method: setuptools - path: . + path: . \ No newline at end of file From 2372c0811dde94f3d6d73819fdc825dcb7ef170a Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Sat, 18 Feb 2023 14:50:43 +0530 Subject: [PATCH 38/48] Updated readthedocs.yml to test --- readthedocs.yml | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/readthedocs.yml b/readthedocs.yml index 2c040e59..0918f89d 100644 --- a/readthedocs.yml +++ b/readthedocs.yml @@ -2,13 +2,13 @@ version: 2 -build: - os: ubuntu-20.04 - tools: - python: "3.8" +sphinx: + configuration: docs/source/conf.py + fail_on_warning: false python: + version: '3.8' install: + - method: pip + path: . - requirements: requirements-doc.txt - - method: setuptools - path: . \ No newline at end of file From 6cdb3e1ef34563d07b3e26db0018f2cef807e59b Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Sat, 18 Feb 2023 16:42:15 +0530 Subject: [PATCH 39/48] Reverting readthedocs.yml --- pylops/avo/__init__.py | 1 + readthedocs.yml | 12 ++++++------ 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/pylops/avo/__init__.py b/pylops/avo/__init__.py index 2e7bd5c8..8c1655ab 100644 --- a/pylops/avo/__init__.py +++ b/pylops/avo/__init__.py @@ -22,6 +22,7 @@ from .poststack import * from .prestack import * +from .avo import * __all__ = [ "AVOLinearModelling", diff --git a/readthedocs.yml b/readthedocs.yml index 0918f89d..8967c8a2 100644 --- a/readthedocs.yml +++ b/readthedocs.yml @@ -2,13 +2,13 @@ version: 2 -sphinx: - configuration: docs/source/conf.py - fail_on_warning: false +build: + os: ubuntu-20.04 + tools: + python: "3.9" python: - version: '3.8' install: - - method: pip - path: . - requirements: requirements-doc.txt + - method: setuptools + path: . From 2f91d801a1e9004ab456155964b8bce28135fa42 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Sat, 18 Feb 2023 16:45:44 +0530 Subject: [PATCH 40/48] Reverting --- examples/plot_fft.py | 1 + pylops/avo/__init__.py | 1 - 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/plot_fft.py b/examples/plot_fft.py index 3017fc46..3bac892e 100644 --- a/examples/plot_fft.py +++ b/examples/plot_fft.py @@ -155,6 +155,7 @@ axs[1][1].axis("tight") fig.tight_layout() + ############################################################################### # Finally can apply the three dimensional FFT to to a three-dimensional signal dt, dx, dy = 0.005, 5, 3 diff --git a/pylops/avo/__init__.py b/pylops/avo/__init__.py index 8c1655ab..2e7bd5c8 100644 --- a/pylops/avo/__init__.py +++ b/pylops/avo/__init__.py @@ -22,7 +22,6 @@ from .poststack import * from .prestack import * -from .avo import * __all__ = [ "AVOLinearModelling", From 492ee68a040caffb308cca0b1ac7a541dd9f44a6 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Sat, 18 Feb 2023 16:49:49 +0530 Subject: [PATCH 41/48] Reverting --- pylops/signalprocessing/fft.py | 4 +--- pylops/signalprocessing/fft2d.py | 4 +--- pylops/signalprocessing/fftnd.py | 6 +----- 3 files changed, 3 insertions(+), 11 deletions(-) diff --git a/pylops/signalprocessing/fft.py b/pylops/signalprocessing/fft.py index 4be01be3..ce037cb9 100644 --- a/pylops/signalprocessing/fft.py +++ b/pylops/signalprocessing/fft.py @@ -396,7 +396,7 @@ def __init__( fftshift_after=fftshift_after, dtype=dtype, ) - self._norm_kwargs = {"norm": None} # equivalent to "backward" in Numpy/Scipy + self._norm_kwargs = {"norm": None} if self.norm is _FFTNorms.ORTHO: self._norm_kwargs["norm"] = "ortho" self._scale = np.sqrt(1 / self.nfft) @@ -413,7 +413,6 @@ def _matvec(self, x: NDArray) -> NDArray: x = np.real(x) if self.real: y = _numpy_fft.rfft(x, n=self.nfft, axis=self.axis, **self._norm_kwargs) - # Apply scaling to obtain a correct adjoint for this operator y = np.swapaxes(y, -1, self.axis) y[..., 1 : 1 + (self.nfft - 1) // 2] *= np.sqrt(2) y = np.swapaxes(y, self.axis, -1) @@ -430,7 +429,6 @@ def _rmatvec(self, x: NDArray) -> NDArray: if self.fftshift_after: x = _scipy_fft_backend.ifftshift(x, axes=self.axis) if self.real: - # Apply scaling to obtain a correct adjoint for this operator x = x.copy() x = np.swapaxes(x, -1, self.axis) x[..., 1 : 1 + (self.nfft - 1) // 2] /= np.sqrt(2) diff --git a/pylops/signalprocessing/fft2d.py b/pylops/signalprocessing/fft2d.py index 788d5f40..30938cda 100644 --- a/pylops/signalprocessing/fft2d.py +++ b/pylops/signalprocessing/fft2d.py @@ -260,7 +260,7 @@ def __init__( self._norm_kwargs: Dict[str, Union[None, str]] = { "norm": None - } # equivalent to "backward" in Numpy/Scipy + } if self.norm is _FFTNorms.ORTHO: self._norm_kwargs["norm"] = "ortho" self._scale = np.sqrt(1 / np.prod(np.sqrt(self.nffts))) @@ -277,7 +277,6 @@ def _matvec(self, x): x = np.real(x) if self.real: y = _FFTND_mklfft.rfftn(x, s=self.nffts, axes=self.axes, **self._norm_kwargs) - # 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) @@ -294,7 +293,6 @@ def _rmatvec(self, x): if self.fftshift_after.any(): x = _scipy_fft_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) diff --git a/pylops/signalprocessing/fftnd.py b/pylops/signalprocessing/fftnd.py index a016a2c0..7162780c 100644 --- a/pylops/signalprocessing/fftnd.py +++ b/pylops/signalprocessing/fftnd.py @@ -159,7 +159,6 @@ def _matvec(self, x: NDArray) -> NDArray: x = np.real(x) if self.real: y = sp_fft.rfftn(x, s=self.nffts, axes=self.axes, **self._norm_kwargs) - # 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) @@ -177,7 +176,6 @@ def _rmatvec(self, x: NDArray) -> NDArray: if self.fftshift_after.any(): x = sp_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() x = np.swapaxes(x, -1, self.axes[-1]) x[..., 1 : 1 + (self.nffts[-1] - 1) // 2] /= np.sqrt(2) @@ -231,7 +229,7 @@ def __init__( dtype=dtype, ) - self._norm_kwargs = {"norm": None} # equivalent to "backward" in Numpy/Scipy + self._norm_kwargs = {"norm": None} if self.norm is _FFTNorms.ORTHO: self._norm_kwargs["norm"] = "ortho" self._scale = np.sqrt(1 / np.prod(self.nffts)) @@ -248,7 +246,6 @@ def _matvec(self, x: NDArray) -> NDArray: x = np.real(x) if self.real: y = self.rfftn(x, s=self.nffts, axes=tuple(self.axes), **self._norm_kwargs) - # 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) @@ -265,7 +262,6 @@ def _rmatvec(self, x: NDArray) -> NDArray: if self.fftshift_after.any(): x = _scipy_fft_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) From b1a9788073b924c0e4d5820e2e8c4f7a15d99fa2 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Sat, 18 Feb 2023 17:01:09 +0530 Subject: [PATCH 42/48] Reverting --- pylops/signalprocessing/fftnd.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pylops/signalprocessing/fftnd.py b/pylops/signalprocessing/fftnd.py index 7162780c..1def1f62 100644 --- a/pylops/signalprocessing/fftnd.py +++ b/pylops/signalprocessing/fftnd.py @@ -159,6 +159,7 @@ def _matvec(self, x: NDArray) -> NDArray: x = np.real(x) if self.real: y = sp_fft.rfftn(x, s=self.nffts, axes=self.axes, **self._norm_kwargs) + # 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) @@ -176,6 +177,7 @@ def _rmatvec(self, x: NDArray) -> NDArray: if self.fftshift_after.any(): x = sp_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() x = np.swapaxes(x, -1, self.axes[-1]) x[..., 1 : 1 + (self.nffts[-1] - 1) // 2] /= np.sqrt(2) From 35fbbcec75192061f50981dc8b02420b77c44f3e Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Sat, 18 Feb 2023 17:06:56 +0530 Subject: [PATCH 43/48] Test plot_fft.py --- examples/plot_fft.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/examples/plot_fft.py b/examples/plot_fft.py index 3bac892e..fc33b40c 100644 --- a/examples/plot_fft.py +++ b/examples/plot_fft.py @@ -73,18 +73,18 @@ D = FFTop * d # Adjoint = 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() +# 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 From de963611a6bb10de4de3b4f23efb7a40dc707248 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Sat, 18 Feb 2023 17:15:06 +0530 Subject: [PATCH 44/48] Updated ffts --- examples/plot_fft.py | 10 +++++----- pylops/__init__.py | 1 + 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/examples/plot_fft.py b/examples/plot_fft.py index fc33b40c..5f168b21 100644 --- a/examples/plot_fft.py +++ b/examples/plot_fft.py @@ -24,7 +24,7 @@ nfft = 2**10 d = np.sin(2 * np.pi * f0 * t) -FFTop = pylops.signalprocessing.FFT(dims=nt, nfft=nfft, sampling=dt, engine="numpy") +FFTop = pylops.FFT(dims=nt, nfft=nfft, sampling=dt, engine="numpy") D = FFTop * d # Adjoint = inverse for FFT @@ -49,7 +49,7 @@ # numpy in many cases but it is not inserted in the mandatory requirements of # PyLops. If interested to use ``FFTW`` backend, read the `fft routines` # section at :ref:`performance`. -FFTop = pylops.signalprocessing.FFT(dims=nt, nfft=nfft, sampling=dt, engine="fftw") +FFTop = pylops.FFT(dims=nt, nfft=nfft, sampling=dt, engine="fftw") D = FFTop * d # Adjoint = inverse for FFT @@ -69,7 +69,7 @@ ############################################################################### # PyLops implements a third engine (``engine='mkl_fft'``) which uses the # well-known `mkl_fft `_ . -FFTop = pylops.signalprocessing.FFT(dims=nt, nfft=nfft, sampling=dt, engine="mkl_fft") +FFTop = pylops.FFT(dims=nt, nfft=nfft, sampling=dt, engine="mkl_fft") D = FFTop * d # Adjoint = inverse for FFT @@ -129,7 +129,7 @@ nfft = 2**10 d = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1) -FFTop = pylops.signalprocessing.FFT2D( +FFTop = pylops.FFT2D( dims=(nt, nx), nffts=(nfft, nfft), sampling=(dt, dx) ) D = FFTop * d.ravel() @@ -170,7 +170,7 @@ d = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1) d = np.tile(d[:, :, np.newaxis], [1, 1, ny]) -FFTop = pylops.signalprocessing.FFTND( +FFTop = pylops.FFTND( dims=(nt, nx, ny), nffts=(nfft, nfftk, nfftk), sampling=(dt, dx, dy) ) D = FFTop * d.ravel() diff --git a/pylops/__init__.py b/pylops/__init__.py index a4623056..770e628f 100755 --- a/pylops/__init__.py +++ b/pylops/__init__.py @@ -57,6 +57,7 @@ waveeqprocessing, ) from .torchoperator import * +from .signalprocessing import * from .avo.poststack import * from .avo.prestack import * from .optimization.basic import * From 591835bb16ac4db8cc57e40c842b3c727dba4977 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Sun, 19 Feb 2023 00:45:18 +0530 Subject: [PATCH 45/48] Test plot_fft.py --- examples/plot_fft.py | 316 +++++++++++++++++++++---------------------- 1 file changed, 158 insertions(+), 158 deletions(-) diff --git a/examples/plot_fft.py b/examples/plot_fft.py index 5f168b21..fe3f180d 100644 --- a/examples/plot_fft.py +++ b/examples/plot_fft.py @@ -6,65 +6,65 @@ 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 matplotlib.pyplot as plt import numpy as np - +# import pylops - -plt.close("all") - -############################################################################### -# 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 -# lenght :math:`nt` and sampling :math:`dt` +# +# plt.close("all") +# +# ############################################################################### +# # 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 +# # lenght :math:`nt` and sampling :math:`dt` dt = 0.005 nt = 100 t = np.arange(nt) * dt f0 = 10 nfft = 2**10 d = np.sin(2 * np.pi * f0 * t) - -FFTop = pylops.FFT(dims=nt, nfft=nfft, sampling=dt, engine="numpy") -D = FFTop * d - -# Adjoint = 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") -axs[1].set_xlim([0, 3 * f0]) -plt.tight_layout() - -############################################################################### -# In this example we used numpy as our engine for the ``fft`` and ``ifft``. -# PyLops implements a second engine (``engine='fftw'``) which uses the -# well-known `FFTW `_ via the python wrapper -# :py:class:`pyfftw.FFTW`. This optimized fft tends to outperform the one from -# numpy in many cases but it is not inserted in the mandatory requirements of -# PyLops. If interested to use ``FFTW`` backend, read the `fft routines` -# section at :ref:`performance`. -FFTop = pylops.FFT(dims=nt, nfft=nfft, sampling=dt, engine="fftw") -D = FFTop * d - -# Adjoint = 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 fftw") -axs[1].set_xlim([0, 3 * f0]) -plt.tight_layout() +# +# FFTop = pylops.FFT(dims=nt, nfft=nfft, sampling=dt, engine="numpy") +# D = FFTop * d +# +# # Adjoint = 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") +# axs[1].set_xlim([0, 3 * f0]) +# plt.tight_layout() +# +# ############################################################################### +# # In this example we used numpy as our engine for the ``fft`` and ``ifft``. +# # PyLops implements a second engine (``engine='fftw'``) which uses the +# # well-known `FFTW `_ via the python wrapper +# # :py:class:`pyfftw.FFTW`. This optimized fft tends to outperform the one from +# # numpy in many cases but it is not inserted in the mandatory requirements of +# # PyLops. If interested to use ``FFTW`` backend, read the `fft routines` +# # section at :ref:`performance`. +# FFTop = pylops.FFT(dims=nt, nfft=nfft, sampling=dt, engine="fftw") +# D = FFTop * d +# +# # Adjoint = 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 fftw") +# axs[1].set_xlim([0, 3 * f0]) +# plt.tight_layout() ############################################################################### # PyLops implements a third engine (``engine='mkl_fft'``) which uses the @@ -89,110 +89,110 @@ ############################################################################### # We can also apply the one dimensional FFT to to a two-dimensional # signal (along one of the first axis) -dt = 0.005 -nt, nx = 100, 20 -t = np.arange(nt) * dt -f0 = 10 -nfft = 2**10 -d = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1) - -FFTop = pylops.signalprocessing.FFT(dims=(nt, nx), axis=0, nfft=nfft, sampling=dt) -D = FFTop * d.ravel() - -# Adjoint = inverse for FFT -dinv = FFTop.H * D -dinv = FFTop / D -dinv = np.real(dinv).reshape(nt, nx) - -fig, axs = plt.subplots(2, 2, figsize=(10, 6)) -axs[0][0].imshow(d, vmin=-20, vmax=20, cmap="bwr") -axs[0][0].set_title("Signal") -axs[0][0].axis("tight") -axs[0][1].imshow(np.abs(D.reshape(nfft, nx)[:200, :]), cmap="bwr") -axs[0][1].set_title("Fourier Transform") -axs[0][1].axis("tight") -axs[1][0].imshow(dinv, vmin=-20, vmax=20, cmap="bwr") -axs[1][0].set_title("Inverted") -axs[1][0].axis("tight") -axs[1][1].imshow(d - dinv, vmin=-20, vmax=20, cmap="bwr") -axs[1][1].set_title("Error") -axs[1][1].axis("tight") -fig.tight_layout() - -############################################################################### -# We can also apply the two dimensional FFT to to a two-dimensional signal -dt, dx = 0.005, 5 -nt, nx = 100, 201 -t = np.arange(nt) * dt -x = np.arange(nx) * dx -f0 = 10 -nfft = 2**10 -d = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1) - -FFTop = pylops.FFT2D( - dims=(nt, nx), nffts=(nfft, nfft), sampling=(dt, dx) -) -D = FFTop * d.ravel() - -dinv = FFTop.H * D -dinv = FFTop / D -dinv = np.real(dinv).reshape(nt, nx) - -fig, axs = plt.subplots(2, 2, figsize=(10, 6)) -axs[0][0].imshow(d, vmin=-100, vmax=100, cmap="bwr") -axs[0][0].set_title("Signal") -axs[0][0].axis("tight") -axs[0][1].imshow( - np.abs(np.fft.fftshift(D.reshape(nfft, nfft), axes=1)[:200, :]), cmap="bwr" -) -axs[0][1].set_title("Fourier Transform") -axs[0][1].axis("tight") -axs[1][0].imshow(dinv, vmin=-100, vmax=100, cmap="bwr") -axs[1][0].set_title("Inverted") -axs[1][0].axis("tight") -axs[1][1].imshow(d - dinv, vmin=-100, vmax=100, cmap="bwr") -axs[1][1].set_title("Error") -axs[1][1].axis("tight") -fig.tight_layout() - - -############################################################################### -# Finally can apply the three dimensional FFT to to a three-dimensional signal -dt, dx, dy = 0.005, 5, 3 -nt, nx, ny = 30, 21, 11 -t = np.arange(nt) * dt -x = np.arange(nx) * dx -y = np.arange(nx) * dy -f0 = 10 -nfft = 2**6 -nfftk = 2**5 - -d = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1) -d = np.tile(d[:, :, np.newaxis], [1, 1, ny]) - -FFTop = pylops.FFTND( - dims=(nt, nx, ny), nffts=(nfft, nfftk, nfftk), sampling=(dt, dx, dy) -) -D = FFTop * d.ravel() - -dinv = FFTop.H * D -dinv = FFTop / D -dinv = np.real(dinv).reshape(nt, nx, ny) - -fig, axs = plt.subplots(2, 2, figsize=(10, 6)) -axs[0][0].imshow(d[:, :, ny // 2], vmin=-20, vmax=20, cmap="bwr") -axs[0][0].set_title("Signal") -axs[0][0].axis("tight") -axs[0][1].imshow( - np.abs(np.fft.fftshift(D.reshape(nfft, nfftk, nfftk), axes=1)[:20, :, nfftk // 2]), - cmap="bwr", -) -axs[0][1].set_title("Fourier Transform") -axs[0][1].axis("tight") -axs[1][0].imshow(dinv[:, :, ny // 2], vmin=-20, vmax=20, cmap="bwr") -axs[1][0].set_title("Inverted") -axs[1][0].axis("tight") -axs[1][1].imshow(d[:, :, ny // 2] - dinv[:, :, ny // 2], vmin=-20, vmax=20, cmap="bwr") -axs[1][1].set_title("Error") -axs[1][1].axis("tight") -fig.tight_layout() +# dt = 0.005 +# nt, nx = 100, 20 +# t = np.arange(nt) * dt +# f0 = 10 +# nfft = 2**10 +# d = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1) +# +# FFTop = pylops.signalprocessing.FFT(dims=(nt, nx), axis=0, nfft=nfft, sampling=dt) +# D = FFTop * d.ravel() +# +# # Adjoint = inverse for FFT +# dinv = FFTop.H * D +# dinv = FFTop / D +# dinv = np.real(dinv).reshape(nt, nx) +# +# fig, axs = plt.subplots(2, 2, figsize=(10, 6)) +# axs[0][0].imshow(d, vmin=-20, vmax=20, cmap="bwr") +# axs[0][0].set_title("Signal") +# axs[0][0].axis("tight") +# axs[0][1].imshow(np.abs(D.reshape(nfft, nx)[:200, :]), cmap="bwr") +# axs[0][1].set_title("Fourier Transform") +# axs[0][1].axis("tight") +# axs[1][0].imshow(dinv, vmin=-20, vmax=20, cmap="bwr") +# axs[1][0].set_title("Inverted") +# axs[1][0].axis("tight") +# axs[1][1].imshow(d - dinv, vmin=-20, vmax=20, cmap="bwr") +# axs[1][1].set_title("Error") +# axs[1][1].axis("tight") +# fig.tight_layout() +# +# ############################################################################### +# # We can also apply the two dimensional FFT to to a two-dimensional signal +# dt, dx = 0.005, 5 +# nt, nx = 100, 201 +# t = np.arange(nt) * dt +# x = np.arange(nx) * dx +# f0 = 10 +# nfft = 2**10 +# d = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1) +# +# FFTop = pylops.FFT2D( +# dims=(nt, nx), nffts=(nfft, nfft), sampling=(dt, dx) +# ) +# D = FFTop * d.ravel() +# +# dinv = FFTop.H * D +# dinv = FFTop / D +# dinv = np.real(dinv).reshape(nt, nx) +# +# fig, axs = plt.subplots(2, 2, figsize=(10, 6)) +# axs[0][0].imshow(d, vmin=-100, vmax=100, cmap="bwr") +# axs[0][0].set_title("Signal") +# axs[0][0].axis("tight") +# axs[0][1].imshow( +# np.abs(np.fft.fftshift(D.reshape(nfft, nfft), axes=1)[:200, :]), cmap="bwr" +# ) +# axs[0][1].set_title("Fourier Transform") +# axs[0][1].axis("tight") +# axs[1][0].imshow(dinv, vmin=-100, vmax=100, cmap="bwr") +# axs[1][0].set_title("Inverted") +# axs[1][0].axis("tight") +# axs[1][1].imshow(d - dinv, vmin=-100, vmax=100, cmap="bwr") +# axs[1][1].set_title("Error") +# axs[1][1].axis("tight") +# fig.tight_layout() +# +# +# ############################################################################### +# # Finally can apply the three dimensional FFT to to a three-dimensional signal +# dt, dx, dy = 0.005, 5, 3 +# nt, nx, ny = 30, 21, 11 +# t = np.arange(nt) * dt +# x = np.arange(nx) * dx +# y = np.arange(nx) * dy +# f0 = 10 +# nfft = 2**6 +# nfftk = 2**5 +# +# d = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1) +# d = np.tile(d[:, :, np.newaxis], [1, 1, ny]) +# +# FFTop = pylops.FFTND( +# dims=(nt, nx, ny), nffts=(nfft, nfftk, nfftk), sampling=(dt, dx, dy) +# ) +# D = FFTop * d.ravel() +# +# dinv = FFTop.H * D +# dinv = FFTop / D +# dinv = np.real(dinv).reshape(nt, nx, ny) +# +# fig, axs = plt.subplots(2, 2, figsize=(10, 6)) +# axs[0][0].imshow(d[:, :, ny // 2], vmin=-20, vmax=20, cmap="bwr") +# axs[0][0].set_title("Signal") +# axs[0][0].axis("tight") +# axs[0][1].imshow( +# np.abs(np.fft.fftshift(D.reshape(nfft, nfftk, nfftk), axes=1)[:20, :, nfftk // 2]), +# cmap="bwr", +# ) +# axs[0][1].set_title("Fourier Transform") +# axs[0][1].axis("tight") +# axs[1][0].imshow(dinv[:, :, ny // 2], vmin=-20, vmax=20, cmap="bwr") +# axs[1][0].set_title("Inverted") +# axs[1][0].axis("tight") +# axs[1][1].imshow(d[:, :, ny // 2] - dinv[:, :, ny // 2], vmin=-20, vmax=20, cmap="bwr") +# axs[1][1].set_title("Error") +# axs[1][1].axis("tight") +# fig.tight_layout() From 0ec7169a5a2b7893f84317545907ce9b691d0509 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Sun, 19 Feb 2023 00:53:28 +0530 Subject: [PATCH 46/48] Test plot_fft.py --- examples/plot_fft.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/plot_fft.py b/examples/plot_fft.py index fe3f180d..7e1bd251 100644 --- a/examples/plot_fft.py +++ b/examples/plot_fft.py @@ -19,10 +19,10 @@ # # lenght :math:`nt` and sampling :math:`dt` dt = 0.005 nt = 100 -t = np.arange(nt) * dt +# t = np.arange(nt) * dt f0 = 10 nfft = 2**10 -d = np.sin(2 * np.pi * f0 * t) +# d = np.sin(2 * np.pi * f0 * t) # # FFTop = pylops.FFT(dims=nt, nfft=nfft, sampling=dt, engine="numpy") # D = FFTop * d @@ -70,7 +70,7 @@ # PyLops implements a third engine (``engine='mkl_fft'``) which uses the # well-known `mkl_fft `_ . FFTop = pylops.FFT(dims=nt, nfft=nfft, sampling=dt, engine="mkl_fft") -D = FFTop * d +# D = FFTop * d # Adjoint = inverse for FFT # dinv = FFTop.H * D From a2a4924b5f789a52228c172cbe8759179fcedcda Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Sun, 19 Feb 2023 01:07:02 +0530 Subject: [PATCH 47/48] Test plot_fft.py --- examples/plot_fft.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/plot_fft.py b/examples/plot_fft.py index 7e1bd251..fe3f180d 100644 --- a/examples/plot_fft.py +++ b/examples/plot_fft.py @@ -19,10 +19,10 @@ # # lenght :math:`nt` and sampling :math:`dt` dt = 0.005 nt = 100 -# t = np.arange(nt) * dt +t = np.arange(nt) * dt f0 = 10 nfft = 2**10 -# d = np.sin(2 * np.pi * f0 * t) +d = np.sin(2 * np.pi * f0 * t) # # FFTop = pylops.FFT(dims=nt, nfft=nfft, sampling=dt, engine="numpy") # D = FFTop * d @@ -70,7 +70,7 @@ # PyLops implements a third engine (``engine='mkl_fft'``) which uses the # well-known `mkl_fft `_ . FFTop = pylops.FFT(dims=nt, nfft=nfft, sampling=dt, engine="mkl_fft") -# D = FFTop * d +D = FFTop * d # Adjoint = inverse for FFT # dinv = FFTop.H * D From afb3885baf7770974cb3ccc679dc3b0116cc4184 Mon Sep 17 00:00:00 2001 From: rohanbabbar04 Date: Sun, 19 Feb 2023 01:14:51 +0530 Subject: [PATCH 48/48] Test plot_fft.py --- examples/plot_fft.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/plot_fft.py b/examples/plot_fft.py index fe3f180d..10f9db47 100644 --- a/examples/plot_fft.py +++ b/examples/plot_fft.py @@ -7,9 +7,10 @@ Transform to the model and the inverse Fourier Transform to the data. """ # import matplotlib.pyplot as plt +import pylops + import numpy as np # -import pylops # # plt.close("all") #