diff --git a/MIGRATION_V1_V2.md b/MIGRATION_V1_V2.md index 04676668..48a6f54f 100644 --- a/MIGRATION_V1_V2.md +++ b/MIGRATION_V1_V2.md @@ -8,6 +8,7 @@ should be used as a checklist when converting a piece of code using PyLops from - Several operators have deprecated `N` as a keyword. To migrate, pass only `dims` if both `N` and `dims` are currently being passed. If only `N` is being passed, ensure it is being passed as a value and not a keyword argument (e.g., change `Flip(N=100)` to `Flip(100)`). +- `dir`, `dirs` and `nodir` have been deprecated in favor of `axis` and `axes`. When previously `nodir` was used, you must now provide the directions along which the operator is applied through `axes`. The default value for `axis` and `axes` are chosen to be -1 and (-2, -1), respectively, whereas the default `dir` and `dirs` was 0 and (0, 1), respectively. Be careful to check all operators where `dir`, `dirs` and `nodir` was not provided explicitly. - `utils.dottest`: Change `tol` into `rtol`. Absolute tolerance is now also supported via the keyword `atol`. When calling it with purely positional arguments, note that after `rtol` comes now first `atol` before `complexflag`. When using `raiseerror=True` it now emits an `AttributeError` instead of a `ValueError`. diff --git a/examples/plot_causalintegration.py b/examples/plot_causalintegration.py index 30524818..c1ae1bac 100755 --- a/examples/plot_causalintegration.py +++ b/examples/plot_causalintegration.py @@ -112,14 +112,14 @@ t = np.arange(nt) * dt + ot x = np.outer(np.sin(t), np.ones(nx)) -Cop = pylops.CausalIntegration((nt, nx), sampling=dt, dir=0, halfcurrent=True) +Cop = pylops.CausalIntegration(dims=(nt, nx), sampling=dt, axis=0, halfcurrent=True) y = Cop * x.ravel() y = y.reshape(nt, nx) yn = y + np.random.normal(0, 4e-1, y.shape) # Numerical derivative -Dop = pylops.FirstDerivative((nt, nx), dir=0, sampling=dt) +Dop = pylops.FirstDerivative(dims=(nt, nx), axis=0, sampling=dt) xder = Dop * yn.ravel() xder = xder.reshape(nt, nx) diff --git a/examples/plot_convolve.py b/examples/plot_convolve.py index fca20ca5..6c73f6a5 100644 --- a/examples/plot_convolve.py +++ b/examples/plot_convolve.py @@ -143,7 +143,7 @@ offset = [1, 2, 1] Cop = pylops.signalprocessing.ConvolveND( - nx * ny * nz, h=h, offset=offset, dims=[ny, nx, nz], dirs=[0, 1, 2], dtype="float32" + dims=(ny, nx, nz), h=h, offset=offset, axes=(0, 1, 2), dtype="float32" ) y = Cop * x.ravel() xinv = lsqr(Cop, y, damp=0, iter_lim=300, show=0)[0] diff --git a/examples/plot_derivative.py b/examples/plot_derivative.py index dac93987..0f36797e 100644 --- a/examples/plot_derivative.py +++ b/examples/plot_derivative.py @@ -83,7 +83,7 @@ A = np.zeros((nx, ny)) A[nx // 2, ny // 2] = 1.0 -D1op = pylops.FirstDerivative((nx, ny), dir=0, dtype="float64") +D1op = pylops.FirstDerivative((nx, ny), axis=0, dtype="float64") B = np.reshape(D1op * A.ravel(), (nx, ny)) fig, axs = plt.subplots(1, 2, figsize=(10, 3)) @@ -106,7 +106,7 @@ A = np.zeros((nx, ny)) A[nx // 2, ny // 2] = 1.0 -D2op = pylops.SecondDerivative((nx, ny), dir=0, dtype="float64") +D2op = pylops.SecondDerivative(dims=(nx, ny), axis=0, dtype="float64") B = np.reshape(D2op * A.ravel(), (nx, ny)) fig, axs = plt.subplots(1, 2, figsize=(10, 3)) @@ -126,8 +126,8 @@ ############################################################################### # We can also apply the second derivative to the second direction of -# our data (``dir=1``) -D2op = pylops.SecondDerivative((nx, ny), dir=1, dtype="float64") +# our data (``axis=1``) +D2op = pylops.SecondDerivative(dims=(nx, ny), axis=1, dtype="float64") B = np.reshape(D2op * np.ndarray.flatten(A), (nx, ny)) fig, axs = plt.subplots(1, 2, figsize=(10, 3)) diff --git a/examples/plot_diagonal.py b/examples/plot_diagonal.py index 31d9da4a..0de27cb5 100644 --- a/examples/plot_diagonal.py +++ b/examples/plot_diagonal.py @@ -95,7 +95,7 @@ # 1st dim d = np.arange(nx) -Dop = pylops.Diagonal(d, dims=(nx, ny), dir=0) +Dop = pylops.Diagonal(d, dims=(nx, ny), axis=0) y = Dop * x.ravel() y1 = Dop.H * x.ravel() @@ -105,7 +105,7 @@ # 2nd dim d = np.arange(ny) -Dop = pylops.Diagonal(d, dims=(nx, ny), dir=1) +Dop = pylops.Diagonal(d, dims=(nx, ny), axis=1) y = Dop * x.ravel() y1 = Dop.H * x.ravel() diff --git a/examples/plot_fft.py b/examples/plot_fft.py index 82a5cb1b..c68e6f69 100644 --- a/examples/plot_fft.py +++ b/examples/plot_fft.py @@ -74,7 +74,7 @@ nfft = 2 ** 10 d = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1) -FFTop = pylops.signalprocessing.FFT(dims=(nt, nx), dir=0, nfft=nfft, sampling=dt) +FFTop = pylops.signalprocessing.FFT(dims=(nt, nx), axis=0, nfft=nfft, sampling=dt) D = FFTop * d.ravel() # Adjoint = inverse for FFT diff --git a/examples/plot_flip.py b/examples/plot_flip.py index 013e12bd..483b36f2 100755 --- a/examples/plot_flip.py +++ b/examples/plot_flip.py @@ -40,7 +40,7 @@ # first flip the model along the first axis and then along the second axis nt, nx = 10, 5 x = np.outer(np.arange(nt), np.ones(nx)) -Fop = pylops.Flip((nt, nx), dir=0) +Fop = pylops.Flip((nt, nx), axis=0) y = Fop * x.ravel() xadj = Fop.H * y.ravel() y = y.reshape(nt, nx) @@ -64,7 +64,7 @@ x = np.outer(np.ones(nt), np.arange(nx)) -Fop = pylops.Flip((nt, nx), dir=1) +Fop = pylops.Flip(dims=(nt, nx), axis=1) y = Fop * x.ravel() xadj = Fop.H * y.ravel() y = y.reshape(nt, nx) diff --git a/examples/plot_restriction.py b/examples/plot_restriction.py index 59b21653..051f9c9a 100644 --- a/examples/plot_restriction.py +++ b/examples/plot_restriction.py @@ -95,7 +95,7 @@ nxsub = int(np.round(nx * perc_subsampling)) iava = np.sort(np.random.permutation(np.arange(nx))[:nxsub]) -Rop = pylops.Restriction((nx, nt), iava, dir=0, dtype="float64") +Rop = pylops.Restriction((nx, nt), iava, axis=0, dtype="float64") y = (Rop * x.ravel()).reshape(nxsub, nt) ymask = Rop.mask(x) diff --git a/examples/plot_roll.py b/examples/plot_roll.py index e339b87d..3429bd70 100644 --- a/examples/plot_roll.py +++ b/examples/plot_roll.py @@ -37,7 +37,7 @@ ny, nx = 10, 5 x = np.arange(ny * nx).reshape(ny, nx) -Rop = pylops.Roll((ny, nx), dir=1, shift=-2) +Rop = pylops.Roll(dims=(ny, nx), axis=1, shift=-2) y = Rop * x.ravel() xadj = Rop.H * y diff --git a/examples/plot_shift.py b/examples/plot_shift.py index b02afd71..2b37f52d 100755 --- a/examples/plot_shift.py +++ b/examples/plot_shift.py @@ -52,10 +52,10 @@ shift = 10.5 * dt -# 1st dir +# 1st axis wav2d = np.outer(wav, np.ones(10)) Op = pylops.signalprocessing.Shift( - (nt, 10), shift, dir=0, sampling=dt, real=True, dtype=np.float64 + (nt, 10), shift, axis=0, sampling=dt, real=True, dtype=np.float64 ) wav2dshift = (Op * wav2d.ravel()).reshape(nt, 10) wav2dshiftback = (Op.H * wav2dshift.ravel()).reshape(nt, 10) @@ -72,10 +72,10 @@ axs[2].axis("tight") fig.tight_layout() -# 2nd dir +# 2nd axis wav2d = np.outer(wav, np.ones(10)).T Op = pylops.signalprocessing.Shift( - (10, nt), shift, dir=1, sampling=dt, real=True, dtype=np.float64 + (10, nt), shift, axis=1, sampling=dt, real=True, dtype=np.float64 ) wav2dshift = (Op * wav2d.ravel()).reshape(10, nt) wav2dshiftback = (Op.H * wav2dshift.ravel()).reshape(10, nt) diff --git a/examples/plot_smoothing1d.py b/examples/plot_smoothing1d.py index c0b3086f..57a6c7f9 100755 --- a/examples/plot_smoothing1d.py +++ b/examples/plot_smoothing1d.py @@ -64,7 +64,7 @@ A = np.zeros((11, 21)) A[5, 10] = 1 -Sop = pylops.Smoothing1D(nsmooth=5, dims=(11, 21), dir=0, dtype="float64") +Sop = pylops.Smoothing1D(nsmooth=5, dims=(11, 21), axis=0, dtype="float64") B = np.reshape(Sop * np.ndarray.flatten(A), (11, 21)) fig, axs = plt.subplots(1, 2, figsize=(10, 3)) diff --git a/examples/plot_stacking.py b/examples/plot_stacking.py index 5f423651..2d8077a4 100755 --- a/examples/plot_stacking.py +++ b/examples/plot_stacking.py @@ -27,8 +27,8 @@ ############################################################################### # Let's start by defining two second derivatives :py:class:`pylops.SecondDerivative` # that we will be using in this example. -D2hop = pylops.SecondDerivative((11, 21), dir=1, dtype="float32") -D2vop = pylops.SecondDerivative((11, 21), dir=0, dtype="float32") +D2hop = pylops.SecondDerivative(dims=(11, 21), axis=1, dtype="float32") +D2vop = pylops.SecondDerivative(dims=(11, 21), axis=0, dtype="float32") ############################################################################### # Chaining of operators represents the simplest concatenation that @@ -268,7 +268,7 @@ # :py:class:`pylops.FirstDerivative` to the second dimension of the model. # # Note that for those operators whose implementation allows their application -# to a single axis via the ``dir`` parameter, using the Kronecker product +# to a single axis via the ``axis`` parameter, using the Kronecker product # would lead to slower performance. Nevertheless, the Kronecker product allows # any other operator to be applied to a single dimension. Nv, Nh = 11, 21 diff --git a/examples/plot_sum.py b/examples/plot_sum.py index bd8f6b15..5dc958a1 100644 --- a/examples/plot_sum.py +++ b/examples/plot_sum.py @@ -19,7 +19,7 @@ ############################################################################### # We can now create the operator and peform forward and adjoint -Sop = pylops.Sum(dims=(ny, nx), dir=0) +Sop = pylops.Sum(dims=(ny, nx), axis=0) y = Sop * x.ravel() xadj = Sop.H * y diff --git a/examples/plot_symmetrize.py b/examples/plot_symmetrize.py index 8f882323..7bd4d6f5 100755 --- a/examples/plot_symmetrize.py +++ b/examples/plot_symmetrize.py @@ -49,7 +49,7 @@ nt, nx = 10, 6 x = np.outer(np.arange(nt), np.ones(nx)) -Sop = pylops.Symmetrize((nt, nx), dir=0) +Sop = pylops.Symmetrize((nt, nx), axis=0) y = Sop * x.ravel() xadj = Sop.H * y.ravel() xinv = Sop / y @@ -75,7 +75,7 @@ x = np.outer(np.ones(nt), np.arange(nx)) -Sop = pylops.Symmetrize((nt, nx), dir=1) +Sop = pylops.Symmetrize((nt, nx), axis=1) y = Sop * x.ravel() xadj = Sop.H * y.ravel() diff --git a/examples/plot_tvreg.py b/examples/plot_tvreg.py index 163328e9..56373d5e 100755 --- a/examples/plot_tvreg.py +++ b/examples/plot_tvreg.py @@ -112,7 +112,7 @@ perc_subsampling = 0.6 nxsub = int(np.round(ny * nx * perc_subsampling)) iava = np.sort(np.random.permutation(np.arange(ny * nx))[:nxsub]) -Rop = pylops.Restriction(ny * nx, iava, dtype=np.complex128) +Rop = pylops.Restriction(ny * nx, iava, axis=0, dtype=np.complex128) Fop = pylops.signalprocessing.FFT2D(dims=(ny, nx)) n = np.random.normal(0, 0.0, (ny, nx)) @@ -148,10 +148,10 @@ Dop = [ pylops.FirstDerivative( - (ny, nx), dir=0, edge=False, kind="backward", dtype=np.complex128 + (ny, nx), axis=0, edge=False, kind="backward", dtype=np.complex128 ), pylops.FirstDerivative( - (ny, nx), dir=1, edge=False, kind="backward", dtype=np.complex128 + (ny, nx), axis=1, edge=False, kind="backward", dtype=np.complex128 ), ] diff --git a/pylops/avo/poststack.py b/pylops/avo/poststack.py index cc70dd8f..4384ff2b 100644 --- a/pylops/avo/poststack.py +++ b/pylops/avo/poststack.py @@ -92,7 +92,12 @@ def _PoststackLinearModelling( # Create wavelet operator if len(wav.shape) == 1: Cop = _Convolve1D( - dims, h=wav, offset=len(wav) // 2, dir=0, dtype=dtype, **args_Convolve1D + dims, + h=wav, + offset=len(wav) // 2, + axis=0, + dtype=dtype, + **args_Convolve1D ) else: Cop = _MatrixMult( @@ -103,7 +108,7 @@ def _PoststackLinearModelling( ) # Create derivative operator Dop = _FirstDerivative( - dims, dir=0, sampling=1.0, kind=kind, dtype=dtype, **args_FirstDerivative + dims, axis=0, sampling=1.0, kind=kind, dtype=dtype, **args_FirstDerivative ) Pop = Cop * Dop return Pop @@ -380,7 +385,7 @@ def PoststackInversion( elif dims == 2: Regop = Laplacian((nt0, nx), dtype=PPop.dtype) else: - Regop = Laplacian((nt0, nx, ny), dirs=(1, 2), dtype=PPop.dtype) + Regop = Laplacian((nt0, nx, ny), axes=(1, 2), dtype=PPop.dtype) minv = RegularizedInversion( PPop, @@ -398,14 +403,14 @@ def PoststackInversion( RegL2op = None elif dims == 2: RegL1op = FirstDerivative( - (nt0, nx), dir=0, kind="forward", dtype=PPop.dtype + (nt0, nx), axis=0, kind="forward", dtype=PPop.dtype ) - RegL2op = SecondDerivative((nt0, nx), dir=1, dtype=PPop.dtype) + RegL2op = SecondDerivative((nt0, nx), axis=1, dtype=PPop.dtype) else: RegL1op = FirstDerivative( - (nt0, nx, ny), dir=0, kind="forward", dtype=PPop.dtype + (nt0, nx, ny), axis=0, kind="forward", dtype=PPop.dtype ) - RegL2op = Laplacian((nt0, nx, ny), dirs=(1, 2), dtype=PPop.dtype) + RegL2op = Laplacian((nt0, nx, ny), axes=(1, 2), dtype=PPop.dtype) if "mu" in kwargs_solver.keys(): mu = kwargs_solver["mu"] diff --git a/pylops/avo/prestack.py b/pylops/avo/prestack.py index c3fa152f..a45780a6 100644 --- a/pylops/avo/prestack.py +++ b/pylops/avo/prestack.py @@ -190,7 +190,7 @@ def PrestackLinearModelling( dims, h=wav, offset=len(wav) // 2, - dir=0, + axis=0, dtype=dtype, ) @@ -202,7 +202,7 @@ def PrestackLinearModelling( # Create derivative operator dimsm = list(dims) dimsm[1] = AVOop.npars - Dop = FirstDerivative(dimsm, dir=0, sampling=1.0, kind=kind, dtype=dtype) + Dop = FirstDerivative(dimsm, axis=0, sampling=1.0, kind=kind, dtype=dtype) Preop = Cop * AVOop * Dop return Preop @@ -593,7 +593,7 @@ def PrestackInversion( if isinstance(epsI, (list, tuple)): if len(epsI) != nm: raise ValueError("epsI must be a scalar or a list of" "size nm") - RegI = Diagonal(np.array(epsI), dims=(nt0, nm, nspatprod), dir=1) + RegI = Diagonal(np.array(epsI), dims=(nt0, nm, nspatprod), axis=1) else: RegI = epsI * Identity(nt0 * nm * nspatprod) @@ -602,9 +602,9 @@ def PrestackInversion( if dims == 1: Regop = SecondDerivative((nt0, nm), dtype=PPop.dtype) elif dims == 2: - Regop = Laplacian((nt0, nm, nx), dirs=(0, 2), dtype=PPop.dtype) + Regop = Laplacian((nt0, nm, nx), axes=(0, 2), dtype=PPop.dtype) else: - Regop = Laplacian((nt0, nm, nx, ny), dirs=(2, 3), dtype=PPop.dtype) + Regop = Laplacian((nt0, nm, nx, ny), axes=(2, 3), dtype=PPop.dtype) if epsI is None: Regop = (Regop,) epsR = (epsR,) @@ -626,11 +626,11 @@ def PrestackInversion( RegL1op = FirstDerivative(nt0 * nm, dtype=PPop.dtype) RegL2op = None elif dims == 2: - RegL1op = FirstDerivative((nt0, nm, nx), dir=0, dtype=PPop.dtype) - RegL2op = SecondDerivative((nt0, nm, nx), dir=2, dtype=PPop.dtype) + RegL1op = FirstDerivative((nt0, nm, nx), axis=0, dtype=PPop.dtype) + RegL2op = SecondDerivative((nt0, nm, nx), axis=2, dtype=PPop.dtype) else: - RegL1op = FirstDerivative((nt0, nm, nx, ny), dir=0, dtype=PPop.dtype) - RegL2op = Laplacian((nt0, nm, nx, ny), dirs=(2, 3), dtype=PPop.dtype) + RegL1op = FirstDerivative((nt0, nm, nx, ny), axis=0, dtype=PPop.dtype) + RegL2op = Laplacian((nt0, nm, nx, ny), axes=(2, 3), dtype=PPop.dtype) if dims == 1: if epsI is not None: RegL2op = (RegI,) diff --git a/pylops/basicoperators/CausalIntegration.py b/pylops/basicoperators/CausalIntegration.py index 01ef05ff..12fcc15b 100644 --- a/pylops/basicoperators/CausalIntegration.py +++ b/pylops/basicoperators/CausalIntegration.py @@ -1,3 +1,5 @@ +import warnings + import numpy as np from pylops import LinearOperator @@ -7,14 +9,16 @@ class CausalIntegration(LinearOperator): r"""Causal integration. - Apply causal integration to a multi-dimensional array along ``dir`` axis. + Apply causal integration to a multi-dimensional array along ``axis``. Parameters ---------- dims : :obj:`list` or :obj:`int` Number of samples for each dimension - dir : :obj:`int`, optional - Direction along which array is integrated. + axis : :obj:`int`, optional + .. versionadded:: 2.0.0 + + Axis along which the model is integrated. sampling : :obj:`float`, optional Sampling step ``dx``. halfcurrent : :obj:`bool`, optional @@ -86,7 +90,7 @@ class CausalIntegration(LinearOperator): def __init__( self, dims, - dir=-1, + axis=-1, sampling=1, halfcurrent=True, dtype="float64", @@ -94,7 +98,7 @@ def __init__( removefirst=False, ): self.dims = _value_or_list_like_to_array(dims) - self.dir = dir + self.axis = axis self.sampling = sampling self.kind = kind if kind == "full" and halfcurrent: # ensure backcompatibility @@ -102,15 +106,15 @@ def __init__( self.removefirst = removefirst self.dimsd = self.dims.copy() if self.removefirst: - self.dimsd[self.dir] -= 1 + self.dimsd[self.axis] -= 1 self.shape = (np.prod(self.dimsd), np.prod(self.dims)) self.dtype = np.dtype(dtype) self.explicit = False def _matvec(self, x): x = np.reshape(x, self.dims) - if self.dir != -1: - x = np.swapaxes(x, self.dir, -1) + if self.axis != -1: + x = np.swapaxes(x, self.axis, -1) y = self.sampling * np.cumsum(x, axis=-1) if self.kind in ("half", "trapezoidal"): y -= self.sampling * x / 2.0 @@ -118,16 +122,16 @@ def _matvec(self, x): y[..., 1:] -= self.sampling * x[..., 0:1] / 2.0 if self.removefirst: y = y[..., 1:] - if self.dir != -1: - y = np.swapaxes(y, -1, self.dir) + if self.axis != -1: + y = np.swapaxes(y, -1, self.axis) return y.ravel() def _rmatvec(self, x): x = np.reshape(x, self.dimsd) if self.removefirst: - x = np.insert(x, 0, 0, axis=self.dir) - if self.dir != -1: - x = np.swapaxes(x, self.dir, -1) + x = np.insert(x, 0, 0, axis=self.axis) + if self.axis != -1: + x = np.swapaxes(x, self.axis, -1) xflip = np.flip(x, axis=-1) if self.kind == "half": y = self.sampling * (np.cumsum(xflip, axis=-1) - xflip / 2.0) @@ -137,6 +141,6 @@ def _rmatvec(self, x): else: y = self.sampling * np.cumsum(xflip, axis=-1) y = np.flip(y, axis=-1) - if self.dir != -1: - y = np.swapaxes(y, -1, self.dir) + if self.axis != -1: + y = np.swapaxes(y, -1, self.axis) return y.ravel() diff --git a/pylops/basicoperators/Diagonal.py b/pylops/basicoperators/Diagonal.py index f1a4cd4d..ebfd20e5 100644 --- a/pylops/basicoperators/Diagonal.py +++ b/pylops/basicoperators/Diagonal.py @@ -1,3 +1,5 @@ +import warnings + import numpy as np from pylops import LinearOperator @@ -12,8 +14,8 @@ class Diagonal(LinearOperator): This operator can also broadcast; in this case the input vector is reshaped into its dimensions ``dims`` and the element-wise multiplication - with ``diag`` is perfomed on the direction ``dir``. Note that the - vector ``diag`` will need to have size equal to ``dims[dir]``. + with ``diag`` is perfomed along ``axis``. Note that the + vector ``diag`` will need to have size equal to ``dims[axis]``. Parameters ---------- @@ -22,8 +24,10 @@ class Diagonal(LinearOperator): dims : :obj:`list`, optional Number of samples for each dimension (``None`` if only one dimension is available) - dir : :obj:`int`, optional - Direction along which multiplication is applied. + axis : :obj:`int`, optional + .. versionadded:: 2.0.0 + + Axis along which multiplication is applied. dtype : :obj:`str`, optional Type of elements in input array. @@ -55,17 +59,18 @@ class Diagonal(LinearOperator): """ - def __init__(self, diag, dims=None, dir=0, dtype="float64"): + def __init__(self, diag, dims=None, axis=-1, dtype="float64"): ncp = get_array_module(diag) self.diag = diag.ravel() self.complex = True if ncp.iscomplexobj(self.diag) else False + if dims is None: self.shape = (len(self.diag), len(self.diag)) self.dims = None self.reshape = False else: diagdims = [1] * len(dims) - diagdims[dir] = dims[dir] + diagdims[axis] = dims[axis] self.diag = self.diag.reshape(diagdims) self.shape = (np.prod(dims), np.prod(dims)) self.dims = dims diff --git a/pylops/basicoperators/DirectionalDerivative.py b/pylops/basicoperators/DirectionalDerivative.py index 7fccfc08..f9c4c728 100644 --- a/pylops/basicoperators/DirectionalDerivative.py +++ b/pylops/basicoperators/DirectionalDerivative.py @@ -8,7 +8,7 @@ def FirstDirectionalDerivative( r"""First Directional derivative. Apply a directional derivative operator to a multi-dimensional array - along either a single common direction or different directions for each + along either a single common axis or different axes for each point of the array. .. note:: At least 2 dimensions are required, consider using @@ -63,10 +63,10 @@ def FirstDirectionalDerivative( """ Gop = Gradient(dims, sampling=sampling, edge=edge, kind=kind, dtype=dtype) if v.ndim == 1: - Dop = Diagonal(v, dims=[len(dims)] + list(dims), dir=0, dtype=dtype) + Dop = Diagonal(v, dims=[len(dims)] + list(dims), axis=0, dtype=dtype) else: Dop = Diagonal(v.ravel(), dtype=dtype) - Sop = Sum(dims=[len(dims)] + list(dims), dir=0, dtype=dtype) + Sop = Sum(dims=[len(dims)] + list(dims), axis=0, dtype=dtype) ddop = Sop * Dop * Gop return LinearOperator(ddop) @@ -75,7 +75,7 @@ def SecondDirectionalDerivative(dims, v, sampling=1, edge=False, dtype="float64" r"""Second Directional derivative. Apply a second directional derivative operator to a multi-dimensional array - along either a single common direction or different directions for each + along either a single common axis or different axes for each point of the array. .. note:: At least 2 dimensions are required, consider using diff --git a/pylops/basicoperators/FirstDerivative.py b/pylops/basicoperators/FirstDerivative.py index c3c82c1e..7ae177a1 100644 --- a/pylops/basicoperators/FirstDerivative.py +++ b/pylops/basicoperators/FirstDerivative.py @@ -1,3 +1,5 @@ +import warnings + import numpy as np from numpy.core.multiarray import normalize_axis_index @@ -16,8 +18,10 @@ class FirstDerivative(LinearOperator): ---------- dims : :obj:`list` or :obj:`int` Number of samples for each dimension - dir : :obj:`int`, optional - Direction along which the derivative is applied. + axis : :obj:`int`, optional + .. versionadded:: 2.0.0 + + Axis along which derivative is applied. sampling : :obj:`float`, optional Sampling step :math:`\Delta x`. edge : :obj:`bool`, optional @@ -63,16 +67,16 @@ class FirstDerivative(LinearOperator): def __init__( self, dims, - dir=0, + axis=-1, sampling=1.0, edge=False, dtype="float64", kind="centered", ): self.dims = _value_or_list_like_to_array(dims) + self.axis = normalize_axis_index(axis, len(self.dims)) self.sampling = sampling self.edge = edge - self.dir = normalize_axis_index(dir, len(self.dims)) self.kind = kind N = np.prod(self.dims) self.shape = (N, N) @@ -95,50 +99,50 @@ def __init__( def _matvec_forward(self, x): ncp = get_array_module(x) x = ncp.reshape(x, self.dims) - if self.dir > 0: # need to bring the dim. to derive to first dim. - x = ncp.swapaxes(x, self.dir, 0) + if self.axis > 0: # need to bring the dim. to derive to first dim. + x = ncp.swapaxes(x, self.axis, 0) y = ncp.zeros(x.shape, self.dtype) y[:-1, ...] = (x[1:, ...] - x[:-1, ...]) / self.sampling - if self.dir > 0: - y = ncp.swapaxes(y, 0, self.dir) + if self.axis > 0: + y = ncp.swapaxes(y, 0, self.axis) y = y.ravel() return y def _rmatvec_forward(self, x): ncp = get_array_module(x) x = ncp.reshape(x, self.dims) - if self.dir > 0: # need to bring the dim. to derive to first dim. - x = ncp.swapaxes(x, self.dir, 0) + if self.axis > 0: # need to bring the dim. to derive to first dim. + x = ncp.swapaxes(x, self.axis, 0) y = ncp.zeros(x.shape, self.dtype) y[:-1, ...] -= x[:-1, ...] y[1:, ...] += x[:-1, ...] y /= self.sampling - if self.dir > 0: - y = ncp.swapaxes(y, 0, self.dir) + if self.axis > 0: + y = ncp.swapaxes(y, 0, self.axis) y = y.ravel() return y def _matvec_centered(self, x): ncp = get_array_module(x) x = ncp.reshape(x, self.dims) - if self.dir > 0: # need to bring the dim. to derive to first dim. - x = ncp.swapaxes(x, self.dir, 0) + if self.axis > 0: # need to bring the dim. to derive to first dim. + x = ncp.swapaxes(x, self.axis, 0) y = ncp.zeros(x.shape, self.dtype) y[1:-1, ...] = 0.5 * x[2:, ...] - 0.5 * x[0:-2, ...] if self.edge: y[0, ...] = x[1, ...] - x[0, ...] y[-1, ...] = x[-1, ...] - x[-2, ...] y /= self.sampling - if self.dir > 0: - y = ncp.swapaxes(y, 0, self.dir) + if self.axis > 0: + y = ncp.swapaxes(y, 0, self.axis) y = y.ravel() return y def _rmatvec_centered(self, x): ncp = get_array_module(x) x = ncp.reshape(x, self.dims) - if self.dir > 0: # need to bring the dim. to derive to first dim. - x = ncp.swapaxes(x, self.dir, 0) + if self.axis > 0: # need to bring the dim. to derive to first dim. + x = ncp.swapaxes(x, self.axis, 0) y = ncp.zeros(x.shape, self.dtype) y[0:-2, ...] -= 0.5 * x[1:-1, ...] y[2:, ...] += 0.5 * x[1:-1, ...] @@ -148,33 +152,33 @@ def _rmatvec_centered(self, x): y[-2, ...] -= x[-1, ...] y[-1, ...] += x[-1, ...] y /= self.sampling - if self.dir > 0: - y = ncp.swapaxes(y, 0, self.dir) + if self.axis > 0: + y = ncp.swapaxes(y, 0, self.axis) y = y.ravel() return y def _matvec_backward(self, x): ncp = get_array_module(x) x = ncp.reshape(x, self.dims) - if self.dir > 0: # need to bring the dim. to derive to first dim. - x = ncp.swapaxes(x, self.dir, 0) + if self.axis > 0: # need to bring the dim. to derive to first dim. + x = ncp.swapaxes(x, self.axis, 0) y = ncp.zeros(x.shape, self.dtype) y[1:, ...] = (x[1:, ...] - x[:-1, ...]) / self.sampling - if self.dir > 0: - y = ncp.swapaxes(y, 0, self.dir) + if self.axis > 0: + y = ncp.swapaxes(y, 0, self.axis) y = y.ravel() return y def _rmatvec_backward(self, x): ncp = get_array_module(x) x = ncp.reshape(x, self.dims) - if self.dir > 0: # need to bring the dim. to derive to first dim. - x = ncp.swapaxes(x, self.dir, 0) + if self.axis > 0: # need to bring the dim. to derive to first dim. + x = ncp.swapaxes(x, self.axis, 0) y = ncp.zeros(x.shape, self.dtype) y[:-1, ...] -= x[1:, ...] y[1:, ...] += x[1:, ...] y /= self.sampling - if self.dir > 0: - y = ncp.swapaxes(y, 0, self.dir) + if self.axis > 0: + y = ncp.swapaxes(y, 0, self.axis) y = y.ravel() return y diff --git a/pylops/basicoperators/Flip.py b/pylops/basicoperators/Flip.py index 17762763..6b816e2e 100644 --- a/pylops/basicoperators/Flip.py +++ b/pylops/basicoperators/Flip.py @@ -1,3 +1,5 @@ +import warnings + import numpy as np from pylops import LinearOperator @@ -7,14 +9,16 @@ class Flip(LinearOperator): r"""Flip along an axis. - Flip a multi-dimensional array along a specified direction ``dir``. + Flip a multi-dimensional array along ``axis``. Parameters ---------- dims : :obj:`list` or :obj:`int` Number of samples for each dimension - dir : :obj:`int`, optional - Direction along which flipping is applied. + axis : :obj:`int`, optional + .. versionadded:: 2.0.0 + + Axis along which model is flipped. dtype : :obj:`str`, optional Type of elements in input array. @@ -35,15 +39,15 @@ class Flip(LinearOperator): .. math:: y[i] = x[N-1-i] \quad \forall i=0,1,2,\ldots,N-1 - where :math:`N` is the dimension of the input model along ``dir``. As this operator is + where :math:`N` is the dimension of the input model along ``axis``. As this operator is self-adjoint, :math:`x` and :math:`y` in the equation above are simply swapped in adjoint mode. """ - def __init__(self, dims, dir=0, dtype="float64"): - self.dir = dir + def __init__(self, dims, axis=-1, dtype="float64"): self.dims = _value_or_list_like_to_array(dims) + self.axis = axis N = np.prod(self.dims) self.shape = (N, N) self.dtype = np.dtype(dtype) @@ -51,7 +55,7 @@ def __init__(self, dims, dir=0, dtype="float64"): def _matvec(self, x): x = np.reshape(x, self.dims) - y = np.flip(x, axis=self.dir) + y = np.flip(x, axis=self.axis) y = y.ravel() return y diff --git a/pylops/basicoperators/Gradient.py b/pylops/basicoperators/Gradient.py index 4f1c12f5..fbef8e60 100644 --- a/pylops/basicoperators/Gradient.py +++ b/pylops/basicoperators/Gradient.py @@ -66,13 +66,13 @@ def Gradient(dims, sampling=1, edge=False, dtype="float64", kind="centered"): [ FirstDerivative( dims, - dir=idir, - sampling=sampling[idir], + axis=iax, + sampling=sampling[iax], edge=edge, kind=kind, dtype=dtype, ) - for idir in range(ndims) + for iax in range(ndims) ] ) return gop diff --git a/pylops/basicoperators/Laplacian.py b/pylops/basicoperators/Laplacian.py index 0ea56074..d854d0df 100644 --- a/pylops/basicoperators/Laplacian.py +++ b/pylops/basicoperators/Laplacian.py @@ -1,4 +1,7 @@ +import warnings + import numpy as np +from numpy.core.multiarray import normalize_axis_index from pylops.basicoperators import SecondDerivative from pylops.LinearOperator import aslinearoperator @@ -6,7 +9,7 @@ def Laplacian( dims, - dirs=(0, 1), + axes=(-2, -1), weights=(1, 1), sampling=(1, 1), edge=False, @@ -24,8 +27,10 @@ def Laplacian( ---------- dims : :obj:`tuple` Number of samples for each dimension. - dirs : :obj:`tuple`, optional - Directions along which laplacian is applied. + axes : :obj:`int`, optional + .. versionadded:: 2.0.0 + + Axes along which the Laplacian is applied. weights : :obj:`tuple`, optional Weight to apply to each direction (real laplacian operator if ``weights=(1, 1)``) @@ -47,7 +52,7 @@ def Laplacian( Raises ------ ValueError - If ``dirs``. ``weights``, and ``sampling`` do not have the same size + If ``axes``. ``weights``, and ``sampling`` do not have the same size Notes ----- @@ -61,16 +66,17 @@ def Laplacian( / (\Delta x \Delta y) """ - if not (len(dirs) == len(weights) == len(sampling)): - raise ValueError("dirs, weights, and sampling have different size") + axes = tuple(normalize_axis_index(ax, len(dims)) for ax in axes) + if not (len(axes) == len(weights) == len(sampling)): + raise ValueError("axes, weights, and sampling have different size") l2op = weights[0] * SecondDerivative( - dims, dir=dirs[0], sampling=sampling[0], edge=edge, kind=kind, dtype=dtype + dims, axis=axes[0], sampling=sampling[0], edge=edge, kind=kind, dtype=dtype ) - for dir, samp, weight in zip(dirs[1:], sampling[1:], weights[1:]): + for ax, samp, weight in zip(axes[1:], sampling[1:], weights[1:]): l2op += weight * SecondDerivative( - dims, dir=dir, sampling=samp, edge=edge, dtype=dtype + dims, axis=ax, sampling=samp, edge=edge, dtype=dtype ) return aslinearoperator(l2op) diff --git a/pylops/basicoperators/Restriction.py b/pylops/basicoperators/Restriction.py index 9c8e6c55..9091574f 100644 --- a/pylops/basicoperators/Restriction.py +++ b/pylops/basicoperators/Restriction.py @@ -1,3 +1,5 @@ +import warnings + import numpy as np import numpy.ma as np_ma from numpy.core.multiarray import normalize_axis_index @@ -7,14 +9,14 @@ from pylops.utils.backend import get_array_module, to_cupy_conditional -def _compute_iavamask(dims, dir, iava, ncp): +def _compute_iavamask(dims, axis, iava, ncp): """Compute restriction mask when using cupy arrays""" otherdims = np.array(dims) - otherdims = np.delete(otherdims, dir) - iavamask = ncp.zeros(dims[dir], dtype=int) + otherdims = np.delete(otherdims, axis) + iavamask = ncp.zeros(dims[axis], dtype=int) iavamask[iava] = 1 iavamask = ncp.moveaxis( - ncp.broadcast_to(iavamask, list(otherdims) + [dims[dir]]), -1, dir + ncp.broadcast_to(iavamask, list(otherdims) + [dims[axis]]), -1, axis ) iavamask = ncp.where(iavamask.ravel() == 1)[0] return iavamask @@ -33,8 +35,10 @@ class Restriction(LinearOperator): Number of samples for each dimension iava : :obj:`list` or :obj:`numpy.ndarray` Integer indices of available samples for data selection. - dir : :obj:`int`, optional - Direction along which restriction is applied. + axis : :obj:`int`, optional + .. versionadded:: 2.0.0 + + Axis along which restriction is applied to model. dtype : :obj:`str`, optional Type of elements in input array. inplace : :obj:`bool`, optional @@ -80,21 +84,21 @@ class Restriction(LinearOperator): """ - def __init__(self, dims, iava, dir=0, dtype="float64", inplace=True): + def __init__(self, dims, iava, axis=-1, dtype="float64", inplace=True): ncp = get_array_module(iava) self.dims = _value_or_list_like_to_array(dims) - self.dir = normalize_axis_index(dir, len(self.dims)) + self.axis = normalize_axis_index(axis, len(self.dims)) self.iava = iava self.dimsd = self.dims.copy() # data dimensions - self.dimsd[self.dir] = len(iava) + self.dimsd[self.axis] = len(iava) self.iavareshape = np.ones(len(self.dims), dtype=int) - self.iavareshape[self.dir] = len(self.iava) + self.iavareshape[self.axis] = len(self.iava) # currently cupy does not support put_along_axis, so we need to # explicitely create a list of indices in the n-dimensional # model space which will be used in _rmatvec to place the input if ncp != np: - self.iavamask = _compute_iavamask(dims, dir, iava, ncp) + self.iavamask = _compute_iavamask(self.dims, self.axis, self.iava, ncp) self.inplace = inplace self.shape = (np.prod(self.dimsd), np.prod(self.dims)) self.dtype = np.dtype(dtype) @@ -105,7 +109,7 @@ def _matvec(self, x): if not self.inplace: x = x.copy() x = ncp.reshape(x, self.dims) - y = ncp.take(x, self.iava, axis=self.dir) + y = ncp.take(x, self.iava, axis=self.axis) y = y.ravel() return y @@ -117,12 +121,12 @@ def _rmatvec(self, x): if ncp == np: y = ncp.zeros(self.dims, dtype=self.dtype) ncp.put_along_axis( - y, ncp.reshape(self.iava, self.iavareshape), x, axis=self.dir + y, ncp.reshape(self.iava, self.iavareshape), x, axis=self.axis ) else: if not hasattr(self, "iavamask"): self.iava = to_cupy_conditional(x, self.iava) - self.iavamask = _compute_iavamask(self.dims, self.dir, self.iava, ncp) + self.iavamask = _compute_iavamask(self.dims, self.axis, self.iava, ncp) y = ncp.zeros(self.shape[-1], dtype=self.dtype) y[self.iavamask] = x.ravel() y = y.ravel() @@ -151,14 +155,14 @@ def mask(self, x): y = np_ma.array(np.zeros(self.dims), mask=np.ones(self.dims), dtype=self.dtype) x = np.reshape(x, self.dims) - if self.dir > 0: - x = np.swapaxes(x, self.dir, 0) - y = np.swapaxes(y, self.dir, 0) + if self.axis > 0: + x = np.swapaxes(x, self.axis, 0) + y = np.swapaxes(y, self.axis, 0) y.mask[iava] = False if ncp == np: y[iava] = x[self.iava] else: y[iava] = ncp.asnumpy(x)[iava] - if self.dir > 0: - y = np.swapaxes(y, 0, self.dir) + if self.axis > 0: + y = np.swapaxes(y, 0, self.axis) return y diff --git a/pylops/basicoperators/Roll.py b/pylops/basicoperators/Roll.py index 94852c15..29ca3681 100644 --- a/pylops/basicoperators/Roll.py +++ b/pylops/basicoperators/Roll.py @@ -1,3 +1,5 @@ +import warnings + import numpy as np from pylops import LinearOperator @@ -7,15 +9,17 @@ class Roll(LinearOperator): r"""Roll along an axis. - Roll a multi-dimensional array along a specified direction ``dir`` for + Roll a multi-dimensional array along ``axis`` for a chosen number of samples (``shift``). Parameters ---------- dims : :obj:`list` or :obj:`int` Number of samples for each dimension - dir : :obj:`int`, optional - Direction along which rolling is applied. + axis : :obj:`int`, optional + .. versionadded:: 2.0.0 + + Axis along which model is rolled. shift : :obj:`int`, optional Number of samples by which elements are shifted dtype : :obj:`str`, optional @@ -37,10 +41,10 @@ class Roll(LinearOperator): """ - def __init__(self, dims, dir=0, shift=1, dtype="float64"): + def __init__(self, dims, axis=-1, shift=1, dtype="float64"): self.dims = _value_or_list_like_to_array(dims) + self.axis = axis N = np.prod(self.dims) - self.dir = dir self.shift = shift self.shape = (N, N) self.dtype = np.dtype(dtype) @@ -48,10 +52,10 @@ def __init__(self, dims, dir=0, shift=1, dtype="float64"): def _matvec(self, x): x = np.reshape(x, self.dims) - y = np.roll(x, shift=self.shift, axis=self.dir) + y = np.roll(x, shift=self.shift, axis=self.axis) return y.ravel() def _rmatvec(self, x): x = np.reshape(x, self.dims) - y = np.roll(x, shift=-self.shift, axis=self.dir) + y = np.roll(x, shift=-self.shift, axis=self.axis) return y.ravel() diff --git a/pylops/basicoperators/SecondDerivative.py b/pylops/basicoperators/SecondDerivative.py index 8689f1a7..f2437aa5 100644 --- a/pylops/basicoperators/SecondDerivative.py +++ b/pylops/basicoperators/SecondDerivative.py @@ -1,3 +1,5 @@ +import warnings + import numpy as np from numpy.core.multiarray import normalize_axis_index @@ -10,14 +12,17 @@ class SecondDerivative(LinearOperator): r"""Second derivative. Apply a second derivative using a three-point stencil finite-difference - approximation. + approximation along ``axis``. Parameters ---------- dims : :obj:`list` or :obj:`int` Number of samples for each dimension - dir : :obj:`int`, optional - Direction along which the derivative is applied. + (``None`` if only one dimension is available) + axis : :obj:`int`, optional + .. versionadded:: 2.0.0 + + Axis along which derivative is applied. sampling : :obj:`float`, optional Sampling step :math:`\Delta x`. edge : :obj:`bool`, optional @@ -60,12 +65,18 @@ class SecondDerivative(LinearOperator): """ def __init__( - self, dims, dir=0, sampling=1, edge=False, dtype="float64", kind="centered" + self, + dims, + axis=-1, + sampling=1, + edge=False, + dtype="float64", + kind="centered", ): self.dims = _value_or_list_like_to_array(dims) + self.axis = normalize_axis_index(axis, len(self.dims)) self.sampling = sampling self.edge = edge - self.dir = normalize_axis_index(dir, len(self.dims)) self.kind = kind N = np.prod(self.dims) self.shape = (N, N) @@ -89,52 +100,52 @@ def __init__( def _matvec_forward(self, x): ncp = get_array_module(x) x = ncp.reshape(x, self.dims) - if self.dir > 0: # need to bring the dim. to derive to first dim. - x = ncp.swapaxes(x, self.dir, 0) + if self.axis > 0: # need to bring the dim. to derive to first dim. + x = ncp.swapaxes(x, self.axis, 0) y = ncp.zeros(x.shape, self.dtype) y[:-2, ...] = x[2:, ...] - 2 * x[1:-1, ...] + x[0:-2, ...] y /= self.sampling ** 2 - if self.dir > 0: - y = ncp.swapaxes(y, 0, self.dir) + if self.axis > 0: + y = ncp.swapaxes(y, 0, self.axis) y = y.ravel() return y def _rmatvec_forward(self, x): ncp = get_array_module(x) x = ncp.reshape(x, self.dims) - if self.dir > 0: # need to bring the dim. to derive to first dim. - x = ncp.swapaxes(x, self.dir, 0) + if self.axis > 0: # need to bring the dim. to derive to first dim. + x = ncp.swapaxes(x, self.axis, 0) y = ncp.zeros(x.shape, self.dtype) y[0:-2, ...] += x[:-2, ...] y[1:-1, ...] -= 2 * x[:-2, ...] y[2:, ...] += x[:-2, ...] y /= self.sampling ** 2 - if self.dir > 0: - y = ncp.swapaxes(y, 0, self.dir) + if self.axis > 0: + y = ncp.swapaxes(y, 0, self.axis) y = y.ravel() return y def _matvec_centered(self, x): ncp = get_array_module(x) x = ncp.reshape(x, self.dims) - if self.dir > 0: # need to bring the dim. to derive to first dim. - x = ncp.swapaxes(x, self.dir, 0) + if self.axis > 0: # need to bring the dim. to derive to first dim. + x = ncp.swapaxes(x, self.axis, 0) y = ncp.zeros(x.shape, self.dtype) y[1:-1, ...] = x[2:, ...] - 2 * x[1:-1, ...] + x[0:-2, ...] if self.edge: y[0, ...] = x[0, ...] - 2 * x[1, ...] + x[2, ...] y[-1, ...] = x[-3, ...] - 2 * x[-2, ...] + x[-1, ...] y /= self.sampling ** 2 - if self.dir > 0: - y = ncp.swapaxes(y, 0, self.dir) + if self.axis > 0: + y = ncp.swapaxes(y, 0, self.axis) y = y.ravel() return y def _rmatvec_centered(self, x): ncp = get_array_module(x) x = ncp.reshape(x, self.dims) - if self.dir > 0: # need to bring the dim. to derive to first dim. - x = ncp.swapaxes(x, self.dir, 0) + if self.axis > 0: # need to bring the dim. to derive to first dim. + x = ncp.swapaxes(x, self.axis, 0) y = ncp.zeros(x.shape, self.dtype) y[0:-2, ...] += x[1:-1, ...] y[1:-1, ...] -= 2 * x[1:-1, ...] @@ -147,7 +158,7 @@ def _rmatvec_centered(self, x): y[-2, ...] -= 2 * x[-1, ...] y[-1, ...] += x[-1, ...] y /= self.sampling ** 2 - if self.dir > 0: - y = ncp.swapaxes(y, 0, self.dir) + if self.axis > 0: + y = ncp.swapaxes(y, 0, self.axis) y = y.ravel() return y diff --git a/pylops/basicoperators/Smoothing1D.py b/pylops/basicoperators/Smoothing1D.py index e742dddf..037550d2 100644 --- a/pylops/basicoperators/Smoothing1D.py +++ b/pylops/basicoperators/Smoothing1D.py @@ -1,22 +1,26 @@ +import warnings + import numpy as np from pylops.signalprocessing import Convolve1D -def Smoothing1D(nsmooth, dims, dir=0, dtype="float64"): +def Smoothing1D(nsmooth, dims, axis=-1, dtype="float64"): r"""1D Smoothing. - Apply smoothing to model (and data) along a specific direction of a - multi-dimensional array depending on the choice of ``dir``. + Apply smoothing to model (and data) to a multi-dimensional array + along ``axis``. Parameters ---------- nsmooth : :obj:`int` - Lenght of smoothing operator (must be odd) + Length of smoothing operator (must be odd) dims : :obj:`tuple` or :obj:`int` Number of samples for each dimension - dir : :obj:`int`, optional - Direction along which smoothing is applied + axis : :obj:`int`, optional + .. versionadded:: 2.0.0 + + Axis along which model (and data) are smoothed. dtype : :obj:`str`, optional Type of elements in input array. @@ -61,11 +65,6 @@ def Smoothing1D(nsmooth, dims, dir=0, dtype="float64"): """ if nsmooth % 2 == 0: nsmooth += 1 - - return Convolve1D( - dims, - np.ones(nsmooth) / float(nsmooth), - dir=dir, - offset=(nsmooth - 1) / 2, - dtype=dtype, - ) + h = np.ones(nsmooth) / float(nsmooth) + offset = (nsmooth - 1) / 2 + return Convolve1D(dims, h, axis=axis, offset=offset, dtype=dtype) diff --git a/pylops/basicoperators/Smoothing2D.py b/pylops/basicoperators/Smoothing2D.py index 13d8c305..1ecebcaf 100644 --- a/pylops/basicoperators/Smoothing2D.py +++ b/pylops/basicoperators/Smoothing2D.py @@ -1,13 +1,16 @@ +import warnings + import numpy as np +from numpy.core.multiarray import normalize_axis_index from pylops.signalprocessing import Convolve2D -def Smoothing2D(nsmooth, dims, nodir=None, dtype="float64"): +def Smoothing2D(nsmooth, dims, axes=(-2, -1), dtype="float64"): r"""2D Smoothing. - Apply smoothing to model (and data) along two directions of a - multi-dimensional array depending on the choice of ``nodir``. + Apply smoothing to model (and data) along two ``axes`` of a + multi-dimensional array. Parameters ---------- @@ -15,9 +18,10 @@ def Smoothing2D(nsmooth, dims, nodir=None, dtype="float64"): Lenght of smoothing operator in 1st and 2nd dimensions (must be odd) dims : :obj:`tuple` Number of samples for each dimension - nodir : :obj:`int`, optional - Direction along which smoothing is **not** applied (set to ``None`` for 2d - arrays) + axes : :obj:`int`, optional + .. versionadded:: 2.0.0 + + Axes along which model (and data) are smoothed. dtype : :obj:`str`, optional Type of elements in input array. @@ -56,12 +60,6 @@ def Smoothing2D(nsmooth, dims, nodir=None, dtype="float64"): nsmooth[0] += 1 if nsmooth[1] % 2 == 0: nsmooth[1] += 1 - h = np.ones((nsmooth[0], nsmooth[1])) / float(nsmooth[0] * nsmooth[1]) - return Convolve2D( - dims, - h=h, - offset=[(nsmooth[0] - 1) / 2, (nsmooth[1] - 1) / 2], - nodir=nodir, - dtype=dtype, - ) + offset = [(nsmooth[0] - 1) / 2, (nsmooth[1] - 1) / 2] + return Convolve2D(dims, h=h, offset=offset, axes=axes, dtype=dtype) diff --git a/pylops/basicoperators/Sum.py b/pylops/basicoperators/Sum.py index c20af8b1..2ed0f9d2 100644 --- a/pylops/basicoperators/Sum.py +++ b/pylops/basicoperators/Sum.py @@ -1,3 +1,5 @@ +import warnings + import numpy as np from pylops import LinearOperator @@ -7,7 +9,7 @@ class Sum(LinearOperator): r"""Sum operator. - Sum along an axis of a multi-dimensional + Sum along ``axis`` of a multi-dimensional array (at least 2 dimensions are required) in forward model, and spread along the same axis in adjoint mode. @@ -15,8 +17,10 @@ class Sum(LinearOperator): ---------- dims : :obj:`tuple` Number of samples for each dimension - dir : :obj:`int` - Direction along which summation is performed. + axis : :obj:`int`, optional + .. versionadded:: 2.0.0 + + Axis along which model is summed. dtype : :obj:`str`, optional Type of elements in input array. @@ -32,7 +36,7 @@ class Sum(LinearOperator): ----- Given a two dimensional array, the *Sum* operator re-arranges the input model into a multi-dimensional array - of size ``dims`` and sums values along direction ``dir``: + of size ``dims`` and sums values along ``axis``: .. math:: @@ -46,17 +50,17 @@ class Sum(LinearOperator): """ - def __init__(self, dims, dir, dtype="float64"): + def __init__(self, dims, axis=-1, dtype="float64"): if len(dims) == 1: dims = (dims[0], 1) # to avoid reducing matvec to a scalar self.dims = dims - self.dir = dir + self.axis = axis # data dimensions self.dims_d = list(dims).copy() - self.dims_d.pop(dir) - # array of ones with dims of model in dir for np.tile in adjoint mode + self.dims_d.pop(self.axis) + # array of ones with dims of model in self.axis for np.tile in adjoint mode self.tile = np.ones(len(dims), dtype=int) - self.tile[dir] = self.dims[dir] + self.tile[self.axis] = self.dims[self.axis] self.dtype = np.dtype(dtype) self.shape = (np.prod(self.dims_d), np.prod(dims)) self.explicit = False @@ -64,12 +68,12 @@ def __init__(self, dims, dir, dtype="float64"): def _matvec(self, x): ncp = get_array_module(x) y = x.reshape(self.dims) - y = ncp.sum(y, axis=self.dir) + y = ncp.sum(y, axis=self.axis) return y.ravel() def _rmatvec(self, x): ncp = get_array_module(x) y = x.reshape(self.dims_d) - y = ncp.expand_dims(y, self.dir) + y = ncp.expand_dims(y, self.axis) y = ncp.tile(y, self.tile) return y.ravel() diff --git a/pylops/basicoperators/Symmetrize.py b/pylops/basicoperators/Symmetrize.py index 33c1eb36..c436fbab 100644 --- a/pylops/basicoperators/Symmetrize.py +++ b/pylops/basicoperators/Symmetrize.py @@ -1,3 +1,5 @@ +import warnings + import numpy as np from pylops import LinearOperator @@ -8,14 +10,17 @@ class Symmetrize(LinearOperator): r"""Symmetrize along an axis. - Symmetrize a multi-dimensional array along a specified direction ``dir``. + Symmetrize a multi-dimensional array along ``axis``. Parameters ---------- dims : :obj:`list` or :obj:`int` Number of samples for each dimension - dir : :obj:`int`, optional - Direction along which symmetrization is applied + (``None`` if only one dimension is available) + axis : :obj:`int`, optional + .. versionadded:: 2.0.0 + + Axis along which model is symmetrized. dtype : :obj:`str`, optional Type of elements in input array @@ -55,12 +60,12 @@ class Symmetrize(LinearOperator): apart from the central sample where :math:`x[0] = y[N-1]`. """ - def __init__(self, dims, dir=0, dtype="float64"): + def __init__(self, dims, axis=-1, dtype="float64"): self.dims = _value_or_list_like_to_array(dims) - self.dir = dir + self.axis = axis self.dimsd = self.dims.copy() - self.dimsd[self.dir] = self.dims[self.dir] * 2 - 1 - self.nsym = self.dims[self.dir] + self.dimsd[self.axis] = self.dims[self.axis] * 2 - 1 + self.nsym = self.dims[self.axis] self.shape = (np.prod(self.dimsd), np.prod(self.dims)) self.dtype = np.dtype(dtype) self.explicit = False @@ -69,24 +74,24 @@ def _matvec(self, x): ncp = get_array_module(x) y = ncp.zeros(self.dimsd, dtype=self.dtype) x = ncp.reshape(x, self.dims) - if self.dir > 0: # bring the dimension to symmetrize to first - x = ncp.swapaxes(x, self.dir, 0) - y = ncp.swapaxes(y, self.dir, 0) + if self.axis > 0: # bring the dimension to symmetrize to first + x = ncp.swapaxes(x, self.axis, 0) + y = ncp.swapaxes(y, self.axis, 0) y[self.nsym - 1 :] = x y[: self.nsym - 1] = x[-1:0:-1] - if self.dir > 0: - y = ncp.swapaxes(y, 0, self.dir) + if self.axis > 0: + y = ncp.swapaxes(y, 0, self.axis) y = y.ravel() return y def _rmatvec(self, x): ncp = get_array_module(x) x = ncp.reshape(x, self.dimsd) - if self.dir > 0: # bring the dimension to symmetrize to first - x = ncp.swapaxes(x, self.dir, 0) + if self.axis > 0: # bring the dimension to symmetrize to first + x = ncp.swapaxes(x, self.axis, 0) y = x[self.nsym - 1 :].copy() y[1:] += x[self.nsym - 2 :: -1] - if self.dir > 0: - y = ncp.swapaxes(y, 0, self.dir) + if self.axis > 0: + y = ncp.swapaxes(y, 0, self.axis) y = y.ravel() return y diff --git a/pylops/signalprocessing/Convolve1D.py b/pylops/signalprocessing/Convolve1D.py index 5093e56e..ff30679a 100644 --- a/pylops/signalprocessing/Convolve1D.py +++ b/pylops/signalprocessing/Convolve1D.py @@ -1,4 +1,7 @@ +import warnings + import numpy as np +from numpy.core.multiarray import normalize_axis_index from pylops import LinearOperator from pylops.utils._internal import _value_or_list_like_to_array @@ -37,8 +40,7 @@ class Convolve1D(LinearOperator): r"""1D convolution operator. Apply one-dimensional convolution with a compact filter to model (and data) - along a specific direction of a multi-dimensional array depending on the - choice of ``dir``. + along an ``axis`` of a multi-dimensional array. Parameters ---------- @@ -48,8 +50,10 @@ class Convolve1D(LinearOperator): 1d compact filter to be convolved to input signal offset : :obj:`int` Index of the center of the compact filter - dir : :obj:`int`, optional - Direction along which convolution is applied + axis : :obj:`int`, optional + .. versionadded:: 2.0.0 + + Axis along which convolution is applied method : :obj:`str`, optional Method used to calculate the convolution (``direct``, ``fft``, or ``overlapadd``). Note that only ``direct`` and ``fft`` are allowed @@ -113,9 +117,9 @@ class Convolve1D(LinearOperator): """ - def __init__(self, dims, h, offset=0, dir=0, dtype="float64", method=None): + def __init__(self, dims, h, offset=0, axis=-1, dtype="float64", method=None): self.dims = _value_or_list_like_to_array(dims) - self.dir = dir + self.axis = axis if offset > len(h) - 1: raise ValueError("offset must be smaller than len(h) - 1") @@ -139,7 +143,7 @@ def __init__(self, dims, h, offset=0, dir=0, dtype="float64", method=None): # add dimensions to filter to match dimensions of model and data hdims = np.ones(len(self.dims), dtype=int) - hdims[self.dir] = len(self.h) + hdims[self.axis] = len(self.h) self.h = self.h.reshape(hdims) self.hstar = self.hstar.reshape(hdims) diff --git a/pylops/signalprocessing/Convolve2D.py b/pylops/signalprocessing/Convolve2D.py index 54a8e3a4..f2232748 100644 --- a/pylops/signalprocessing/Convolve2D.py +++ b/pylops/signalprocessing/Convolve2D.py @@ -1,13 +1,16 @@ +import warnings + +from numpy.core.multiarray import normalize_axis_index + from pylops.signalprocessing import ConvolveND -from pylops.utils._internal import _value_or_list_like_to_array -def Convolve2D(dims, h, offset=(0, 0), nodir=None, dtype="float64", method="fft"): +def Convolve2D(dims, h, offset=(0, 0), axes=(-2, -1), dtype="float64", method="fft"): r"""2D convolution operator. Apply two-dimensional convolution with a compact filter to model - (and data) along a pair of specific directions of a two or - three-dimensional array depending on the choice of ``nodir``. + (and data) along a pair of ``axes`` of a two or + three-dimensional array. Parameters ---------- @@ -16,10 +19,11 @@ def Convolve2D(dims, h, offset=(0, 0), nodir=None, dtype="float64", method="fft" h : :obj:`numpy.ndarray` 2d compact filter to be convolved to input signal offset : :obj:`tuple`, optional - Indeces of the center of the compact filter - nodir : :obj:`int`, optional - Direction along which convolution is NOT applied - (set to ``None`` for 2d arrays) + Indices of the center of the compact filter + axes : :obj:`int`, optional + .. versionadded:: 2.0.0 + + Axes along which convolution is applied dtype : :obj:`str`, optional Type of elements in input array. method : :obj:`str`, optional @@ -72,14 +76,5 @@ def Convolve2D(dims, h, offset=(0, 0), nodir=None, dtype="float64", method="fft" """ if h.ndim != 2: raise ValueError("h must be 2-dimensional") - if nodir is None: - dirs = (0, 1) - elif nodir == 0: - dirs = (1, 2) - elif nodir == 1: - dirs = (0, 2) - else: - dirs = (0, 1) - - cop = ConvolveND(dims, h, offset=offset, dirs=dirs, method=method, dtype=dtype) + cop = ConvolveND(dims, h, offset=offset, axes=axes, method=method, dtype=dtype) return cop diff --git a/pylops/signalprocessing/ConvolveND.py b/pylops/signalprocessing/ConvolveND.py index 5d999377..b99b8bdd 100644 --- a/pylops/signalprocessing/ConvolveND.py +++ b/pylops/signalprocessing/ConvolveND.py @@ -1,4 +1,7 @@ +import warnings + import numpy as np +from numpy.core.multiarray import normalize_axis_index from pylops import LinearOperator from pylops.utils._internal import _value_or_list_like_to_array @@ -14,8 +17,7 @@ class ConvolveND(LinearOperator): r"""ND convolution operator. Apply n-dimensional convolution with a compact filter to model - (and data) along a set of directions ``dirs`` of a n-dimensional - array. + (and data) along the ``axes`` of a n-dimensional array. Parameters ---------- @@ -25,9 +27,10 @@ class ConvolveND(LinearOperator): nd compact filter to be convolved to input signal offset : :obj:`tuple`, optional Indices of the center of the compact filter - dirs : :obj:`tuple`, optional - Directions along which convolution is applied - (set to ``None`` for filter of same dimension as input vector) + axes : :obj:`int`, optional + .. versionadded:: 2.0.0 + + Axes along which convolution is applied method : :obj:`str`, optional Method used to calculate the convolution (``direct`` or ``fft``). dtype : :obj:`str`, optional @@ -51,12 +54,24 @@ class ConvolveND(LinearOperator): """ - def __init__(self, dims, h, offset=None, dirs=None, method="fft", dtype="float64"): + def __init__( + self, + dims, + h, + offset=None, + axes=(-2, -1), + method="fft", + dtype="float64", + ): ncp = get_array_module(h) self.dims = _value_or_list_like_to_array(dims) + self.axes = ( + np.arange(len(self.dims)) + if axes is None + else np.array([normalize_axis_index(ax, len(self.dims)) for ax in axes]) + ) self.h = h self.nh = np.array(self.h.shape) - self.dirs = np.arange(len(self.dims)) if dirs is None else np.array(dirs) # padding if offset is None: @@ -82,8 +97,8 @@ def __init__(self, dims, h, offset=None, dirs=None, method="fft", dtype="float64 # find out which directions are used for convolution and define offsets if len(self.dims) != len(self.nh): dimsh = np.ones(len(self.dims), dtype=int) - for idir, dir in enumerate(self.dirs): - dimsh[dir] = self.nh[idir] + for iax, ax in enumerate(self.axes): + dimsh[ax] = self.nh[iax] self.h = self.h.reshape(dimsh) # convolve and correlate functions diff --git a/pylops/signalprocessing/DWT.py b/pylops/signalprocessing/DWT.py index 0698fdcd..3253ce5e 100644 --- a/pylops/signalprocessing/DWT.py +++ b/pylops/signalprocessing/DWT.py @@ -1,4 +1,5 @@ import logging +import warnings from math import ceil, log import numpy as np @@ -42,7 +43,7 @@ def _adjointwavelet(wavelet): class DWT(LinearOperator): """One dimensional Wavelet operator. - Apply 1D-Wavelet Transform along a specific direction ``dir`` of a + Apply 1D-Wavelet Transform along an ``axis`` of a multi-dimensional array of size ``dims``. Note that the Wavelet operator is an overload of the ``pywt`` @@ -54,8 +55,10 @@ class DWT(LinearOperator): ---------- dims : :obj:`int` or :obj:`tuple` Number of samples for each dimension - dir : :obj:`int`, optional - Direction along which DWT is applied. + axis : :obj:`int`, optional + .. versionadded:: 2.0.0 + + Axis along which DWT is applied wavelet : :obj:`str`, optional Name of wavelet type. Use :func:`pywt.wavelist(kind='discrete')` for a list of @@ -93,7 +96,7 @@ class DWT(LinearOperator): """ - def __init__(self, dims, dir=0, wavelet="haar", level=1, dtype="float64"): + def __init__(self, dims, axis=-1, wavelet="haar", level=1, dtype="float64"): if pywt is None: raise ModuleNotFoundError(pywt_message) _checkwavelet(wavelet) @@ -102,14 +105,14 @@ def __init__(self, dims, dir=0, wavelet="haar", level=1, dtype="float64"): dims = (dims,) # define padding for length to be power of 2 - ndimpow2 = max(2 ** ceil(log(dims[dir], 2)), 2 ** level) + ndimpow2 = max(2 ** ceil(log(dims[axis], 2)), 2 ** level) pad = [(0, 0)] * len(dims) - pad[dir] = (0, ndimpow2 - dims[dir]) + pad[axis] = (0, ndimpow2 - dims[axis]) self.pad = Pad(dims, pad) self.dims = dims - self.dir = dir + self.axis = axis self.dimsd = list(dims) - self.dimsd[self.dir] = ndimpow2 + self.dimsd[self.axis] = ndimpow2 # apply transform to find out slices _, self.sl = pywt.coeffs_to_array( @@ -118,9 +121,9 @@ def __init__(self, dims, dir=0, wavelet="haar", level=1, dtype="float64"): wavelet=wavelet, level=level, mode="periodization", - axes=(self.dir,), + axes=(self.axis,), ), - axes=(self.dir,), + axes=(self.axis,), ) self.wavelet = wavelet @@ -141,9 +144,9 @@ def _matvec(self, x): wavelet=self.wavelet, level=self.level, mode="periodization", - axes=(self.dir,), + axes=(self.axis,), ), - axes=(self.dir,), + axes=(self.axis,), )[0] return y.ravel() @@ -152,7 +155,7 @@ def _rmatvec(self, x): x = np.reshape(x, self.dimsd) x = pywt.array_to_coeffs(x, self.sl, output_format="wavedecn") y = pywt.waverecn( - x, wavelet=self.waveletadj, mode="periodization", axes=(self.dir,) + x, wavelet=self.waveletadj, mode="periodization", axes=(self.axis,) ) y = self.pad.rmatvec(y.ravel()) return y diff --git a/pylops/signalprocessing/DWT2D.py b/pylops/signalprocessing/DWT2D.py index f310a529..c6996ca0 100644 --- a/pylops/signalprocessing/DWT2D.py +++ b/pylops/signalprocessing/DWT2D.py @@ -1,4 +1,5 @@ import logging +import warnings from math import ceil, log import numpy as np @@ -12,6 +13,14 @@ import pywt except ModuleNotFoundError: pywt = None + pywt_message = ( + "Pywt package not installed. " + 'Run "pip install PyWavelets" or ' + 'conda install pywavelets".' + ) +except Exception as e: + pywt = None + pywt_message = "Failed to import pywt (error:%s)." % e logging.basicConfig(format="%(levelname)s: %(message)s", level=logging.WARNING) @@ -19,7 +28,7 @@ class DWT2D(LinearOperator): """Two dimensional Wavelet operator. - Apply 2D-Wavelet Transform along two directions ``dirs`` of a + Apply 2D-Wavelet Transform along two ``axes`` of a multi-dimensional array of size ``dims``. Note that the Wavelet operator is an overload of the ``pywt`` @@ -31,8 +40,10 @@ class DWT2D(LinearOperator): ---------- dims : :obj:`tuple` Number of samples for each dimension - dirs : :obj:`tuple`, optional - Direction along which DWT2D is applied. + axes : :obj:`int`, optional + .. versionadded:: 2.0.0 + + Axis along which DWT2D is applied wavelet : :obj:`str`, optional Name of wavelet type. Use :func:`pywt.wavelist(kind='discrete')` for a list of available wavelets. @@ -64,27 +75,22 @@ class DWT2D(LinearOperator): """ - def __init__(self, dims, dirs=(0, 1), wavelet="haar", level=1, dtype="float64"): + def __init__(self, dims, axes=(-2, -1), wavelet="haar", level=1, dtype="float64"): if pywt is None: - raise ModuleNotFoundError( - "The wavelet operator requires " - "the pywt package t be installed. " - 'Run "pip install PyWavelets" or ' - '"conda install pywavelets".' - ) + raise ModuleNotFoundError(pywt_message) _checkwavelet(wavelet) # define padding for length to be power of 2 - ndimpow2 = [max(2 ** ceil(log(dims[dir], 2)), 2 ** level) for dir in dirs] + ndimpow2 = [max(2 ** ceil(log(dims[ax], 2)), 2 ** level) for ax in axes] pad = [(0, 0)] * len(dims) - for i, dir in enumerate(dirs): - pad[dir] = (0, ndimpow2[i] - dims[dir]) + for i, ax in enumerate(axes): + pad[ax] = (0, ndimpow2[i] - dims[ax]) self.pad = Pad(dims, pad) self.dims = dims - self.dirs = dirs + self.axes = axes self.dimsd = list(dims) - for i, dir in enumerate(dirs): - self.dimsd[dir] = ndimpow2[i] + for i, ax in enumerate(axes): + self.dimsd[ax] = ndimpow2[i] # apply transform once again to find out slices _, self.sl = pywt.coeffs_to_array( @@ -93,9 +99,9 @@ def __init__(self, dims, dirs=(0, 1), wavelet="haar", level=1, dtype="float64"): wavelet=wavelet, level=level, mode="periodization", - axes=self.dirs, + axes=self.axes, ), - axes=self.dirs, + axes=self.axes, ) self.wavelet = wavelet self.waveletadj = _adjointwavelet(wavelet) @@ -113,9 +119,9 @@ def _matvec(self, x): wavelet=self.wavelet, level=self.level, mode="periodization", - axes=self.dirs, + axes=self.axes, ), - axes=(self.dirs), + axes=(self.axes), )[0] return y.ravel() @@ -123,7 +129,7 @@ def _rmatvec(self, x): x = np.reshape(x, self.dimsd) x = pywt.array_to_coeffs(x, self.sl, output_format="wavedec2") y = pywt.waverec2( - x, wavelet=self.waveletadj, mode="periodization", axes=self.dirs + x, wavelet=self.waveletadj, mode="periodization", axes=self.axes ) y = self.pad.rmatvec(y.ravel()) return y diff --git a/pylops/signalprocessing/FFT.py b/pylops/signalprocessing/FFT.py index cdd2e66e..ddc13441 100644 --- a/pylops/signalprocessing/FFT.py +++ b/pylops/signalprocessing/FFT.py @@ -28,7 +28,7 @@ class _FFT_numpy(_BaseFFT): def __init__( self, dims, - dir=0, + axis=-1, nfft=None, sampling=1.0, norm="ortho", @@ -39,7 +39,7 @@ def __init__( ): super().__init__( dims=dims, - dir=dir, + axis=axis, nfft=nfft, sampling=sampling, norm=norm, @@ -64,21 +64,21 @@ def __init__( def _matvec(self, x): x = np.reshape(x, self.dims) if self.ifftshift_before: - x = np.fft.ifftshift(x, axes=self.dir) + x = np.fft.ifftshift(x, axes=self.axis) if not self.clinear: x = np.real(x) if self.real: - y = np.fft.rfft(x, n=self.nfft, axis=self.dir, **self._norm_kwargs) + y = np.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.dir) + y = np.swapaxes(y, -1, self.axis) y[..., 1 : 1 + (self.nfft - 1) // 2] *= np.sqrt(2) - y = np.swapaxes(y, self.dir, -1) + y = np.swapaxes(y, self.axis, -1) else: - y = np.fft.fft(x, n=self.nfft, axis=self.dir, **self._norm_kwargs) + y = np.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 = np.fft.fftshift(y, axes=self.dir) + y = np.fft.fftshift(y, axes=self.axis) y = y.ravel() y = y.astype(self.cdtype) return y @@ -86,28 +86,28 @@ def _matvec(self, x): def _rmatvec(self, x): x = np.reshape(x, self.dims_fft) if self.fftshift_after: - x = np.fft.ifftshift(x, axes=self.dir) + x = np.fft.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.dir) + x = np.swapaxes(x, -1, self.axis) x[..., 1 : 1 + (self.nfft - 1) // 2] /= np.sqrt(2) - x = np.swapaxes(x, self.dir, -1) - y = np.fft.irfft(x, n=self.nfft, axis=self.dir, **self._norm_kwargs) + x = np.swapaxes(x, self.axis, -1) + y = np.fft.irfft(x, n=self.nfft, axis=self.axis, **self._norm_kwargs) else: - y = np.fft.ifft(x, n=self.nfft, axis=self.dir, **self._norm_kwargs) + y = np.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.dir]: - y = np.take(y, range(0, self.dims[self.dir]), axis=self.dir) - elif self.nfft < self.dims[self.dir]: + 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 = np.fft.fftshift(y, axes=self.dir) + y = np.fft.fftshift(y, axes=self.axis) y = y.ravel() y = y.astype(self.rdtype) return y @@ -124,7 +124,7 @@ class _FFT_scipy(_BaseFFT): def __init__( self, dims, - dir=0, + axis=-1, nfft=None, sampling=1.0, norm="ortho", @@ -135,7 +135,7 @@ def __init__( ): super().__init__( dims=dims, - dir=dir, + axis=axis, nfft=nfft, sampling=sampling, norm=norm, @@ -156,49 +156,49 @@ def __init__( def _matvec(self, x): x = np.reshape(x, self.dims) if self.ifftshift_before: - x = scipy.fft.ifftshift(x, axes=self.dir) + x = scipy.fft.ifftshift(x, axes=self.axis) if not self.clinear: x = np.real(x) if self.real: - y = scipy.fft.rfft(x, n=self.nfft, axis=self.dir, **self._norm_kwargs) + y = scipy.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.dir) + y = np.swapaxes(y, -1, self.axis) y[..., 1 : 1 + (self.nfft - 1) // 2] *= np.sqrt(2) - y = np.swapaxes(y, self.dir, -1) + y = np.swapaxes(y, self.axis, -1) else: - y = scipy.fft.fft(x, n=self.nfft, axis=self.dir, **self._norm_kwargs) + y = scipy.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 = scipy.fft.fftshift(y, axes=self.dir) + y = scipy.fft.fftshift(y, axes=self.axis) y = y.ravel() return y def _rmatvec(self, x): x = np.reshape(x, self.dims_fft) if self.fftshift_after: - x = scipy.fft.ifftshift(x, axes=self.dir) + x = scipy.fft.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.dir) + x = np.swapaxes(x, -1, self.axis) x[..., 1 : 1 + (self.nfft - 1) // 2] /= np.sqrt(2) - x = np.swapaxes(x, self.dir, -1) - y = scipy.fft.irfft(x, n=self.nfft, axis=self.dir, **self._norm_kwargs) + x = np.swapaxes(x, self.axis, -1) + y = scipy.fft.irfft(x, n=self.nfft, axis=self.axis, **self._norm_kwargs) else: - y = scipy.fft.ifft(x, n=self.nfft, axis=self.dir, **self._norm_kwargs) + y = scipy.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.dir]: - y = np.take(y, range(0, self.dims[self.dir]), axis=self.dir) - elif self.nfft < self.dims[self.dir]: + 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 = scipy.fft.fftshift(y, axes=self.dir) + y = scipy.fft.fftshift(y, axes=self.axis) y = y.ravel() return y @@ -214,7 +214,7 @@ class _FFT_fftw(_BaseFFT): def __init__( self, dims, - dir=0, + axis=-1, nfft=None, sampling=1.0, norm="ortho", @@ -241,7 +241,7 @@ def __init__( super().__init__( dims=dims, - dir=dir, + axis=axis, nfft=nfft, sampling=sampling, norm=norm, @@ -256,21 +256,21 @@ def __init__( ) self.dims_t = self.dims.copy() - self.dims_t[self.dir] = self.nfft + self.dims_t[self.axis] = self.nfft # define padding(fftw requires the user to provide padded input signal) self.pad = np.zeros((self.ndim, 2), dtype=int) if self.real: if self.nfft % 2: - self.pad[self.dir, 1] = ( - 2 * (self.dims_fft[self.dir] - 1) + 1 - self.dims[self.dir] + self.pad[self.axis, 1] = ( + 2 * (self.dims_fft[self.axis] - 1) + 1 - self.dims[self.axis] ) else: - self.pad[self.dir, 1] = ( - 2 * (self.dims_fft[self.dir] - 1) - self.dims[self.dir] + self.pad[self.axis, 1] = ( + 2 * (self.dims_fft[self.axis] - 1) - self.dims[self.axis] ) else: - self.pad[self.dir, 1] = self.dims_fft[self.dir] - self.dims[self.dir] + self.pad[self.axis, 1] = self.dims_fft[self.axis] - self.dims[self.axis] self.dopad = True if np.sum(self.pad) > 0 else False # create empty arrays and plans for fft/ifft @@ -292,22 +292,22 @@ def __init__( self._scale = 1.0 / self.nfft self.fftplan = pyfftw.FFTW( - self.x, self.y, axes=(self.dir,), direction="FFTW_FORWARD", **kwargs_fftw + self.x, self.y, axes=(self.axis,), direction="FFTW_FORWARD", **kwargs_fftw ) self.ifftplan = pyfftw.FFTW( - self.y, self.x, axes=(self.dir,), direction="FFTW_BACKWARD", **kwargs_fftw + self.y, self.x, axes=(self.axis,), direction="FFTW_BACKWARD", **kwargs_fftw ) def _matvec(self, x): x = np.reshape(x, self.dims) if self.ifftshift_before: - x = np.fft.ifftshift(x, axes=self.dir) + x = np.fft.ifftshift(x, axes=self.axis) if not self.clinear: x = np.real(x) if self.dopad: x = np.pad(x, self.pad, "constant", constant_values=0) elif self.doifftpad: - x = np.take(x, range(0, self.nfft), axis=self.dir) + x = np.take(x, range(0, self.nfft), axis=self.axis) # self.fftplan() always uses byte-alligned self.x as input array and # returns self.y as output array. As such, self.x must be copied so as @@ -319,17 +319,17 @@ def _matvec(self, x): if self.real: # Apply scaling to obtain a correct adjoint for this operator - y = np.swapaxes(y, -1, self.dir) + y = np.swapaxes(y, -1, self.axis) y[..., 1 : 1 + (self.nfft - 1) // 2] *= np.sqrt(2) - y = np.swapaxes(y, self.dir, -1) + y = np.swapaxes(y, self.axis, -1) if self.fftshift_after: - y = np.fft.fftshift(y, axes=self.dir) + y = np.fft.fftshift(y, axes=self.axis) return y.ravel() def _rmatvec(self, x): x = np.reshape(x, self.dims_fft) if self.fftshift_after: - x = np.fft.ifftshift(x, axes=self.dir) + x = np.fft.ifftshift(x, axes=self.axis) # self.ifftplan() always uses byte-alligned self.y as input array. # We copy here so we don't need to copy again in the case of `real=True`, @@ -339,9 +339,9 @@ def _rmatvec(self, x): if self.real: # Apply scaling to obtain a correct adjoint for this operator - x = np.swapaxes(x, -1, self.dir) + x = np.swapaxes(x, -1, self.axis) x[..., 1 : 1 + (self.nfft - 1) // 2] /= np.sqrt(2) - x = np.swapaxes(x, self.dir, -1) + x = np.swapaxes(x, self.axis, -1) # self.ifftplan() always returns self.x, which must be copied so as not # to be overwritten on a subsequent call to _rmatvec. @@ -351,13 +351,13 @@ def _rmatvec(self, x): elif self.norm is _FFTNorms.NONE: y *= self._scale - if self.nfft > self.dims[self.dir]: - y = np.take(y, range(0, self.dims[self.dir]), axis=self.dir) - elif self.nfft < self.dims[self.dir]: + 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 self.ifftshift_before: - y = np.fft.fftshift(y, axes=self.dir) + y = np.fft.fftshift(y, axes=self.axis) if not self.clinear: y = np.real(y) return y.ravel() @@ -370,7 +370,7 @@ def __truediv__(self, y): def FFT( dims, - dir=0, + axis=-1, nfft=None, sampling=1.0, norm="ortho", @@ -384,7 +384,7 @@ def FFT( ): r"""One dimensional Fast-Fourier Transform. - Apply Fast-Fourier Transform (FFT) along a specific direction ``dir`` of a + Apply Fast-Fourier Transform (FFT) along an ``axis`` of a multi-dimensional array of size ``dim``. Using the default NumPy engine, the FFT operator is an overload to either the NumPy @@ -410,8 +410,10 @@ def FFT( ---------- dims : :obj:`tuple` Number of samples for each dimension - dir : :obj:`int`, optional - Direction along which FFT is applied. + axis : :obj:`int`, optional + .. versionadded:: 2.0.0 + + Axis along which FFT is applied nfft : :obj:`int`, optional Number of samples in Fourier Transform (same as input if ``nfft=None``) sampling : :obj:`float`, optional @@ -503,7 +505,7 @@ def FFT( Raises ------ ValueError - - If ``dims`` is provided and ``dir`` is bigger than ``len(dims)``. + - 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``. @@ -559,7 +561,7 @@ def FFT( if engine == "fftw" and pyfftw is not None: f = _FFT_fftw( dims, - dir=dir, + axis=axis, nfft=nfft, sampling=sampling, norm=norm, @@ -574,7 +576,7 @@ def FFT( logging.warning(pyfftw_message) f = _FFT_numpy( dims, - dir=dir, + axis=axis, nfft=nfft, sampling=sampling, norm=norm, @@ -586,7 +588,7 @@ def FFT( elif engine == "scipy": f = _FFT_scipy( dims, - dir=dir, + axis=axis, nfft=nfft, sampling=sampling, norm=norm, diff --git a/pylops/signalprocessing/FFT2D.py b/pylops/signalprocessing/FFT2D.py index b936abff..b18a9023 100644 --- a/pylops/signalprocessing/FFT2D.py +++ b/pylops/signalprocessing/FFT2D.py @@ -15,7 +15,7 @@ class _FFT2D_numpy(_BaseFFTND): def __init__( self, dims, - dirs=(0, 1), + axes=(-2, -1), nffts=None, sampling=1.0, norm="ortho", @@ -26,7 +26,7 @@ def __init__( ): super().__init__( dims=dims, - dirs=dirs, + axes=axes, nffts=nffts, sampling=sampling, norm=norm, @@ -43,7 +43,7 @@ def __init__( # checks if self.ndim < 2: raise ValueError("FFT2D requires at least two input dimensions") - if self.ndirs != 2: + if self.naxes != 2: raise ValueError("FFT2D must be applied along exactly two dimensions") self.f1, self.f2 = self.fs @@ -60,50 +60,50 @@ def __init__( def _matvec(self, x): x = np.reshape(x, self.dims) if self.ifftshift_before.any(): - x = np.fft.ifftshift(x, axes=self.dirs[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: - y = np.fft.rfft2(x, s=self.nffts, axes=self.dirs, **self._norm_kwargs) + y = np.fft.rfft2(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.dirs[-1]) + y = np.swapaxes(y, -1, self.axes[-1]) y[..., 1 : 1 + (self.nffts[-1] - 1) // 2] *= np.sqrt(2) - y = np.swapaxes(y, self.dirs[-1], -1) + y = np.swapaxes(y, self.axes[-1], -1) else: - y = np.fft.fft2(x, s=self.nffts, axes=self.dirs, **self._norm_kwargs) + y = np.fft.fft2(x, s=self.nffts, axes=self.axes, **self._norm_kwargs) if self.norm is _FFTNorms.ONE_OVER_N: y *= self._scale y = y.astype(self.cdtype) if self.fftshift_after.any(): - y = np.fft.fftshift(y, axes=self.dirs[self.fftshift_after]) + y = np.fft.fftshift(y, axes=self.axes[self.fftshift_after]) return y.ravel() def _rmatvec(self, x): x = np.reshape(x, self.dims_fft) if self.fftshift_after.any(): - x = np.fft.ifftshift(x, axes=self.dirs[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() - x = np.swapaxes(x, -1, self.dirs[-1]) + x = np.swapaxes(x, -1, self.axes[-1]) x[..., 1 : 1 + (self.nffts[-1] - 1) // 2] /= np.sqrt(2) - x = np.swapaxes(x, self.dirs[-1], -1) - y = np.fft.irfft2(x, s=self.nffts, axes=self.dirs, **self._norm_kwargs) + x = np.swapaxes(x, self.axes[-1], -1) + y = np.fft.irfft2(x, s=self.nffts, axes=self.axes, **self._norm_kwargs) else: - y = np.fft.ifft2(x, s=self.nffts, axes=self.dirs, **self._norm_kwargs) + y = np.fft.ifft2(x, s=self.nffts, axes=self.axes, **self._norm_kwargs) if self.norm is _FFTNorms.NONE: y *= self._scale - if self.nffts[0] > self.dims[self.dirs[0]]: - y = np.take(y, range(self.dims[self.dirs[0]]), axis=self.dirs[0]) - if self.nffts[1] > self.dims[self.dirs[1]]: - y = np.take(y, range(self.dims[self.dirs[1]]), axis=self.dirs[1]) + if self.nffts[0] > self.dims[self.axes[0]]: + y = np.take(y, range(self.dims[self.axes[0]]), axis=self.axes[0]) + if self.nffts[1] > self.dims[self.axes[1]]: + y = np.take(y, range(self.dims[self.axes[1]]), axis=self.axes[1]) if self.doifftpad: y = np.pad(y, self.ifftpad) if not self.clinear: y = np.real(y) y = y.astype(self.rdtype) if self.ifftshift_before.any(): - y = np.fft.fftshift(y, axes=self.dirs[self.ifftshift_before]) + y = np.fft.fftshift(y, axes=self.axes[self.ifftshift_before]) return y.ravel() def __truediv__(self, y): @@ -118,7 +118,7 @@ class _FFT2D_scipy(_BaseFFTND): def __init__( self, dims, - dirs=(0, 1), + axes=(-2, -1), nffts=None, sampling=1.0, norm="ortho", @@ -129,7 +129,7 @@ def __init__( ): super().__init__( dims=dims, - dirs=dirs, + axes=axes, nffts=nffts, sampling=sampling, norm=norm, @@ -142,7 +142,7 @@ def __init__( # checks if self.ndim < 2: raise ValueError("FFT2D requires at least two input dimensions") - if self.ndirs != 2: + if self.naxes != 2: raise ValueError("FFT2D must be applied along exactly two dimensions") self.f1, self.f2 = self.fs @@ -159,44 +159,44 @@ def __init__( def _matvec(self, x): x = np.reshape(x, self.dims) if self.ifftshift_before.any(): - x = scipy.fft.ifftshift(x, axes=self.dirs[self.ifftshift_before]) + x = scipy.fft.ifftshift(x, axes=self.axes[self.ifftshift_before]) if not self.clinear: x = np.real(x) if self.real: - y = scipy.fft.rfft2(x, s=self.nffts, axes=self.dirs, **self._norm_kwargs) + y = scipy.fft.rfft2(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.dirs[-1]) + y = np.swapaxes(y, -1, self.axes[-1]) y[..., 1 : 1 + (self.nffts[-1] - 1) // 2] *= np.sqrt(2) - y = np.swapaxes(y, self.dirs[-1], -1) + y = np.swapaxes(y, self.axes[-1], -1) else: - y = scipy.fft.fft2(x, s=self.nffts, axes=self.dirs, **self._norm_kwargs) + y = scipy.fft.fft2(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 = scipy.fft.fftshift(y, axes=self.dirs[self.fftshift_after]) + y = scipy.fft.fftshift(y, axes=self.axes[self.fftshift_after]) return y.ravel() def _rmatvec(self, x): x = np.reshape(x, self.dims_fft) if self.fftshift_after.any(): - x = scipy.fft.ifftshift(x, axes=self.dirs[self.fftshift_after]) + x = scipy.fft.ifftshift(x, axes=self.axes[self.fftshift_after]) if self.real: # Apply scaling to obtain a correct adjoint for this operator x = x.copy() - x = np.swapaxes(x, -1, self.dirs[-1]) + x = np.swapaxes(x, -1, self.axes[-1]) x[..., 1 : 1 + (self.nffts[-1] - 1) // 2] /= np.sqrt(2) - x = np.swapaxes(x, self.dirs[-1], -1) - y = scipy.fft.irfft2(x, s=self.nffts, axes=self.dirs, **self._norm_kwargs) + x = np.swapaxes(x, self.axes[-1], -1) + y = scipy.fft.irfft2(x, s=self.nffts, axes=self.axes, **self._norm_kwargs) else: - y = scipy.fft.ifft2(x, s=self.nffts, axes=self.dirs, **self._norm_kwargs) + y = scipy.fft.ifft2(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.dirs[0]]), axis=self.dirs[0]) - y = np.take(y, range(self.dims[self.dirs[1]]), axis=self.dirs[1]) + 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 = scipy.fft.fftshift(y, axes=self.dirs[self.ifftshift_before]) + y = scipy.fft.fftshift(y, axes=self.axes[self.ifftshift_before]) return y.ravel() def __truediv__(self, y): @@ -207,7 +207,7 @@ def __truediv__(self, y): def FFT2D( dims, - dirs=(0, 1), + axes=(-2, -1), nffts=None, sampling=1.0, norm="ortho", @@ -219,8 +219,8 @@ def FFT2D( ): r"""Two dimensional Fast-Fourier Transform. - Apply two dimensional Fast-Fourier Transform (FFT) to any pair of axes of a - multi-dimensional array depending on the choice of ``dirs``. + Apply two dimensional Fast-Fourier Transform (FFT) to any pair of ``axes`` of a + multi-dimensional array. Using the default NumPy engine, the FFT operator is an overload to either the NumPy :py:func:`numpy.fft.fft2` (or :py:func:`numpy.fft.rfft2` for real models) in @@ -236,25 +236,27 @@ def FFT2D( the adjoint is multiplied by :math:`1 / \sqrt{2}` for the same frequencies. For a real valued input signal, it is advised to use the flag ``real=True`` - as it stores the values of the Fourier transform of the last direction at positive + as it stores the values of the Fourier transform of the last axis in ``axes`` at positive frequencies only as values at negative frequencies are simply their complex conjugates. Parameters ---------- dims : :obj:`tuple` Number of samples for each dimension - dirs : :obj:`tuple`, optional - Pair of directions along which FFT2D is applied + axes : :obj:`int`, optional + .. versionadded:: 2.0.0 + + Pair of axes along which FFT2D is applied nffts : :obj:`tuple` or :obj:`int`, optional - Number of samples in Fourier Transform for each direction. In case only one + Number of samples in Fourier Transform for each axis in ``axes``. In case only one dimension needs to be specified, use ``None`` for the other dimension in the - tuple. The direction with None will use ``dims[dir]`` as ``nfft``. When - supplying a tuple, the order must agree with that of ``dirs``. When a single - value is passed, it will be used for both directions. As such the default is - equivalent to ``nffts=(None, None)``. + tuple. An axis with ``None`` will use ``dims[axis]`` as ``nfft``. + When supplying a tuple, the length must be 2. + When a single value is passed, it will be used for both + axes. As such the default is equivalent to ``nffts=(None, None)``. sampling : :obj:`tuple` or :obj:`float`, optional - Sampling steps for each direction. When supplied a single value, it is used - for both directions. Unlike ``nffts``, any ``None`` will not be converted to the + Sampling steps for each axis in ``axes``. When supplied a single value, it is used + for both axes. Unlike ``nffts``, any ``None`` will not be converted to the default value. norm : `{"ortho", "none", "1/n"}`, optional .. versionadded:: 1.17.0 @@ -276,7 +278,7 @@ def FFT2D( (``False``). Used to enforce that the output of adjoint of a real model is real. Note that the real FFT is applied only to the first dimension to which the FFT2D operator is applied (last element of - ``dirs``) + ``axes``) ifftshift_before : :obj:`tuple` or :obj:`bool`, optional Apply ifftshift (``True``) or not (``False``) to model vector (before FFT). Consider using this option when the model vector's respective axis is symmetric @@ -284,7 +286,7 @@ def FFT2D( coincide with the zero index sample. With such an arrangement, FFT will not introduce a sample-dependent phase-shift when compared to the continuous Fourier Transform. - When passing a single value, the shift will the same for every direction. Pass + When passing a single value, the shift will the same for every axis in ``axes``. Pass a tuple to specify which dimensions are shifted. fftshift_after : :obj:`tuple` or :obj:`bool`, optional Apply fftshift (``True``) or not (``False``) to data vector (after FFT). @@ -292,7 +294,7 @@ def FFT2D( naturally, from negative to positive. When not applying fftshift after FFT, frequencies are arranged from zero to largest positive, and then from negative Nyquist to the frequency bin before zero. - When passing a single value, the shift will the same for every direction. Pass + When passing a single value, the shift will the same for every axis in ``axes``. Pass a tuple to specify which dimensions are shifted. engine : :obj:`str`, optional .. versionadded:: 1.17.0 @@ -315,9 +317,9 @@ def FFT2D( For example, ``y_reshaped = (Op * x.ravel()).reshape(Op.dims_fft)``. f1 : :obj:`numpy.ndarray` - Discrete Fourier Transform sample frequencies along ``dir[0]`` + Discrete Fourier Transform sample frequencies along ``axes[0]`` f2 : :obj:`numpy.ndarray` - Discrete Fourier Transform sample frequencies along ``dir[1]`` + Discrete Fourier Transform sample frequencies along ``axes[1]`` real : :obj:`bool` When ``True``, uses ``rfft2``/``irfft2`` rdtype : :obj:`bool` @@ -339,7 +341,7 @@ def FFT2D( ------ ValueError - If ``dims`` has less than two elements. - - If ``dirs`` does not have exactly two elements. + - If ``axes`` does not have exactly two elements. - If ``nffts`` or ``sampling`` are not either a single value or a tuple with two elements. - If ``norm`` is not one of "ortho", "none", or "1/n". @@ -380,7 +382,7 @@ def FFT2D( if engine == "numpy": f = _FFT2D_numpy( dims=dims, - dirs=dirs, + axes=axes, nffts=nffts, sampling=sampling, norm=norm, @@ -392,7 +394,7 @@ def FFT2D( elif engine == "scipy": f = _FFT2D_scipy( dims=dims, - dirs=dirs, + axes=axes, nffts=nffts, sampling=sampling, norm=norm, diff --git a/pylops/signalprocessing/FFTND.py b/pylops/signalprocessing/FFTND.py index 87f1fe42..349bbb4d 100644 --- a/pylops/signalprocessing/FFTND.py +++ b/pylops/signalprocessing/FFTND.py @@ -15,7 +15,7 @@ class _FFTND_numpy(_BaseFFTND): def __init__( self, dims, - dirs=(0, 1, 2), + axes=(-3, -2, -1), nffts=None, sampling=1.0, norm="ortho", @@ -26,7 +26,7 @@ def __init__( ): super().__init__( dims=dims, - dirs=dirs, + axes=axes, nffts=nffts, sampling=sampling, norm=norm, @@ -51,49 +51,49 @@ def __init__( def _matvec(self, x): x = np.reshape(x, self.dims) if self.ifftshift_before.any(): - x = np.fft.ifftshift(x, axes=self.dirs[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: - y = np.fft.rfftn(x, s=self.nffts, axes=self.dirs, **self._norm_kwargs) + y = np.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.dirs[-1]) + y = np.swapaxes(y, -1, self.axes[-1]) y[..., 1 : 1 + (self.nffts[-1] - 1) // 2] *= np.sqrt(2) - y = np.swapaxes(y, self.dirs[-1], -1) + y = np.swapaxes(y, self.axes[-1], -1) else: - y = np.fft.fftn(x, s=self.nffts, axes=self.dirs, **self._norm_kwargs) + y = np.fft.fftn(x, s=self.nffts, axes=self.axes, **self._norm_kwargs) if self.norm is _FFTNorms.ONE_OVER_N: y *= self._scale y = y.astype(self.cdtype) if self.fftshift_after.any(): - y = np.fft.fftshift(y, axes=self.dirs[self.fftshift_after]) + y = np.fft.fftshift(y, axes=self.axes[self.fftshift_after]) return y.ravel() def _rmatvec(self, x): x = np.reshape(x, self.dims_fft) if self.fftshift_after.any(): - x = np.fft.ifftshift(x, axes=self.dirs[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() - x = np.swapaxes(x, -1, self.dirs[-1]) + x = np.swapaxes(x, -1, self.axes[-1]) x[..., 1 : 1 + (self.nffts[-1] - 1) // 2] /= np.sqrt(2) - x = np.swapaxes(x, self.dirs[-1], -1) - y = np.fft.irfftn(x, s=self.nffts, axes=self.dirs, **self._norm_kwargs) + x = np.swapaxes(x, self.axes[-1], -1) + y = np.fft.irfftn(x, s=self.nffts, axes=self.axes, **self._norm_kwargs) else: - y = np.fft.ifftn(x, s=self.nffts, axes=self.dirs, **self._norm_kwargs) + y = np.fft.ifftn(x, s=self.nffts, axes=self.axes, **self._norm_kwargs) if self.norm is _FFTNorms.NONE: y *= self._scale - for direction, nfft in zip(self.dirs, self.nffts): - if nfft > self.dims[direction]: - y = np.take(y, range(self.dims[direction]), axis=direction) + 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) y = y.astype(self.rdtype) if self.ifftshift_before.any(): - y = np.fft.fftshift(y, axes=self.dirs[self.ifftshift_before]) + y = np.fft.fftshift(y, axes=self.axes[self.ifftshift_before]) return y.ravel() def __truediv__(self, y): @@ -108,7 +108,7 @@ class _FFTND_scipy(_BaseFFTND): def __init__( self, dims, - dirs=(0, 1, 2), + axes=(-3, -2, -1), nffts=None, sampling=1.0, norm="ortho", @@ -119,7 +119,7 @@ def __init__( ): super().__init__( dims=dims, - dirs=dirs, + axes=axes, nffts=nffts, sampling=sampling, norm=norm, @@ -140,47 +140,47 @@ def __init__( def _matvec(self, x): x = np.reshape(x, self.dims) if self.ifftshift_before.any(): - x = scipy.fft.ifftshift(x, axes=self.dirs[self.ifftshift_before]) + x = scipy.fft.ifftshift(x, axes=self.axes[self.ifftshift_before]) if not self.clinear: x = np.real(x) if self.real: - y = scipy.fft.rfftn(x, s=self.nffts, axes=self.dirs, **self._norm_kwargs) + y = scipy.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.dirs[-1]) + y = np.swapaxes(y, -1, self.axes[-1]) y[..., 1 : 1 + (self.nffts[-1] - 1) // 2] *= np.sqrt(2) - y = np.swapaxes(y, self.dirs[-1], -1) + y = np.swapaxes(y, self.axes[-1], -1) else: - y = scipy.fft.fftn(x, s=self.nffts, axes=self.dirs, **self._norm_kwargs) + y = scipy.fft.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 = scipy.fft.fftshift(y, axes=self.dirs[self.fftshift_after]) + y = scipy.fft.fftshift(y, axes=self.axes[self.fftshift_after]) return y.ravel() def _rmatvec(self, x): x = np.reshape(x, self.dims_fft) if self.fftshift_after.any(): - x = scipy.fft.ifftshift(x, axes=self.dirs[self.fftshift_after]) + x = scipy.fft.ifftshift(x, axes=self.axes[self.fftshift_after]) if self.real: # Apply scaling to obtain a correct adjoint for this operator x = x.copy() - x = np.swapaxes(x, -1, self.dirs[-1]) + x = np.swapaxes(x, -1, self.axes[-1]) x[..., 1 : 1 + (self.nffts[-1] - 1) // 2] /= np.sqrt(2) - x = np.swapaxes(x, self.dirs[-1], -1) - y = scipy.fft.irfftn(x, s=self.nffts, axes=self.dirs, **self._norm_kwargs) + x = np.swapaxes(x, self.axes[-1], -1) + y = scipy.fft.irfftn(x, s=self.nffts, axes=self.axes, **self._norm_kwargs) else: - y = scipy.fft.ifftn(x, s=self.nffts, axes=self.dirs, **self._norm_kwargs) + y = scipy.fft.ifftn(x, s=self.nffts, axes=self.axes, **self._norm_kwargs) if self.norm is _FFTNorms.NONE: y *= self._scale - for direction, nfft in zip(self.dirs, self.nffts): - if nfft > self.dims[direction]: - y = np.take(y, range(self.dims[direction]), axis=direction) + 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 = scipy.fft.fftshift(y, axes=self.dirs[self.ifftshift_before]) + y = scipy.fft.fftshift(y, axes=self.axes[self.ifftshift_before]) return y.ravel() def __truediv__(self, y): @@ -191,7 +191,7 @@ def __truediv__(self, y): def FFTND( dims, - dirs=(0, 1, 2), + axes=(-3, -2, -1), nffts=None, sampling=1.0, norm="ortho", @@ -203,8 +203,8 @@ def FFTND( ): r"""N-dimensional Fast-Fourier Transform. - Apply N-dimensional Fast-Fourier Transform (FFT) to any n axes - of a multi-dimensional array depending on the choice of ``dirs``. + Apply N-dimensional Fast-Fourier Transform (FFT) to any n ``axes`` + of a multi-dimensional array. Using the default NumPy engine, the FFT operator is an overload to either the NumPy :py:func:`numpy.fft.fftn` (or :py:func:`numpy.fft.rfftn` for real models) in @@ -217,26 +217,28 @@ def FFTND( 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 - direction of ``dirs``, and the input of the adjoint is multiplied by + ``axes``, and the input of the adjoint is multiplied by :math:`1 / \sqrt{2}` for the same frequencies. For a real valued input signal, it is advised to use the flag ``real=True`` - as it stores the values of the Fourier transform of the last direction at positive + as it stores the values of the Fourier transform of the last axis in ``axes`` at positive frequencies only as values at negative frequencies are simply their complex conjugates. Parameters ---------- dims : :obj:`tuple` Number of samples for each dimension - dirs : :obj:`tuple` or :obj:`int`, optional - Direction(s) along which FFTND is applied + axes : :obj:`int`, optional + .. versionadded:: 2.0.0 + + Axes (or axis) along which FFTND is applied nffts : :obj:`tuple` or :obj:`int`, optional - Number of samples in Fourier Transform for each direction. In case only one + Number of samples in Fourier Transform for each axis in ``axes``. In case only one dimension needs to be specified, use ``None`` for the other dimension in the - tuple. The direction with None will use ``dims[dir]`` as ``nfft``. When - supplying a tuple, the order must agree with that of ``dirs``. When a single - value is passed, it will be used for both directions. As such the default is - equivalent to ``nffts=(None,..., None)``. + tuple. An axis with ``None`` will use ``dims[axis]`` as ``nfft``. + When supplying a tuple, the length must agree with that + of ``axes``. When a single value is passed, it will be used for all + ``axes`. As such the default is equivalent to ``nffts=(None, ..., None)``. sampling : :obj:`tuple` or :obj:`float`, optional Sampling steps for each direction. When supplied a single value, it is used for all directions. Unlike ``nffts``, any ``None`` will not be converted to the @@ -261,7 +263,7 @@ def FFTND( (``False``). Used to enforce that the output of adjoint of a real model is real. Note that the real FFT is applied only to the first dimension to which the FFTND operator is applied (last element of - ``dirs``) + ``axes``) ifftshift_before : :obj:`tuple` or :obj:`bool`, optional .. versionadded:: 1.17.0 @@ -305,7 +307,7 @@ def FFTND( For example, ``y_reshaped = (Op * x.ravel()).reshape(Op.dims_fft)``. fs : :obj:`tuple` Each element of the tuple corresponds to the Discrete Fourier Transform - sample frequencies along the respective direction given by ``dirs``. + sample frequencies along the respective direction given by ``axes``. real : :obj:`bool` When ``True``, uses ``rfftn``/``irfftn`` rdtype : :obj:`bool` @@ -332,7 +334,7 @@ def FFTND( ------ ValueError - If ``nffts`` or ``sampling`` are not either a single value or tuple with - the same dimension ``dirs``. + the same dimension ``axes``. - If ``norm`` is not one of "ortho", "none", or "1/n". NotImplementedError If ``engine`` is neither ``numpy``, nor ``scipy``. @@ -370,11 +372,10 @@ def FFTND( for real input signals. """ - if engine == "numpy": f = _FFTND_numpy( dims=dims, - dirs=dirs, + axes=axes, nffts=nffts, sampling=sampling, norm=norm, @@ -386,7 +387,7 @@ def FFTND( elif engine == "scipy": f = _FFTND_scipy( dims=dims, - dirs=dirs, + axes=axes, nffts=nffts, sampling=sampling, norm=norm, diff --git a/pylops/signalprocessing/Interp.py b/pylops/signalprocessing/Interp.py index 529ea7dd..3e3654d0 100644 --- a/pylops/signalprocessing/Interp.py +++ b/pylops/signalprocessing/Interp.py @@ -1,4 +1,5 @@ import logging +import warnings import numpy as np @@ -16,23 +17,23 @@ def _checkunique(iava): raise ValueError("Repeated values in iava array") -def _nearestinterp(dims, iava, dir=0, dtype="float64"): +def _nearestinterp(dims, iava, axis=-1, dtype="float64"): """Nearest neighbour interpolation.""" iava = np.round(iava).astype(int) _checkunique(iava) - return Restriction(dims, iava, dir=dir, dtype=dtype), iava + return Restriction(dims, iava, axis=axis, dtype=dtype), iava -def _linearinterp(dims, iava, dir=0, dtype="float64"): +def _linearinterp(dims, iava, axis=-1, dtype="float64"): """Linear interpolation.""" ncp = get_array_module(iava) if np.issubdtype(iava.dtype, np.integer): iava = iava.astype(np.float64) - lastsample = dims[dir] + lastsample = dims[axis] dimsd = list(dims) - dimsd[dir] = len(iava) + dimsd[axis] = len(iava) dimsd = tuple(dimsd) # ensure that samples are not beyond the last sample, in that case set to @@ -52,47 +53,47 @@ def _linearinterp(dims, iava, dir=0, dtype="float64"): weights = iava - iva_l # create operators - Op = Diagonal(1 - weights, dims=dimsd, dir=dir, dtype=dtype) * Restriction( - dims, iva_l, dir=dir, dtype=dtype - ) + Diagonal(weights, dims=dimsd, dir=dir, dtype=dtype) * Restriction( - dims, iva_r, dir=dir, dtype=dtype + Op = Diagonal(1 - weights, dims=dimsd, axis=axis, dtype=dtype) * Restriction( + dims, iva_l, axis=axis, dtype=dtype + ) + Diagonal(weights, dims=dimsd, axis=axis, dtype=dtype) * Restriction( + dims, iva_r, axis=axis, dtype=dtype ) return Op, iava -def _sincinterp(dims, iava, dir=0, dtype="float64"): +def _sincinterp(dims, iava, axis=0, dtype="float64"): """Sinc interpolation.""" ncp = get_array_module(iava) _checkunique(iava) # create sinc interpolation matrix - nreg = dims[dir] + nreg = dims[axis] ireg = ncp.arange(nreg) sinc = ncp.tile(iava[:, np.newaxis], (1, nreg)) - ncp.tile(ireg, (len(iava), 1)) sinc = ncp.sinc(sinc) # identify additional dimensions and create MatrixMult operator otherdims = ncp.array(dims) - otherdims = ncp.roll(otherdims, -dir) + otherdims = ncp.roll(otherdims, -axis) otherdims = otherdims[1:] Op = MatrixMult(sinc, dims=otherdims, dtype=dtype) - # create Transpose operator that brings dir to first dimension - if dir > 0: + # create Transpose operator that brings axis to first dimension + if axis > 0: axes = np.arange(len(dims), dtype=int) - axes = np.roll(axes, -dir) + axes = np.roll(axes, -axis) dimsd = list(dims) - dimsd[dir] = len(iava) + dimsd[axis] = len(iava) Top = Transpose(dims, axes=axes, dtype=dtype) T1op = Transpose(dimsd, axes=axes, dtype=dtype) Op = T1op.H * Op * Top return Op -def Interp(dims, iava, dir=0, kind="linear", dtype="float64"): +def Interp(dims, iava, axis=-1, kind="linear", dtype="float64"): r"""Interpolation operator. - Apply interpolation along direction ``dir`` + Apply interpolation along ``axis`` from regularly sampled input vector into fractionary positions ``iava`` using one of the following algorithms: @@ -123,8 +124,10 @@ def Interp(dims, iava, dir=0, kind="linear", dtype="float64"): Number of samples for each dimension iava : :obj:`list` or :obj:`numpy.ndarray` Floating indices of locations of available samples for interpolation. - dir : :obj:`int`, optional - Direction along which restriction is applied. + axis : :obj:`int`, optional + .. versionadded:: 2.0.0 + + Axis along which interpolation is applied. kind : :obj:`str`, optional Kind of interpolation (``nearest``, ``linear``, and ``sinc`` are currently supported) @@ -188,11 +191,11 @@ def Interp(dims, iava, dir=0, kind="linear", dtype="float64"): dims = _value_or_list_like_to_array(dims) if kind == "nearest": - interpop, iava = _nearestinterp(dims, iava, dir=dir, dtype=dtype) + interpop, iava = _nearestinterp(dims, iava, axis=axis, dtype=dtype) elif kind == "linear": - interpop, iava = _linearinterp(dims, iava, dir=dir, dtype=dtype) + interpop, iava = _linearinterp(dims, iava, axis=axis, dtype=dtype) elif kind == "sinc": - interpop = _sincinterp(dims, iava, dir=dir, dtype=dtype) + interpop = _sincinterp(dims, iava, axis=axis, dtype=dtype) else: raise NotImplementedError("kind is not correct...") return LinearOperator(interpop), iava diff --git a/pylops/signalprocessing/Patch2D.py b/pylops/signalprocessing/Patch2D.py index 824e7691..3d422313 100644 --- a/pylops/signalprocessing/Patch2D.py +++ b/pylops/signalprocessing/Patch2D.py @@ -167,7 +167,7 @@ def Patch2D(Op, dims, dimsd, nwin, nover, nop, tapertype="hanning", design=False hstack = HStack( [ Restriction( - (nwin[0], dimsd[1]), range(win_in, win_end), dir=1, dtype=Op.dtype + (nwin[0], dimsd[1]), range(win_in, win_end), axis=1, dtype=Op.dtype ).H for win_in, win_end in zip(dwin1_ins, dwin1_ends) ] @@ -176,7 +176,7 @@ def Patch2D(Op, dims, dimsd, nwin, nover, nop, tapertype="hanning", design=False combining1 = BlockDiag([hstack] * nwins0) combining0 = HStack( [ - Restriction(dimsd, range(win_in, win_end), dir=0, dtype=Op.dtype).H + Restriction(dimsd, range(win_in, win_end), axis=0, dtype=Op.dtype).H for win_in, win_end in zip(dwin0_ins, dwin0_ends) ] ) diff --git a/pylops/signalprocessing/Shift.py b/pylops/signalprocessing/Shift.py index 2aa1df76..e89e030e 100644 --- a/pylops/signalprocessing/Shift.py +++ b/pylops/signalprocessing/Shift.py @@ -1,3 +1,5 @@ +import warnings + import numpy as np from pylops.basicoperators import Diagonal @@ -7,7 +9,7 @@ def Shift( dims, shift, - dir=0, + axis=-1, nfft=None, sampling=1.0, real=False, @@ -17,8 +19,8 @@ def Shift( ): r"""Shift operator - Apply fractional shift in the frequency domain along a specific direction - ``dir`` of a multi-dimensional array of size ``dim``. + Apply fractional shift in the frequency domain along an ``axis`` + of a multi-dimensional array of size ``dims``. Parameters ---------- @@ -26,8 +28,10 @@ def Shift( Number of samples for each dimension shift : :obj:`float` Fractional shift to apply in the same unit as ``sampling``. - dir : :obj:`int`, optional - Direction along which FFT is applied. + axis : :obj:`int`, optional + .. versionadded:: 2.0.0 + + Axis along which shift is applied nfft : :obj:`int`, optional Number of samples in Fourier Transform (same as input if ``nfft=None``) sampling : :obj:`float`, optional @@ -56,7 +60,7 @@ def Shift( Raises ------ ValueError - If ``dims`` is provided and ``dir`` is bigger than ``len(dims)`` + If ``dims`` is provided and ``axis`` is bigger than ``len(dims)`` NotImplementedError If ``engine`` is neither ``numpy``, ``scipy``, nor ``fftw`` @@ -73,17 +77,23 @@ def Shift( chosen ``shift``. """ - # TODO: Use offer the same keywords as new FFT Fop = FFT( - dims, dir, nfft, sampling, real=real, engine=engine, dtype=dtype, **kwargs_fftw + dims, + axis=axis, + nfft=nfft, + sampling=sampling, + real=real, + engine=engine, + dtype=dtype, + **kwargs_fftw ) if isinstance(dims, int): dimsdiag = None else: dimsdiag = list(dims) - dimsdiag[dir] = len(Fop.f) + dimsdiag[axis] = len(Fop.f) shift = np.exp(-1j * 2 * np.pi * Fop.f * shift) - Sop = Diagonal(shift, dims=dimsdiag, dir=dir, dtype=Fop.cdtype) + Sop = Diagonal(shift, dims=dimsdiag, axis=axis, dtype=Fop.cdtype) Op = Fop.H * Sop * Fop # force dtype to that of input (FFT always upcasts it to complex) Op.dtype = dtype diff --git a/pylops/signalprocessing/Sliding2D.py b/pylops/signalprocessing/Sliding2D.py index 6d50ccb7..c55ff826 100644 --- a/pylops/signalprocessing/Sliding2D.py +++ b/pylops/signalprocessing/Sliding2D.py @@ -135,7 +135,7 @@ def Sliding2D(Op, dims, dimsd, nwin, nover, tapertype="hanning", design=False): combining = HStack( [ - Restriction(dimsd, range(win_in, win_end), dtype=Op.dtype).H + Restriction(dimsd, range(win_in, win_end), axis=0, dtype=Op.dtype).H for win_in, win_end in zip(dwin_ins, dwin_ends) ] ) diff --git a/pylops/signalprocessing/Sliding3D.py b/pylops/signalprocessing/Sliding3D.py index c7dab54d..470c76d5 100644 --- a/pylops/signalprocessing/Sliding3D.py +++ b/pylops/signalprocessing/Sliding3D.py @@ -133,7 +133,7 @@ def Sliding3D( Restriction( (nwin[0], dimsd[1], dimsd[2]), range(win_in, win_end), - dir=1, + axis=1, dtype=Op.dtype, ).H for win_in, win_end in zip(dwin1_ins, dwin1_ends) @@ -146,7 +146,7 @@ def Sliding3D( Restriction( dimsd, range(win_in, win_end), - dir=0, + axis=0, dtype=Op.dtype, ).H for win_in, win_end in zip(dwin0_ins, dwin0_ends) diff --git a/pylops/signalprocessing/_BaseFFTs.py b/pylops/signalprocessing/_BaseFFTs.py index 00724caa..95003ccb 100644 --- a/pylops/signalprocessing/_BaseFFTs.py +++ b/pylops/signalprocessing/_BaseFFTs.py @@ -24,7 +24,7 @@ class _BaseFFT(LinearOperator): def __init__( self, dims, - dir=0, + axis=-1, nfft=None, sampling=1.0, norm="ortho", @@ -38,12 +38,12 @@ def __init__( self.ndim = len(self.dims) - dirs = _value_or_list_like_to_array(dir) - _raise_on_wrong_dtype(dirs, np.integer, "dir") - self.dir = normalize_axis_index(dirs[0], self.ndim) + axes = _value_or_list_like_to_array(axis) + _raise_on_wrong_dtype(axes, np.integer, "axis") + self.axis = normalize_axis_index(axes[0], self.ndim) if nfft is None: - nfft = self.dims[self.dir] + nfft = self.dims[self.axis] else: nffts = _value_or_list_like_to_array(nfft) _raise_on_wrong_dtype(nffts, np.integer, "nfft") @@ -56,12 +56,12 @@ def __init__( # before applying fft, which is lost forever) and set a flag such that # a padding is applied after ifft self.doifftpad = False - if self.nfft < self.dims[dir]: + if self.nfft < self.dims[self.axis]: self.doifftpad = True self.ifftpad = [(0, 0)] * self.ndim - self.ifftpad[self.dir] = (0, self.dims[dir] - self.nfft) + self.ifftpad[self.axis] = (0, self.dims[self.axis] - self.nfft) warnings.warn( - f"nfft={self.nfft} has been selected to be smaller than the size of the original signal (n=self.dims[dir]). " + f"nfft={self.nfft} has been selected to be smaller than the size of the original signal (self.dims[axis]={self.dims[axis]}). " f"This is rarely intended behavior as the original signal will be truncated prior to applying fft, " f"if this is the required behaviour ignore this message." ) @@ -100,7 +100,7 @@ def __init__( self.f = np.fft.fftshift(self.f) self.dims_fft = self.dims.copy() - self.dims_fft[self.dir] = self.nfft // 2 + 1 if self.real else self.nfft + self.dims_fft[self.axis] = self.nfft // 2 + 1 if self.real else self.nfft self.shape = (int(np.prod(self.dims_fft)), int(np.prod(self.dims))) # Find types to enforce to forward and adjoint outputs. This is @@ -131,7 +131,7 @@ class _BaseFFTND(LinearOperator): def __init__( self, dims, - dirs=None, + axes=None, nffts=None, sampling=1.0, norm="ortho", @@ -145,51 +145,51 @@ def __init__( self.ndim = len(self.dims) - dirs = _value_or_list_like_to_array(dirs) - _raise_on_wrong_dtype(dirs, np.integer, "dirs") - self.dirs = np.array([normalize_axis_index(d, self.ndim) for d in dirs]) - self.ndirs = len(self.dirs) - if self.ndirs != len(np.unique(self.dirs)): + axes = _value_or_list_like_to_array(axes) + _raise_on_wrong_dtype(axes, np.integer, "axes") + self.axes = np.array([normalize_axis_index(d, self.ndim) for d in axes]) + self.naxes = len(self.axes) + if self.naxes != len(np.unique(self.axes)): warnings.warn( "At least one direction is repeated. This may cause unexpected results." ) - nffts = _value_or_list_like_to_array(nffts, repeat=self.ndirs) + nffts = _value_or_list_like_to_array(nffts, repeat=self.naxes) if len(nffts[np.equal(nffts, None)]) > 0: # Found None(s) in nffts nffts[np.equal(nffts, None)] = np.array( - [self.dims[d] for d, n in zip(dirs, nffts) if n is None] + [self.dims[d] for d, n in zip(axes, nffts) if n is None] ) nffts = nffts.astype(self.dims.dtype) self.nffts = nffts _raise_on_wrong_dtype(self.nffts, np.integer, "nffts") - sampling = _value_or_list_like_to_array(sampling, repeat=self.ndirs) + sampling = _value_or_list_like_to_array(sampling, repeat=self.naxes) if np.issubdtype(sampling.dtype, np.integer): # Promote to float64 if integer sampling = sampling.astype(np.float64) self.sampling = sampling _raise_on_wrong_dtype(self.sampling, np.floating, "sampling") self.ifftshift_before = _value_or_list_like_to_array( - ifftshift_before, repeat=self.ndirs + ifftshift_before, repeat=self.naxes ) _raise_on_wrong_dtype(self.ifftshift_before, bool, "ifftshift_before") self.fftshift_after = _value_or_list_like_to_array( - fftshift_after, repeat=self.ndirs + fftshift_after, repeat=self.naxes ) _raise_on_wrong_dtype(self.fftshift_after, bool, "fftshift_after") if ( - self.ndirs != len(self.nffts) - or self.ndirs != len(self.sampling) - or self.ndirs != len(self.ifftshift_before) - or self.ndirs != len(self.fftshift_after) + self.naxes != len(self.nffts) + or self.naxes != len(self.sampling) + or self.naxes != len(self.ifftshift_before) + or self.naxes != len(self.fftshift_after) ): raise ValueError( ( - "`dirs`, `nffts`, `sampling`, `ifftshift_before` and " + "`axes`, `nffts`, `sampling`, `ifftshift_before` and " "`fftshift_after` must the have same number of elements. Received " - f"{self.ndirs}, {len(self.nffts)}, {len(self.sampling)}, " + f"{self.naxes}, {len(self.nffts)}, {len(self.sampling)}, " f"{len(self.ifftshift_before)} and {len(self.fftshift_after)}, " "respectively." ) @@ -199,12 +199,12 @@ def __init__( # details nfftshort = [ nfft < self.dims[direction] - for direction, nfft in zip(self.dirs, self.nffts) + for direction, nfft in zip(self.axes, self.nffts) ] self.doifftpad = any(nfftshort) if self.doifftpad: self.ifftpad = [(0, 0)] * self.ndim - for idir, (direction, nfshort) in enumerate(zip(self.dirs, nfftshort)): + for idir, (direction, nfshort) in enumerate(zip(self.axes, nfftshort)): if nfshort: self.ifftpad[direction] = ( 0, @@ -257,9 +257,9 @@ def __init__( fs[-1] = np.fft.fftshift(fs[-1]) self.fs = tuple(fs) self.dims_fft = self.dims.copy() - self.dims_fft[self.dirs] = self.nffts + self.dims_fft[self.axes] = self.nffts if self.real: - self.dims_fft[self.dirs[-1]] = self.nffts[-1] // 2 + 1 + self.dims_fft[self.axes[-1]] = self.nffts[-1] // 2 + 1 self.shape = (int(np.prod(self.dims_fft)), int(np.prod(self.dims))) self.rdtype = get_real_dtype(dtype) if self.real else np.dtype(dtype) self.cdtype = get_complex_dtype(dtype) diff --git a/pylops/waveeqprocessing/lsm.py b/pylops/waveeqprocessing/lsm.py index cbbe00f9..99553ddf 100644 --- a/pylops/waveeqprocessing/lsm.py +++ b/pylops/waveeqprocessing/lsm.py @@ -266,7 +266,7 @@ def Demigration( dims=dims, dimsd=(ns * nr, nt), table=itrav, dtable=travd, engine="numba" ) - cop = Convolve1D((ns * nr, nt), h=wav, offset=wavcenter, dir=1) + cop = Convolve1D((ns * nr, nt), h=wav, offset=wavcenter, axis=1) demop = cop * sop else: raise NotImplementedError("method must be analytic or eikonal") diff --git a/pylops/waveeqprocessing/marchenko.py b/pylops/waveeqprocessing/marchenko.py index fdfb32ec..327e819b 100644 --- a/pylops/waveeqprocessing/marchenko.py +++ b/pylops/waveeqprocessing/marchenko.py @@ -402,7 +402,7 @@ def apply_onepoint( ) Rollop = Roll( (self.nt2, self.ns), - dir=0, + axis=0, shift=-1, dtype=self.dtype, ) @@ -611,7 +611,7 @@ def apply_multiplepoints( ) Rollop = Roll( (self.nt2, self.ns, nvs), - dir=0, + axis=0, shift=-1, dtype=self.dtype, ) diff --git a/pylops/waveeqprocessing/mdd.py b/pylops/waveeqprocessing/mdd.py index 87b66e2d..5c6ce650 100644 --- a/pylops/waveeqprocessing/mdd.py +++ b/pylops/waveeqprocessing/mdd.py @@ -94,7 +94,7 @@ def _MDC( Fop = _FFT( dims=(nt, nr, nv), - dir=0, + axis=0, real=True, ifftshift_before=twosided, dtype=rdtype, @@ -102,7 +102,7 @@ def _MDC( ) F1op = _FFT( dims=(nt, ns, nv), - dir=0, + axis=0, real=True, ifftshift_before=False, dtype=rdtype, diff --git a/pylops/waveeqprocessing/oneway.py b/pylops/waveeqprocessing/oneway.py index 4f446348..b762112f 100644 --- a/pylops/waveeqprocessing/oneway.py +++ b/pylops/waveeqprocessing/oneway.py @@ -165,14 +165,14 @@ def PhaseShift(vel, dz, nt, freq, kx, ky=None, dtype="float64"): else: dims = (nt, kx.size, ky.size) dimsfft = (freq.size, kx.size, ky.size) - Fop = FFT(dims, dir=0, nfft=nt, real=True, dtype=dtype) + Fop = FFT(dims, axis=0, nfft=nt, real=True, dtype=dtype) Kxop = FFT( - dimsfft, dir=1, nfft=kx.size, real=False, fftshift_after=True, dtype=dtypefft + dimsfft, axis=1, nfft=kx.size, real=False, fftshift_after=True, dtype=dtypefft ) if ky is not None: Kyop = FFT( dimsfft, - dir=2, + axis=2, nfft=ky.size, real=False, fftshift_after=True, diff --git a/pylops/waveeqprocessing/seismicinterpolation.py b/pylops/waveeqprocessing/seismicinterpolation.py index c53d8e62..32032a00 100644 --- a/pylops/waveeqprocessing/seismicinterpolation.py +++ b/pylops/waveeqprocessing/seismicinterpolation.py @@ -226,16 +226,16 @@ def SeismicInterpolation( # create restriction/interpolation operator if iava.dtype == float: - Rop = Interp(dims, iava, dir=0, kind="linear", dtype=dtype) + Rop = Interp(dims, iava, axis=0, kind="linear", dtype=dtype) if ndims == 3 and iava1 is not None: dims1 = (len(iava), nrec[1], dimsd[2]) - Rop1 = Interp(dims1, iava1, dir=1, kind="linear", dtype=dtype) + Rop1 = Interp(dims1, iava1, axis=1, kind="linear", dtype=dtype) Rop = Rop1 * Rop else: - Rop = Restriction(dims, iava, dir=0, dtype=dtype) + Rop = Restriction(dims, iava, axis=0, dtype=dtype) if ndims == 3 and iava1 is not None: dims1 = (len(iava), nrec[1], dimsd[2]) - Rop1 = Restriction(dims1, iava1, dir=1, dtype=dtype) + Rop1 = Restriction(dims1, iava1, axis=1, dtype=dtype) Rop = Rop1 * Rop # create other operators for inversion @@ -243,9 +243,9 @@ def SeismicInterpolation( prec = False dotcflag = 0 if ndims == 3 and iava1 is not None: - Regop = Laplacian(dims=dims, dirs=(0, 1), dtype=dtype) + Regop = Laplacian(dims=dims, axes=(0, 1), dtype=dtype) else: - Regop = SecondDerivative(dims, dir=0, dtype=dtype) + Regop = SecondDerivative(dims, axis=0, dtype=dtype) SIop = Rop elif kind == "fk": prec = True diff --git a/pytests/test_basicoperators.py b/pytests/test_basicoperators.py index 4b0a6fcf..8f9ca54b 100755 --- a/pytests/test_basicoperators.py +++ b/pytests/test_basicoperators.py @@ -220,18 +220,18 @@ def test_Flip2D(par): "imag" ] * np.outer(np.ones(par["ny"]), np.arange(par["nx"])) - for dir in [0, 1]: + for axis in [0, 1]: Fop = Flip( (par["ny"], par["nx"]), - dir=dir, + axis=axis, dtype=par["dtype"], ) assert dottest(Fop, par["ny"] * par["nx"], par["ny"] * par["nx"]) - y = Fop * x[str(dir)].ravel() + y = Fop * x[str(axis)].ravel() xadj = Fop.H * y.ravel() xadj = xadj.reshape(par["ny"], par["nx"]) - assert_array_equal(x[str(dir)], xadj) + assert_array_equal(x[str(axis)], xadj) @pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j)]) @@ -268,20 +268,20 @@ def test_Flip3D(par): par["nx"] ) - for dir in [0, 1, 2]: + for axis in [0, 1, 2]: Fop = Flip( (par["ny"], par["nx"], par["nx"]), - dir=dir, + axis=axis, dtype=par["dtype"], ) assert dottest( Fop, par["ny"] * par["nx"] * par["nx"], par["ny"] * par["nx"] * par["nx"] ) - y = Fop * x[str(dir)].ravel() + y = Fop * x[str(axis)].ravel() xadj = Fop.H * y.ravel() xadj = xadj.reshape(par["ny"], par["nx"], par["nx"]) - assert_array_equal(x[str(dir)], xadj) + assert_array_equal(x[str(axis)], xadj) @pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j), (par3)]) @@ -310,17 +310,17 @@ def test_Symmetrize2D(par): "imag" ] * np.outer(np.ones(par["ny"]), np.arange(par["nx"])) - for dir in [0, 1]: + for axis in [0, 1]: Sop = Symmetrize( (par["ny"], par["nx"]), - dir=dir, + axis=axis, dtype=par["dtype"], ) - y = Sop * x[str(dir)].ravel() + y = Sop * x[str(axis)].ravel() assert dottest(Sop, y.size, par["ny"] * par["nx"]) xinv = Sop / y - assert_array_almost_equal(x[str(dir)].ravel(), xinv, decimal=3) + assert_array_almost_equal(x[str(axis)].ravel(), xinv, decimal=3) @pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j), (par3)]) @@ -357,17 +357,17 @@ def test_Symmetrize3D(par): par["nx"] ) - for dir in [0, 1, 2]: + for axis in [0, 1, 2]: Sop = Symmetrize( (par["ny"], par["nx"], par["nx"]), - dir=dir, + axis=axis, dtype=par["dtype"], ) - y = Sop * x[str(dir)].ravel() + y = Sop * x[str(axis)].ravel() assert dottest(Sop, y.size, par["ny"] * par["nx"] * par["nx"]) xinv = Sop / y - assert_array_almost_equal(x[str(dir)].ravel(), xinv, decimal=3) + assert_array_almost_equal(x[str(axis)].ravel(), xinv, decimal=3) @pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j), (par3)]) @@ -396,18 +396,18 @@ def test_Roll2D(par): "imag" ] * np.outer(np.ones(par["ny"]), np.arange(par["nx"])) - for dir in [0, 1]: + for axis in [0, 1]: Rop = Roll( (par["ny"], par["nx"]), - dir=dir, + axis=axis, shift=-2, dtype=par["dtype"], ) - y = Rop * x[str(dir)].ravel() + y = Rop * x[str(axis)].ravel() assert dottest(Rop, par["ny"] * par["nx"], par["ny"] * par["nx"]) xadj = Rop.H * y - assert_array_almost_equal(x[str(dir)].ravel(), xadj, decimal=3) + assert_array_almost_equal(x[str(axis)].ravel(), xadj, decimal=3) @pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j), (par3)]) @@ -444,39 +444,39 @@ def test_Roll3D(par): par["nx"] ) - for dir in [0, 1, 2]: + for axis in [0, 1, 2]: Rop = Roll( (par["ny"], par["nx"], par["nx"]), - dir=dir, + axis=axis, shift=3, dtype=par["dtype"], ) - y = Rop * x[str(dir)].ravel() + y = Rop * x[str(axis)].ravel() assert dottest( Rop, par["ny"] * par["nx"] * par["nx"], par["ny"] * par["nx"] * par["nx"] ) xinv = Rop.H * y - assert_array_almost_equal(x[str(dir)].ravel(), xinv, decimal=3) + assert_array_almost_equal(x[str(axis)].ravel(), xinv, decimal=3) @pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j), (par3)]) def test_Sum2D(par): """Dot-test for Sum operator on 2d signal""" - for dir in [0, 1]: + for axis in [0, 1]: dim_d = [par["ny"], par["nx"]] - dim_d.pop(dir) - Sop = Sum(dims=(par["ny"], par["nx"]), dir=dir, dtype=par["dtype"]) + dim_d.pop(axis) + Sop = Sum(dims=(par["ny"], par["nx"]), axis=axis, dtype=par["dtype"]) assert dottest(Sop, np.prod(dim_d), par["ny"] * par["nx"]) @pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j), (par3)]) def test_Sum3D(par): """Dot-test, forward and adjoint for Sum operator on 3d signal""" - for dir in [0, 1, 2]: + for axis in [0, 1, 2]: dim_d = [par["ny"], par["nx"], par["nx"]] - dim_d.pop(dir) - Sop = Sum(dims=(par["ny"], par["nx"], par["nx"]), dir=dir, dtype=par["dtype"]) + dim_d.pop(axis) + Sop = Sum(dims=(par["ny"], par["nx"], par["nx"]), axis=axis, dtype=par["dtype"]) assert dottest(Sop, np.prod(dim_d), par["ny"] * par["nx"] * par["nx"]) diff --git a/pytests/test_causalintegration.py b/pytests/test_causalintegration.py index 4698ea57..b01aebdb 100755 --- a/pytests/test_causalintegration.py +++ b/pytests/test_causalintegration.py @@ -131,7 +131,7 @@ def test_CausalIntegration2d(par): Cop = CausalIntegration( (par["nt"], par["nx"]), sampling=dt, - dir=0, + axis=0, kind=kind, halfcurrent=False, removefirst=rf, @@ -166,7 +166,7 @@ def test_CausalIntegration2d(par): # numerical derivative Dop = FirstDerivative( - (par["nt"], par["nx"]), dir=0, sampling=dt, dtype=par["dtype"] + (par["nt"], par["nx"]), axis=0, sampling=dt, dtype=par["dtype"] ) xder = Dop * y.ravel() xder = xder.reshape(par["nt"], par["nx"]) diff --git a/pytests/test_convolve.py b/pytests/test_convolve.py index 5d0556e6..d55dd2e8 100755 --- a/pytests/test_convolve.py +++ b/pytests/test_convolve.py @@ -21,42 +21,42 @@ "ny": 51, "nx": 31, "offset": nfilt[0] // 2, - "dir": 0, + "axis": 0, } # zero phase, first direction par2_1d = { "nz": 21, "ny": 61, "nx": 31, "offset": 0, - "dir": 0, + "axis": 0, } # non-zero phase, first direction par3_1d = { "nz": 21, "ny": 51, "nx": 31, "offset": nfilt[0] // 2, - "dir": 1, + "axis": 1, } # zero phase, second direction par4_1d = { "nz": 21, "ny": 61, "nx": 31, "offset": nfilt[0] // 2 - 1, - "dir": 1, + "axis": 1, } # non-zero phase, second direction par5_1d = { "nz": 21, "ny": 51, "nx": 31, "offset": nfilt[0] // 2, - "dir": 2, + "axis": 2, } # zero phase, third direction par6_1d = { "nz": 21, "ny": 61, "nx": 31, "offset": nfilt[0] // 2 - 1, - "dir": 2, + "axis": 2, } # non-zero phase, third direction par1_2d = { @@ -64,42 +64,42 @@ "ny": 51, "nx": 31, "offset": (nfilt[0] // 2, nfilt[1] // 2), - "dir": 0, + "axis": 0, } # zero phase, first direction par2_2d = { "nz": 21, "ny": 61, "nx": 31, "offset": (nfilt[0] // 2 - 1, nfilt[1] // 2 + 1), - "dir": 0, + "axis": 0, } # non-zero phase, first direction par3_2d = { "nz": 21, "ny": 51, "nx": 31, "offset": (nfilt[0] // 2, nfilt[1] // 2), - "dir": 1, + "axis": 1, } # zero phase, second direction par4_2d = { "nz": 21, "ny": 61, "nx": 31, "offset": (nfilt[0] // 2 - 1, nfilt[1] // 2 + 1), - "dir": 1, + "axis": 1, } # non-zero phase, second direction par5_2d = { "nz": 21, "ny": 51, "nx": 31, "offset": (nfilt[0] // 2, nfilt[1] // 2), - "dir": 2, + "axis": 2, } # zero phase, third direction par6_2d = { "nz": 21, "ny": 61, "nx": 31, "offset": (nfilt[0] // 2 - 1, nfilt[1] // 2 + 1), - "dir": 2, + "axis": 2, } # non-zero phase, third direction par1_3d = { @@ -108,7 +108,7 @@ "nx": 31, "nt": 5, "offset": (nfilt[0] // 2, nfilt[1] // 2, nfilt[2] // 2), - "dir": 0, + "axis": 0, } # zero phase, all directions par2_3d = { "nz": 21, @@ -116,7 +116,7 @@ "nx": 31, "nt": 5, "offset": (nfilt[0] // 2 - 1, nfilt[1] // 2 + 1, nfilt[2] // 2 + 1), - "dir": 0, + "axis": 0, } # non-zero phase, first direction @@ -127,7 +127,7 @@ def test_Convolve1D(par): """Dot-test and inversion for Convolve1D operator""" np.random.seed(10) # 1D - if par["dir"] == 0: + if par["axis"] == 0: Cop = Convolve1D(par["nx"], h=h1, offset=par["offset"], dtype="float64") assert dottest(Cop, par["nx"], par["nx"]) @@ -137,12 +137,12 @@ def test_Convolve1D(par): assert_array_almost_equal(x, xlsqr, decimal=1) # 1D on 2D - if par["dir"] < 2: + if par["axis"] < 2: Cop = Convolve1D( (par["ny"], par["nx"]), h=h1, offset=par["offset"], - dir=par["dir"], + axis=par["axis"], dtype="float64", ) assert dottest(Cop, par["ny"] * par["nx"], par["ny"] * par["nx"]) @@ -161,7 +161,7 @@ def test_Convolve1D(par): (par["nz"], par["ny"], par["nx"]), h=h1, offset=par["offset"], - dir=par["dir"], + axis=par["axis"], dtype="float64", ) assert dottest( @@ -185,7 +185,7 @@ def test_Convolve1D(par): def test_Convolve2D(par): """Dot-test and inversion for Convolve2D operator""" # 2D on 2D - if par["dir"] == 2: + if par["axis"] == 2: Cop = Convolve2D( (par["ny"], par["nx"]), h=h2, @@ -204,11 +204,13 @@ def test_Convolve2D(par): assert_array_almost_equal(x, xlsqr, decimal=1) # 2D on 3D + axes = list(range(3)) + axes.remove(par["axis"]) Cop = Convolve2D( (par["nz"], par["ny"], par["nx"]), h=h2, offset=par["offset"], - nodir=par["dir"], + axes=axes, dtype="float64", ) assert dottest( @@ -258,7 +260,7 @@ def test_Convolve3D(par): (par["nz"], par["ny"], par["nx"], par["nt"]), h=h3, offset=par["offset"], - dirs=[0, 1, 2], + axes=[0, 1, 2], dtype="float64", ) assert dottest( diff --git a/pytests/test_derivative.py b/pytests/test_derivative.py index 2270977b..63643c09 100755 --- a/pytests/test_derivative.py +++ b/pytests/test_derivative.py @@ -107,7 +107,7 @@ def test_FirstDerivative_centered(par): # 2d - derivative on 1st direction D1op = FirstDerivative( (par["ny"], par["nx"]), - dir=0, + axis=0, sampling=par["dy"], edge=par["edge"], dtype="float32", @@ -123,7 +123,7 @@ def test_FirstDerivative_centered(par): # 2d - derivative on 2nd direction D1op = FirstDerivative( (par["ny"], par["nx"]), - dir=1, + axis=1, sampling=par["dx"], edge=par["edge"], dtype="float32", @@ -139,7 +139,7 @@ def test_FirstDerivative_centered(par): # 3d - derivative on 1st direction D1op = FirstDerivative( (par["nz"], par["ny"], par["nx"]), - dir=0, + axis=0, sampling=par["dz"], edge=par["edge"], dtype="float32", @@ -164,7 +164,7 @@ def test_FirstDerivative_centered(par): # 3d - derivative on 2nd direction D1op = FirstDerivative( (par["nz"], par["ny"], par["nx"]), - dir=1, + axis=1, sampling=par["dy"], edge=par["edge"], dtype="float32", @@ -187,7 +187,7 @@ def test_FirstDerivative_centered(par): # 3d - derivative on 3rd direction D1op = FirstDerivative( (par["nz"], par["ny"], par["nx"]), - dir=2, + axis=2, sampling=par["dx"], edge=par["edge"], dtype="float32", @@ -223,7 +223,7 @@ def test_FirstDerivative_forwaback(par): # 2d - derivative on 1st direction D1op = FirstDerivative( (par["ny"], par["nx"]), - dir=0, + axis=0, sampling=par["dy"], edge=par["edge"], kind=kind, @@ -234,7 +234,7 @@ def test_FirstDerivative_forwaback(par): # 2d - derivative on 2nd direction D1op = FirstDerivative( (par["ny"], par["nx"]), - dir=1, + axis=1, sampling=par["dx"], edge=par["edge"], kind=kind, @@ -245,7 +245,7 @@ def test_FirstDerivative_forwaback(par): # 3d - derivative on 1st direction D1op = FirstDerivative( (par["nz"], par["ny"], par["nx"]), - dir=0, + axis=0, sampling=par["dz"], edge=par["edge"], kind=kind, @@ -261,7 +261,7 @@ def test_FirstDerivative_forwaback(par): # 3d - derivative on 2nd direction D1op = FirstDerivative( (par["nz"], par["ny"], par["nx"]), - dir=1, + axis=1, sampling=par["dy"], edge=par["edge"], kind=kind, @@ -277,7 +277,7 @@ def test_FirstDerivative_forwaback(par): # 3d - derivative on 3rd direction D1op = FirstDerivative( (par["nz"], par["ny"], par["nx"]), - dir=2, + axis=2, sampling=par["dx"], edge=par["edge"], kind=kind, @@ -322,7 +322,7 @@ def test_SecondDerivative_centered(par): # 2d - derivative on 1st direction D2op = SecondDerivative( (par["ny"], par["nx"]), - dir=0, + axis=0, sampling=par["dy"], edge=par["edge"], dtype="float32", @@ -340,7 +340,7 @@ def test_SecondDerivative_centered(par): # 2d - derivative on 2nd direction D2op = SecondDerivative( (par["ny"], par["nx"]), - dir=1, + axis=1, sampling=par["dx"], edge=par["edge"], dtype="float32", @@ -358,7 +358,7 @@ def test_SecondDerivative_centered(par): # 3d - derivative on 1st direction D2op = SecondDerivative( (par["ny"], par["nx"], par["nz"]), - dir=0, + axis=0, sampling=par["dy"], edge=par["edge"], dtype="float32", @@ -382,7 +382,7 @@ def test_SecondDerivative_centered(par): # 3d - derivative on 2nd direction D2op = SecondDerivative( (par["ny"], par["nx"], par["nz"]), - dir=1, + axis=1, sampling=par["dx"], edge=par["edge"], dtype="float32", @@ -406,7 +406,7 @@ def test_SecondDerivative_centered(par): # 3d - derivative on 3rd direction D2op = SecondDerivative( (par["ny"], par["nx"], par["nz"]), - dir=2, + axis=2, sampling=par["dz"], edge=par["edge"], dtype="float32", @@ -451,7 +451,7 @@ def test_SecondDerivative_forward(par): # 2d - derivative on 1st direction D2op = SecondDerivative( (par["ny"], par["nx"]), - dir=0, + axis=0, sampling=par["dy"], edge=par["edge"], kind="forward", @@ -463,7 +463,7 @@ def test_SecondDerivative_forward(par): # 2d - derivative on 2nd direction D2op = SecondDerivative( (par["ny"], par["nx"]), - dir=1, + axis=1, sampling=par["dx"], edge=par["edge"], kind="forward", @@ -475,7 +475,7 @@ def test_SecondDerivative_forward(par): # 3d - derivative on 1st direction D2op = SecondDerivative( (par["ny"], par["nx"], par["nz"]), - dir=0, + axis=0, sampling=par["dy"], edge=par["edge"], kind="forward", @@ -492,7 +492,7 @@ def test_SecondDerivative_forward(par): # 3d - derivative on 2nd direction D2op = SecondDerivative( (par["ny"], par["nx"], par["nz"]), - dir=1, + axis=1, sampling=par["dx"], edge=par["edge"], kind="forward", @@ -509,7 +509,7 @@ def test_SecondDerivative_forward(par): # 3d - derivative on 3rd direction D2op = SecondDerivative( (par["ny"], par["nx"], par["nz"]), - dir=2, + axis=2, sampling=par["dz"], edge=par["edge"], kind="forward", @@ -532,7 +532,7 @@ def test_Laplacian(par): # 2d - symmetrical Dlapop = Laplacian( (par["ny"], par["nx"]), - dirs=(0, 1), + axes=(0, 1), weights=(1, 1), sampling=(par["dy"], par["dx"]), edge=par["edge"], @@ -543,7 +543,7 @@ def test_Laplacian(par): # 2d - asymmetrical Dlapop = Laplacian( (par["ny"], par["nx"]), - dirs=(0, 1), + axes=(0, 1), weights=(1, 2), sampling=(par["dy"], par["dx"]), edge=par["edge"], @@ -554,7 +554,7 @@ def test_Laplacian(par): # 3d - symmetrical on 1st and 2nd direction Dlapop = Laplacian( (par["nz"], par["ny"], par["nx"]), - dirs=(0, 1), + axes=(0, 1), weights=(1, 1), sampling=(par["dy"], par["dx"]), edge=par["edge"], @@ -570,7 +570,7 @@ def test_Laplacian(par): # 3d - symmetrical on 1st and 2nd direction Dlapop = Laplacian( (par["nz"], par["ny"], par["nx"]), - dirs=(0, 1), + axes=(0, 1), weights=(1, 1), sampling=(par["dy"], par["dx"]), edge=par["edge"], @@ -586,7 +586,7 @@ def test_Laplacian(par): # 3d - symmetrical on all directions Dlapop = Laplacian( (par["nz"], par["ny"], par["nx"]), - dirs=(0, 1, 2), + axes=(0, 1, 2), weights=(1, 1, 1), sampling=(par["dz"], par["dx"], par["dx"]), edge=par["edge"], @@ -705,7 +705,7 @@ def test_SecondDirectionalDerivative_verticalderivative(par): and SecondDerivative """ Fop = FirstDerivative( - (par["ny"], par["nx"]), dir=0, edge=par["edge"], dtype="float32" + (par["ny"], par["nx"]), axis=0, edge=par["edge"], dtype="float32" ) F2op = Fop.H * Fop diff --git a/pytests/test_diagonal.py b/pytests/test_diagonal.py index 65625dee..6bdf93f5 100755 --- a/pytests/test_diagonal.py +++ b/pytests/test_diagonal.py @@ -33,7 +33,7 @@ def test_Diagonal_2dsignal(par): for idim, ddim in enumerate((par["nx"], par["nt"])): d = np.arange(ddim) + 1.0 + par["imag"] * (np.arange(ddim) + 1.0) - Dop = Diagonal(d, dims=(par["nx"], par["nt"]), dir=idim, dtype=par["dtype"]) + Dop = Diagonal(d, dims=(par["nx"], par["nt"]), axis=idim, dtype=par["dtype"]) assert dottest( Dop, par["nx"] * par["nt"], @@ -56,7 +56,7 @@ def test_Diagonal_3dsignal(par): d = np.arange(ddim) + 1.0 + par["imag"] * (np.arange(ddim) + 1.0) Dop = Diagonal( - d, dims=(par["ny"], par["nx"], par["nt"]), dir=idim, dtype=par["dtype"] + d, dims=(par["ny"], par["nx"], par["nt"]), axis=idim, dtype=par["dtype"] ) assert dottest( Dop, diff --git a/pytests/test_dwts.py b/pytests/test_dwts.py index a4881430..4f0fd407 100755 --- a/pytests/test_dwts.py +++ b/pytests/test_dwts.py @@ -22,7 +22,7 @@ def test_unknown_wavelet(par): @pytest.mark.parametrize("par", [(par1), (par2)]) def test_DWT_1dsignal(par): """Dot-test and inversion for DWT operator for 1d signal""" - DWTop = DWT(dims=[par["nt"]], dir=0, wavelet="haar", level=3) + DWTop = DWT(dims=[par["nt"]], axis=0, wavelet="haar", level=3) x = np.random.normal(0.0, 1.0, par["nt"]) + par["imag"] * np.random.normal( 0.0, 1.0, par["nt"] ) @@ -42,8 +42,8 @@ def test_DWT_1dsignal(par): @pytest.mark.parametrize("par", [(par1), (par2)]) def test_DWT_2dsignal(par): """Dot-test and inversion for DWT operator for 2d signal""" - for dir in [0, 1]: - DWTop = DWT(dims=(par["nt"], par["nx"]), dir=dir, wavelet="haar", level=3) + for axis in [0, 1]: + DWTop = DWT(dims=(par["nt"], par["nx"]), axis=axis, wavelet="haar", level=3) x = np.random.normal(0.0, 1.0, (par["nt"], par["nx"])) + par[ "imag" ] * np.random.normal(0.0, 1.0, (par["nt"], par["nx"])) @@ -66,9 +66,9 @@ def test_DWT_2dsignal(par): @pytest.mark.parametrize("par", [(par1), (par2)]) def test_DWT_3dsignal(par): """Dot-test and inversion for DWT operator for 3d signal""" - for dir in [0, 1, 2]: + for axis in [0, 1, 2]: DWTop = DWT( - dims=(par["nt"], par["nx"], par["ny"]), dir=dir, wavelet="haar", level=3 + dims=(par["nt"], par["nx"], par["ny"]), axis=axis, wavelet="haar", level=3 ) x = np.random.normal(0.0, 1.0, (par["nt"], par["nx"], par["ny"])) + par[ "imag" @@ -92,7 +92,7 @@ def test_DWT_3dsignal(par): @pytest.mark.parametrize("par", [(par1), (par2)]) def test_DWT2D_2dsignal(par): """Dot-test and inversion for DWT2D operator for 2d signal""" - DWTop = DWT2D(dims=(par["nt"], par["nx"]), dirs=(0, 1), wavelet="haar", level=3) + DWTop = DWT2D(dims=(par["nt"], par["nx"]), axes=(0, 1), wavelet="haar", level=3) x = np.random.normal(0.0, 1.0, (par["nt"], par["nx"])) + par[ "imag" ] * np.random.normal(0.0, 1.0, (par["nt"], par["nx"])) @@ -112,9 +112,9 @@ def test_DWT2D_2dsignal(par): @pytest.mark.parametrize("par", [(par1), (par2)]) def test_DWT2D_3dsignal(par): """Dot-test and inversion for DWT operator for 3d signal""" - for dirs in [(0, 1), (0, 2), (1, 2)]: + for axes in [(0, 1), (0, 2), (1, 2)]: DWTop = DWT2D( - dims=(par["nt"], par["nx"], par["ny"]), dirs=dirs, wavelet="haar", level=3 + dims=(par["nt"], par["nx"], par["ny"]), axes=axes, wavelet="haar", level=3 ) x = np.random.normal(0.0, 1.0, (par["nt"], par["nx"], par["ny"])) + par[ "imag" diff --git a/pytests/test_ffts.py b/pytests/test_ffts.py index 214b5ba7..fea72128 100755 --- a/pytests/test_ffts.py +++ b/pytests/test_ffts.py @@ -186,7 +186,7 @@ def test_FFT_small_real(par): FFTop = FFT( dims=x.shape, - dir=0, + axis=0, norm=norm, real=True, ifftshift_before=ifftshift_before, @@ -259,7 +259,7 @@ def test_FFT_random_real(par): FFTop = FFT( dims=x.shape, - dir=axis, + axis=axis, ifftshift_before=ifftshift_before, real=True, dtype=dtype, @@ -303,7 +303,7 @@ def test_FFT_small_complex(par): FFTop = FFT( dims=x.shape, - dir=0, + axis=0, norm=norm, ifftshift_before=ifftshift_before, fftshift_after=fftshift_after, @@ -377,7 +377,7 @@ def test_FFT_random_complex(par): FFTop = FFT( dims=x.shape, - dir=axis, + axis=axis, ifftshift_before=ifftshift_before, fftshift_after=fftshift_after, dtype=dtype, @@ -453,7 +453,7 @@ def test_FFT2D_random_real(par): FFTop = FFT2D( dims=x.shape, - dirs=axes, + axes=axes, ifftshift_before=ifftshift_before, real=True, dtype=dtype, @@ -519,7 +519,7 @@ def test_FFT2D_random_complex(par): FFTop = FFT2D( dims=x.shape, - dirs=axes, + axes=axes, ifftshift_before=ifftshift_before, fftshift_after=fftshift_after, dtype=dtype, @@ -595,7 +595,7 @@ def test_FFTND_random_real(par): FFTop = FFTND( dims=x.shape, - dirs=axes, + axes=axes, ifftshift_before=ifftshift_before, real=True, dtype=dtype, @@ -661,7 +661,7 @@ def test_FFTND_random_complex(par): FFTop = FFTND( dims=x.shape, - dirs=axes, + axes=axes, ifftshift_before=ifftshift_before, fftshift_after=fftshift_after, dtype=dtype, @@ -726,7 +726,7 @@ def test_FFT2D_small_complex(par): FFTop = FFT2D( dims=x.shape, - dirs=(0, 1), + axes=(0, 1), norm=norm, dtype=dtype, ) @@ -773,7 +773,7 @@ def test_FFTND_small_complex(par): FFTop = FFTND( dims=x.shape, - dirs=(0, 1), + axes=(0, 1), norm=norm, dtype=dtype, ) @@ -912,7 +912,7 @@ def test_FFT_2dsignal(par): nfft = par["nt"] if par["nfft"] is None else par["nfft"] FFTop = FFT( dims=(nt, nx), - dir=0, + axis=0, nfft=nfft, sampling=dt, real=par["real"], @@ -943,7 +943,7 @@ def test_FFT_2dsignal(par): if not par["real"]: FFTop_fftshift = FFT( dims=(nt, nx), - dir=0, + axis=0, nfft=nfft, sampling=dt, real=par["real"], @@ -970,7 +970,7 @@ def test_FFT_2dsignal(par): nfft = par["nx"] if par["nfft"] is None else par["nfft"] FFTop = FFT( dims=(nt, nx), - dir=1, + axis=1, nfft=nfft, sampling=dt, real=par["real"], @@ -1001,7 +1001,7 @@ def test_FFT_2dsignal(par): if not par["real"]: FFTop_fftshift = FFT( dims=(nt, nx), - dir=1, + axis=1, nfft=nfft, sampling=dt, real=par["real"], @@ -1059,7 +1059,7 @@ def test_FFT_3dsignal(par): nfft = par["nt"] if par["nfft"] is None else par["nfft"] FFTop = FFT( dims=(nt, nx, ny), - dir=0, + axis=0, nfft=nfft, sampling=dt, real=par["real"], @@ -1099,7 +1099,7 @@ def test_FFT_3dsignal(par): nfft = par["nx"] if par["nfft"] is None else par["nfft"] FFTop = FFT( dims=(nt, nx, ny), - dir=1, + axis=1, nfft=nfft, sampling=dt, real=par["real"], @@ -1140,7 +1140,7 @@ def test_FFT_3dsignal(par): nfft = par["ny"] if par["nfft"] is None else par["nfft"] FFTop = FFT( dims=(nt, nx, ny), - dir=2, + axis=2, nfft=nfft, sampling=dt, real=par["real"], @@ -1179,7 +1179,7 @@ def test_FFT_3dsignal(par): if not par["real"]: FFTop_fftshift = FFT( dims=(nt, nx, ny), - dir=2, + axis=2, nfft=nfft, sampling=dt, real=par["real"], @@ -1216,13 +1216,13 @@ def test_FFT2D(par): d = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(par["nx"]) + 1) d = d.astype(par["dtype"]) - # first fft on dir 1 + # first fft on axis 1 FFTop = FFT2D( dims=(par["nt"], par["nx"]), nffts=(nfft1, nfft2), sampling=(dt, dx), real=par["real"], - dirs=(0, 1), + axes=(0, 1), ) if par["real"]: @@ -1262,13 +1262,13 @@ def test_FFT2D(par): assert_array_almost_equal(d[:imax1, :imax2], dadj[:imax1, :imax2], decimal=decimal) assert_array_almost_equal(d[:imax1, :imax2], dinv[:imax1, :imax2], decimal=decimal) - # first fft on dir 0 + # first fft on axis 0 FFTop = FFT2D( dims=(par["nt"], par["nx"]), nffts=(nfft2, nfft1), sampling=(dx, dt), real=par["real"], - dirs=(1, 0), + axes=(1, 0), ) if par["real"]: @@ -1322,11 +1322,11 @@ def test_FFT3D(par): d = np.tile(d[:, :, np.newaxis], [1, 1, par["ny"]]) d = d.astype(par["dtype"]) - # first fft on dir 2 + # first fft on axis 2 FFTop = FFTND( dims=(par["nt"], par["nx"], par["ny"]), nffts=(nfft1, nfft2, nfft3), - dirs=(0, 1, 2), + axes=(0, 1, 2), sampling=(dt, dx, dy), real=par["real"], ) @@ -1373,11 +1373,11 @@ def test_FFT3D(par): d[:imax1, :imax2, :imax3], dinv[:imax1, :imax2, :imax3], decimal=decimal ) - # first fft on dir 1 + # first fft on axis 1 FFTop = FFTND( dims=(par["nt"], par["nx"], par["ny"]), nffts=(nfft1, nfft3, nfft2), - dirs=(0, 2, 1), + axes=(0, 2, 1), sampling=(dt, dy, dx), real=par["real"], ) @@ -1420,11 +1420,11 @@ def test_FFT3D(par): d[:imax1, :imax2, :imax3], dinv[:imax1, :imax2, :imax3], decimal=decimal ) - # first fft on dir 0 + # first fft on axis 0 FFTop = FFTND( dims=(par["nt"], par["nx"], par["ny"]), nffts=(nfft2, nfft3, nfft1), - dirs=(1, 2, 0), + axes=(1, 2, 0), sampling=(dx, dy, dt), real=par["real"], ) diff --git a/pytests/test_interpolation.py b/pytests/test_interpolation.py index b8054bc9..593521d1 100755 --- a/pytests/test_interpolation.py +++ b/pytests/test_interpolation.py @@ -140,7 +140,7 @@ def test_Interp_2dsignal(par): Iop, _ = Interp( (par["nx"], par["nt"]), iava, - dir=0, + axis=0, kind=par["kind"], dtype=par["dtype"], ) @@ -155,7 +155,7 @@ def test_Interp_2dsignal(par): Idecop, _ = Interp( (par["nx"], par["nt"]), iava + 0.3, - dir=0, + axis=0, kind=par["kind"], dtype=par["dtype"], ) @@ -168,7 +168,7 @@ def test_Interp_2dsignal(par): _, _ = Interp( (par["nx"], par["nt"]), iava_rep + 0.3, - dir=0, + axis=0, kind=par["kind"], dtype=par["dtype"], ) @@ -188,7 +188,7 @@ def test_Interp_2dsignal(par): Iop, _ = Interp( (par["nx"], par["nt"]), iava, - dir=1, + axis=1, kind=par["kind"], dtype=par["dtype"], ) @@ -203,7 +203,7 @@ def test_Interp_2dsignal(par): Idecop, _ = Interp( (par["nx"], par["nt"]), iava + 0.3, - dir=1, + axis=1, kind=par["kind"], dtype=par["dtype"], ) @@ -238,7 +238,7 @@ def test_Interp_3dsignal(par): Iop, _ = Interp( (par["ny"], par["nx"], par["nt"]), iava, - dir=0, + axis=0, kind=par["kind"], dtype=par["dtype"], ) @@ -253,7 +253,7 @@ def test_Interp_3dsignal(par): Idecop, _ = Interp( (par["ny"], par["nx"], par["nt"]), iava + 0.3, - dir=0, + axis=0, kind=par["kind"], dtype=par["dtype"], ) @@ -272,7 +272,7 @@ def test_Interp_3dsignal(par): _, _ = Interp( (par["ny"], par["nx"], par["nt"]), iava_rep + 0.3, - dir=0, + axis=0, kind=par["kind"], dtype=par["dtype"], ) @@ -292,7 +292,7 @@ def test_Interp_3dsignal(par): Iop, _ = Interp( (par["ny"], par["nx"], par["nt"]), iava, - dir=1, + axis=1, kind=par["kind"], dtype=par["dtype"], ) @@ -307,7 +307,7 @@ def test_Interp_3dsignal(par): Idecop, _ = Interp( (par["ny"], par["nx"], par["nt"]), iava + 0.3, - dir=1, + axis=1, kind=par["kind"], dtype=par["dtype"], ) @@ -333,7 +333,7 @@ def test_Interp_3dsignal(par): Iop, _ = Interp( (par["ny"], par["nx"], par["nt"]), iava, - dir=2, + axis=2, kind=par["kind"], dtype=par["dtype"], ) @@ -348,7 +348,7 @@ def test_Interp_3dsignal(par): Idecop, _ = Interp( (par["ny"], par["nx"], par["nt"]), iava + 0.3, - dir=2, + axis=2, kind=par["kind"], dtype=par["dtype"], ) diff --git a/pytests/test_kronecker.py b/pytests/test_kronecker.py index beb8b973..08a342b9 100755 --- a/pytests/test_kronecker.py +++ b/pytests/test_kronecker.py @@ -39,11 +39,11 @@ def test_Kroneker(par): @pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j)]) def test_Kroneker_Derivative(par): """Use Kronecker operator to apply the Derivative operator over one axis - and compare with FirstDerivative(... dir=axis) + and compare with FirstDerivative(... axis=axis) """ Dop = FirstDerivative(par["ny"], sampling=1, edge=True, dtype="float32") D2op = FirstDerivative( - (par["ny"], par["nx"]), dir=0, sampling=1, edge=True, dtype="float32" + (par["ny"], par["nx"]), axis=0, sampling=1, edge=True, dtype="float32" ) Kop = Kronecker(Dop, Identity(par["nx"], dtype=par["dtype"]), dtype=par["dtype"]) diff --git a/pytests/test_restriction.py b/pytests/test_restriction.py index 21b820c3..0b132ff6 100755 --- a/pytests/test_restriction.py +++ b/pytests/test_restriction.py @@ -74,7 +74,7 @@ def test_Restriction_2dsignal(par): iava = np.sort(np.random.permutation(np.arange(par["nx"]))[:Nsub]) Rop = Restriction( - (par["nx"], par["nt"]), iava, dir=0, inplace=par["dtype"], dtype=par["dtype"] + (par["nx"], par["nt"]), iava, axis=0, inplace=par["dtype"], dtype=par["dtype"] ) assert dottest( Rop, @@ -97,7 +97,7 @@ def test_Restriction_2dsignal(par): iava = np.sort(np.random.permutation(np.arange(par["nt"]))[:Nsub]) Rop = Restriction( - (par["nx"], par["nt"]), iava, dir=1, inplace=par["dtype"], dtype=par["dtype"] + (par["nx"], par["nt"]), iava, axis=1, inplace=par["dtype"], dtype=par["dtype"] ) assert dottest( Rop, @@ -132,7 +132,7 @@ def test_Restriction_3dsignal(par): Rop = Restriction( (par["ny"], par["nx"], par["nt"]), iava, - dir=0, + axis=0, inplace=par["dtype"], dtype=par["dtype"], ) @@ -161,7 +161,7 @@ def test_Restriction_3dsignal(par): Rop = Restriction( (par["ny"], par["nx"], par["nt"]), iava, - dir=1, + axis=1, inplace=par["dtype"], dtype=par["dtype"], ) @@ -188,7 +188,7 @@ def test_Restriction_3dsignal(par): Rop = Restriction( (par["ny"], par["nx"], par["nt"]), iava, - dir=2, + axis=2, inplace=par["dtype"], dtype=par["dtype"], ) diff --git a/pytests/test_seismicinterpolation.py b/pytests/test_seismicinterpolation.py index 1442e58f..d98a8403 100755 --- a/pytests/test_seismicinterpolation.py +++ b/pytests/test_seismicinterpolation.py @@ -48,10 +48,10 @@ _, x3d = linear3d(xaxis, yaxis, taxis, v, t0, theta, phi, amp, wav) # Create restriction operator -Rop2d = Restriction((par["ny"], par["nt"]), iava, dir=0, dtype="float64") +Rop2d = Restriction((par["ny"], par["nt"]), iava, axis=0, dtype="float64") y2d = Rop2d * x2d.ravel() y2d = y2d.reshape(nysub, par["nt"]) -Rop3d = Restriction((par["ny"], par["nx"], par["nt"]), iava, dir=0, dtype="float64") +Rop3d = Restriction((par["ny"], par["nx"], par["nt"]), iava, axis=0, dtype="float64") y3d = Rop3d * x3d.ravel() y3d = y3d.reshape(nysub, par["nx"], par["nt"]) diff --git a/pytests/test_shift.py b/pytests/test_shift.py index 999b34b9..06bd79b5 100755 --- a/pytests/test_shift.py +++ b/pytests/test_shift.py @@ -56,7 +56,7 @@ def test_Shift2D(par): np.random.seed(0) shift = 5.5 - # 1st dir + # 1st axis x = ( gaussian(np.arange(par["nt"] // 2 + 1), 2.0)[0] + par["imag"] * gaussian(np.arange(par["nt"] // 2 + 1), 2.0)[0] @@ -65,7 +65,7 @@ def test_Shift2D(par): Sop = Shift( (par["nt"], par["nx"]), shift, - dir=0, + axis=0, real=True if par["imag"] == 0 else False, dtype=par["dtype"], ) @@ -78,7 +78,7 @@ def test_Shift2D(par): xlsqr = lsqr(Sop, Sop * x.ravel(), damp=1e-20, iter_lim=200, show=0)[0] assert_array_almost_equal(x.ravel(), xlsqr, decimal=1) - # 2nd dir + # 2nd axis x = ( gaussian(np.arange(par["nt"] // 2 + 1), 2.0)[0] + par["imag"] * gaussian(np.arange(par["nt"] // 2 + 1), 2.0)[0] @@ -87,7 +87,7 @@ def test_Shift2D(par): Sop = Shift( (par["nx"], par["nt"]), shift, - dir=1, + axis=1, real=True if par["imag"] == 0 else False, dtype=par["dtype"], ) diff --git a/pytests/test_smoothing.py b/pytests/test_smoothing.py index 615ddee7..bf32b0fd 100755 --- a/pytests/test_smoothing.py +++ b/pytests/test_smoothing.py @@ -6,12 +6,12 @@ from pylops.basicoperators import Smoothing1D, Smoothing2D from pylops.utils import dottest -par1 = {"nz": 10, "ny": 30, "nx": 20, "dir": 0} # even, first direction -par2 = {"nz": 11, "ny": 51, "nx": 31, "dir": 0} # odd, first direction -par3 = {"nz": 10, "ny": 30, "nx": 20, "dir": 1} # even, second direction -par4 = {"nz": 11, "ny": 51, "nx": 31, "dir": 1} # odd, second direction -par5 = {"nz": 10, "ny": 30, "nx": 20, "dir": 2} # even, third direction -par6 = {"nz": 11, "ny": 51, "nx": 31, "dir": 2} # odd, third direction +par1 = {"nz": 10, "ny": 30, "nx": 20, "axis": 0} # even, first direction +par2 = {"nz": 11, "ny": 51, "nx": 31, "axis": 0} # odd, first direction +par3 = {"nz": 10, "ny": 30, "nx": 20, "axis": 1} # even, second direction +par4 = {"nz": 11, "ny": 51, "nx": 31, "axis": 1} # odd, second direction +par5 = {"nz": 10, "ny": 30, "nx": 20, "axis": 2} # even, third direction +par6 = {"nz": 11, "ny": 51, "nx": 31, "axis": 2} # odd, third direction np.random.seed(0) @@ -30,7 +30,7 @@ def test_Smoothing1D(par): # 1d kernel on 2d signal D1op = Smoothing1D( - nsmooth=5, dims=(par["ny"], par["nx"]), dir=par["dir"], dtype="float64" + nsmooth=5, dims=(par["ny"], par["nx"]), axis=par["axis"], dtype="float64" ) assert dottest(D1op, par["ny"] * par["nx"], par["ny"] * par["nx"]) @@ -43,7 +43,7 @@ def test_Smoothing1D(par): D1op = Smoothing1D( nsmooth=5, dims=(par["nz"], par["ny"], par["nx"]), - dir=par["dir"], + axis=par["axis"], dtype="float64", ) assert dottest( @@ -63,7 +63,7 @@ def test_Smoothing1D(par): def test_Smoothing2D(par): """Dot-test for smoothing""" # 2d kernel on 2d signal - if par["dir"] < 2: + if par["axis"] < 2: D2op = Smoothing2D(nsmooth=(5, 5), dims=(par["ny"], par["nx"]), dtype="float64") assert dottest(D2op, par["ny"] * par["nx"], par["ny"] * par["nx"], rtol=1e-3) @@ -85,10 +85,12 @@ def test_Smoothing2D(par): assert_array_almost_equal(x, xlsqr, decimal=1) # 2d kernel on 3d signal + axes = list(range(3)) + axes.remove(par["axis"]) D2op = Smoothing2D( nsmooth=(5, 5), dims=(par["nz"], par["ny"], par["nx"]), - nodir=par["dir"], + axes=axes, dtype="float64", ) assert dottest( @@ -102,7 +104,7 @@ def test_Smoothing2D(par): y = D2op * x y = y.reshape(par["nz"], par["ny"], par["nx"]) - if par["dir"] == 0: + if par["axis"] == 0: assert_array_almost_equal( y[par["nz"] // 2, par["ny"] // 2 - 2 : par["ny"] // 2 + 3, par["nx"] // 2], np.ones(5) / 25, @@ -111,7 +113,7 @@ def test_Smoothing2D(par): y[par["nz"] // 2, par["ny"] // 2, par["nx"] // 2 - 2 : par["nx"] // 2 + 3], np.ones(5) / 25, ) - elif par["dir"] == 1: + elif par["axis"] == 1: assert_array_almost_equal( y[par["nz"] // 2 - 2 : par["nz"] // 2 + 3, par["ny"] // 2, par["nx"] // 2], np.ones(5) / 25, @@ -120,7 +122,7 @@ def test_Smoothing2D(par): y[par["nz"] // 2, par["ny"] // 2, par["nx"] // 2 - 2 : par["nx"] // 2 + 3], np.ones(5) / 25, ) - elif par["dir"] == 2: + elif par["axis"] == 2: assert_array_almost_equal( y[par["nz"] // 2 - 2 : par["nz"] // 2 + 3, par["ny"] // 2, par["nx"] // 2], np.ones(5) / 25, diff --git a/tutorials/ctscan.py b/tutorials/ctscan.py index adca860e..44baa2b5 100755 --- a/tutorials/ctscan.py +++ b/tutorials/ctscan.py @@ -93,10 +93,10 @@ def radoncurve(x, r, theta): # modelling operator both in a least-squares sense and using TV-reg. Dop = [ pylops.FirstDerivative( - (nx, ny), dir=0, edge=True, kind="backward", dtype=np.float64 + (nx, ny), axis=0, edge=True, kind="backward", dtype=np.float64 ), pylops.FirstDerivative( - (nx, ny), dir=1, edge=True, kind="backward", dtype=np.float64 + (nx, ny), axis=1, edge=True, kind="backward", dtype=np.float64 ), ] D2op = pylops.Laplacian(dims=(nx, ny), edge=True, dtype=np.float64) diff --git a/tutorials/deblurring.py b/tutorials/deblurring.py index 959b6444..95f8e7d8 100755 --- a/tutorials/deblurring.py +++ b/tutorials/deblurring.py @@ -57,8 +57,8 @@ Wop = pylops.signalprocessing.DWT2D((Nz, Nx), wavelet="haar", level=3) Dop = [ - pylops.FirstDerivative((Nz, Nx), dir=0, edge=False), - pylops.FirstDerivative((Nz, Nx), dir=1, edge=False), + pylops.FirstDerivative((Nz, Nx), axis=0, edge=False), + pylops.FirstDerivative((Nz, Nx), axis=1, edge=False), ] DWop = Dop + [ Wop, diff --git a/tutorials/seismicinterpolation.py b/tutorials/seismicinterpolation.py index ef52b5cf..a2fa1ec0 100755 --- a/tutorials/seismicinterpolation.py +++ b/tutorials/seismicinterpolation.py @@ -62,7 +62,7 @@ iava = np.sort(np.random.permutation(np.arange(par["nx"]))[:nxsub]) # restriction operator -Rop = pylops.Restriction((par["nx"], par["nt"]), iava, dir=0, dtype="float64") +Rop = pylops.Restriction((par["nx"], par["nt"]), iava, axis=0, dtype="float64") # data y = Rop * x.ravel() @@ -110,7 +110,7 @@ # \epsilon \|\mathbf{F}^H \mathbf{x}\|_1 # smooth inversion -D2op = pylops.SecondDerivative((par["nx"], par["nt"]), dir=0, dtype="float64") +D2op = pylops.SecondDerivative((par["nx"], par["nt"]), axis=0, dtype="float64") xsmooth, _, _ = pylops.waveeqprocessing.SeismicInterpolation( y, @@ -314,7 +314,7 @@ iava = np.sort(np.random.permutation(np.arange(par["nx"]))[:Nsub]) # restriction operator -Rop = pylops.Restriction((par["nx"], par["nt"]), iava, dir=0, dtype="float64") +Rop = pylops.Restriction((par["nx"], par["nt"]), iava, axis=0, dtype="float64") y = Rop * x.ravel() xadj = Rop.H * y.ravel()