From 493e4d485c618bea64f323bc1249b8c582292d80 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Wed, 8 May 2024 22:12:27 +0300 Subject: [PATCH 01/11] feature: improved Sliding1d Changed sliding1d from using other operators to being directly implemented (including also an option to apply Op simultaneously to all windows). --- pylops/signalprocessing/sliding1d.py | 173 ++++++++++++++++++--------- pylops/utils/backend.py | 24 ++++ 2 files changed, 138 insertions(+), 59 deletions(-) diff --git a/pylops/signalprocessing/sliding1d.py b/pylops/signalprocessing/sliding1d.py index 5cb8c213..c9db1373 100644 --- a/pylops/signalprocessing/sliding1d.py +++ b/pylops/signalprocessing/sliding1d.py @@ -6,10 +6,18 @@ import logging from typing import Tuple, Union +import numpy as np +from numpy.lib.stride_tricks import sliding_window_view + from pylops import LinearOperator -from pylops.basicoperators import BlockDiag, Diagonal, HStack, Restriction from pylops.signalprocessing.sliding2d import _slidingsteps from pylops.utils._internal import _value_or_sized_to_tuple +from pylops.utils.backend import ( + get_array_module, + get_sliding_window_view, + to_cupy_conditional, +) +from pylops.utils.decorators import reshaped from pylops.utils.tapers import taper from pylops.utils.typing import InputDimsLike, NDArray @@ -77,15 +85,7 @@ def sliding1d_design( return nwins, dim, mwins_inends, dwins_inends -def Sliding1D( - Op: LinearOperator, - dim: Union[int, InputDimsLike], - dimd: Union[int, InputDimsLike], - nwin: int, - nover: int, - tapertype: str = "hanning", - name: str = "S", -) -> LinearOperator: +class Sliding1D(LinearOperator): r"""1D Sliding transform operator. Apply a transform operator ``Op`` repeatedly to slices of the model @@ -103,6 +103,12 @@ def Sliding1D( ``nover``, it is recommended to first run ``sliding1d_design`` to obtain the corresponding ``dims`` and number of windows. + .. note:: Two kind of operators ``Op`` can be provided: the first + applies a single transformation to each window separately; the second + applies the transformation to all of the windows at the same time. This + is directly inferred during initialization when the following condition + holds ``Op.shape[1] == dim[0]``. + .. warning:: Depending on the choice of `nwin` and `nover` as well as the size of the data, sliding windows may not cover the entire data. The start and end indices of each window will be displayed and returned @@ -127,11 +133,6 @@ def Sliding1D( Name of operator (to be used by :func:`pylops.utils.describe.describe`) - Returns - ------- - Sop : :obj:`pylops.LinearOperator` - Sliding operator - Raises ------ ValueError @@ -139,50 +140,104 @@ def Sliding1D( shape (``dims``). """ - dim: Tuple[int, ...] = _value_or_sized_to_tuple(dim) - dimd: Tuple[int, ...] = _value_or_sized_to_tuple(dimd) - # data windows - dwin_ins, dwin_ends = _slidingsteps(dimd[0], nwin, nover) - nwins = len(dwin_ins) - - # check windows - if nwins * Op.shape[1] != dim[0]: - raise ValueError( - f"Model shape (dim={dim}) is not consistent with chosen " - f"number of windows. Run sliding1d_design to identify the " - f"correct number of windows for the current " - "model size..." + def __init__( + self, + Op: LinearOperator, + dim: Union[int, InputDimsLike], + dimd: Union[int, InputDimsLike], + nwin: int, + nover: int, + tapertype: str = "hanning", + name: str = "S", + ) -> None: + + dim: Tuple[int, ...] = _value_or_sized_to_tuple(dim) + dimd: Tuple[int, ...] = _value_or_sized_to_tuple(dimd) + + # data windows + dwin_ins, dwin_ends = _slidingsteps(dimd[0], nwin, nover) + self.dwin_inends = (dwin_ins, dwin_ends) + nwins = len(dwin_ins) + self.nwin = nwin + self.nover = nover + + # check windows + if nwins * Op.shape[1] != dim[0] and Op.shape[1] != dim[0]: + raise ValueError( + f"Model shape (dim={dim}) is not consistent with chosen " + f"number of windows. Run sliding1d_design to identify the " + f"correct number of windows for the current " + "model size..." + ) + + # create tapers + self.tapertype = tapertype + if self.tapertype is not None: + tap = taper(nwin, nover, tapertype=self.tapertype) + tapin = tap.copy() + tapin[:nover] = 1 + tapend = tap.copy() + tapend[-nover:] = 1 + self.taps = [ + tapin, + ] + for i in range(1, nwins - 1): + self.taps.append(tap) + self.taps.append(tapend) + self.taps = np.vstack(self.taps) + + # check if operator is applied to all windows simultaneously + self.simOp = False + if Op.shape[1] == dim[0]: + self.simOp = True + self.Op = Op + + # create temporary shape and strides for cpy + self.shape_wins = None + self.strides_wins = None + + super().__init__( + dtype=Op.dtype, + dims=(nwins, int(dim[0] // nwins)), + dimsd=dimd, + clinear=False, + name=name, ) - # create tapers - if tapertype is not None: - tap = taper(nwin, nover, tapertype=tapertype).astype(Op.dtype) - tapin = tap.copy() - tapin[:nover] = 1 - tapend = tap.copy() - tapend[-nover:] = 1 - taps = {} - taps[0] = tapin - for i in range(1, nwins - 1): - taps[i] = tap - taps[nwins - 1] = tapend - - # transform to apply - if tapertype is None: - OOp = BlockDiag([Op for _ in range(nwins)]) - else: - OOp = BlockDiag( - [Diagonal(taps[itap].ravel(), dtype=Op.dtype) * Op for itap in range(nwins)] - ) - - combining = HStack( - [ - Restriction(dimd, range(win_in, win_end), dtype=Op.dtype).H - for win_in, win_end in zip(dwin_ins, dwin_ends) - ] - ) - Sop = LinearOperator(combining * OOp) - Sop.dims, Sop.dimsd = (nwins, int(dim[0] // nwins)), dimd - Sop.name = name - return Sop + @reshaped + def _matvec(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + if self.tapertype is not None: + self.taps = to_cupy_conditional(x, self.taps) + y = ncp.zeros(self.dimsd, dtype=self.dtype) + if self.simOp: + x = self.Op @ x + for iwin0 in range(self.dims[0]): + if self.simOp: + xx = x[iwin0] + else: + xx = self.Op.matvec(x[iwin0]) + if self.tapertype is not None: + xxwin = self.taps[iwin0] * xx + else: + xxwin = xx + y[self.dwin_inends[0][iwin0] : self.dwin_inends[1][iwin0]] += xxwin + return y + + @reshaped + def _rmatvec(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + ncp_sliding_window_view = get_sliding_window_view(x) + if self.tapertype is not None: + self.taps = to_cupy_conditional(x, self.taps) + ywins = ncp_sliding_window_view(x, self.nwin)[:: self.nwin - self.nover] + if self.tapertype is not None: + ywins = ywins * self.taps + if self.simOp: + y = self.Op.H @ ywins + else: + y = ncp.zeros(self.dims, dtype=self.dtype) + for iwin0 in range(self.dims[0]): + y[iwin0] = self.Op.rmatvec(ywins[iwin0]) + return y diff --git a/pylops/utils/backend.py b/pylops/utils/backend.py index 1cef8bdb..00bdb064 100644 --- a/pylops/utils/backend.py +++ b/pylops/utils/backend.py @@ -7,6 +7,7 @@ "get_oaconvolve", "get_correlate", "get_add_at", + "get_sliding_window_view", "get_block_diag", "get_toeplitz", "get_csc_matrix", @@ -228,6 +229,29 @@ def get_add_at(x: npt.ArrayLike) -> Callable: return cupyx.scatter_add +def get_sliding_window_view(x: npt.ArrayLike) -> Callable: + """Returns correct sliding_window_view module based on input + + Parameters + ---------- + x : :obj:`numpy.ndarray` + Array + + Returns + ------- + mod : :obj:`func` + Module to be used to process array (:mod:`numpy` or :mod:`cupy`) + + """ + if not deps.cupy_enabled: + return np.lib.stride_tricks.sliding_window_view + + if cp.get_array_module(x) == np: + return np.lib.stride_tricks.sliding_window_view + else: + return cp.lib.stride_tricks.sliding_window_view + + def get_block_diag(x: npt.ArrayLike) -> Callable: """Returns correct block_diag module based on input From 2df59f08a4e7abc4e14c821b2e3a62af52d10948 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sun, 12 May 2024 22:28:27 +0300 Subject: [PATCH 02/11] feature: improved Sliding2d and Sliding3d Changed sliding2d and sliding3d from using other operators to being directly implemented (including also an option to apply Op simultaneously to all windows). --- docs/source/gpu.rst | 2 +- pylops/signalprocessing/sliding1d.py | 5 - pylops/signalprocessing/sliding2d.py | 163 +++++++++++------ pylops/signalprocessing/sliding3d.py | 251 +++++++++++++++++++-------- 4 files changed, 291 insertions(+), 130 deletions(-) diff --git a/docs/source/gpu.rst b/docs/source/gpu.rst index 864451c4..e861c567 100755 --- a/docs/source/gpu.rst +++ b/docs/source/gpu.rst @@ -5,7 +5,7 @@ GPU Support Overview -------- -PyLops supports computations on GPUs powered by `CuPy `_ (``cupy-cudaXX>=10.6.0``). +PyLops supports computations on GPUs powered by `CuPy `_ (``cupy-cudaXX>=v13.0.0``). This library must be installed *before* PyLops is installed. .. note:: diff --git a/pylops/signalprocessing/sliding1d.py b/pylops/signalprocessing/sliding1d.py index c9db1373..76ebc92e 100644 --- a/pylops/signalprocessing/sliding1d.py +++ b/pylops/signalprocessing/sliding1d.py @@ -7,7 +7,6 @@ from typing import Tuple, Union import numpy as np -from numpy.lib.stride_tricks import sliding_window_view from pylops import LinearOperator from pylops.signalprocessing.sliding2d import _slidingsteps @@ -193,10 +192,6 @@ def __init__( self.simOp = True self.Op = Op - # create temporary shape and strides for cpy - self.shape_wins = None - self.strides_wins = None - super().__init__( dtype=Op.dtype, dims=(nwins, int(dim[0] // nwins)), diff --git a/pylops/signalprocessing/sliding2d.py b/pylops/signalprocessing/sliding2d.py index c97802c6..71d7ebc9 100644 --- a/pylops/signalprocessing/sliding2d.py +++ b/pylops/signalprocessing/sliding2d.py @@ -9,7 +9,13 @@ import numpy as np from pylops import LinearOperator -from pylops.basicoperators import BlockDiag, Diagonal, HStack, Restriction +from pylops.utils._internal import _value_or_sized_to_tuple +from pylops.utils.backend import ( + get_array_module, + get_sliding_window_view, + to_cupy_conditional, +) +from pylops.utils.decorators import reshaped from pylops.utils.tapers import taper2d from pylops.utils.typing import InputDimsLike, NDArray @@ -110,15 +116,7 @@ def sliding2d_design( return nwins, dims, mwins_inends, dwins_inends -def Sliding2D( - Op: LinearOperator, - dims: InputDimsLike, - dimsd: InputDimsLike, - nwin: int, - nover: int, - tapertype: str = "hanning", - name: str = "S", -) -> LinearOperator: +class Sliding2D(LinearOperator): """2D Sliding transform operator. Apply a transform operator ``Op`` repeatedly to slices of the model @@ -139,6 +137,12 @@ def Sliding2D( ``nover``, it is recommended to first run ``sliding2d_design`` to obtain the corresponding ``dims`` and number of windows. + .. note:: Two kind of operators ``Op`` can be provided: the first + applies a single transformation to each window separately; the second + applies the transformation to all of the windows at the same time. This + is directly inferred during initialization when the following condition + holds ``Op.shape[1] == np.prod(dims)``. + .. warning:: Depending on the choice of `nwin` and `nover` as well as the size of the data, sliding windows may not cover the entire data. The start and end indices of each window will be displayed and returned @@ -176,47 +180,104 @@ def Sliding2D( shape (``dims``). """ - # data windows - dwin_ins, dwin_ends = _slidingsteps(dimsd[0], nwin, nover) - nwins = len(dwin_ins) - - # check patching - if nwins * Op.shape[1] // dims[1] != dims[0]: - raise ValueError( - f"Model shape (dims={dims}) is not consistent with chosen " - f"number of windows. Run sliding2d_design to identify the " - f"correct number of windows for the current " - "model size..." - ) - # create tapers - if tapertype is not None: - tap = taper2d(dimsd[1], nwin, nover, tapertype=tapertype).astype(Op.dtype) - tapin = tap.copy() - tapin[:nover] = 1 - tapend = tap.copy() - tapend[-nover:] = 1 - taps = {} - taps[0] = tapin - for i in range(1, nwins - 1): - taps[i] = tap - taps[nwins - 1] = tapend - - # transform to apply - if tapertype is None: - OOp = BlockDiag([Op for _ in range(nwins)]) - else: - OOp = BlockDiag( - [Diagonal(taps[itap].ravel(), dtype=Op.dtype) * Op for itap in range(nwins)] + def __init__( + self, + Op: LinearOperator, + dims: InputDimsLike, + dimsd: InputDimsLike, + nwin: int, + nover: int, + tapertype: str = "hanning", + name: str = "S", + ) -> None: + + dims: Tuple[int, ...] = _value_or_sized_to_tuple(dims) + dimsd: Tuple[int, ...] = _value_or_sized_to_tuple(dimsd) + + # data windows + dwin_ins, dwin_ends = _slidingsteps(dimsd[0], nwin, nover) + self.dwin_inends = (dwin_ins, dwin_ends) + nwins = len(dwin_ins) + self.nwin = nwin + self.nover = nover + + # check patching + if nwins * Op.shape[1] // dims[1] != dims[0] and Op.shape[1] != np.prod(dims): + raise ValueError( + f"Model shape (dims={dims}) is not consistent with chosen " + f"number of windows. Run sliding2d_design to identify the " + f"correct number of windows for the current " + "model size..." + ) + + # create tapers + self.tapertype = tapertype + if self.tapertype is not None: + tap = taper2d(dimsd[1], nwin, nover, tapertype=self.tapertype) + tapin = tap.copy() + tapin[:nover] = 1 + tapend = tap.copy() + tapend[-nover:] = 1 + self.taps = [ + tapin[np.newaxis, :], + ] + for i in range(1, nwins - 1): + self.taps.append(tap[np.newaxis, :]) + self.taps.append(tapend[np.newaxis, :]) + self.taps = np.concatenate(self.taps, axis=0) + + # check if operator is applied to all windows simultaneously + self.simOp = False + if Op.shape[1] == np.prod(dims): + self.simOp = True + self.Op = Op + + super().__init__( + dtype=Op.dtype, + dims=(nwins, int(dims[0] // nwins), dims[1]), + dimsd=dimsd, + clinear=False, + name=name, ) - combining = HStack( - [ - Restriction(dimsd, range(win_in, win_end), axis=0, dtype=Op.dtype).H - for win_in, win_end in zip(dwin_ins, dwin_ends) - ] - ) - Sop = LinearOperator(combining * OOp) - Sop.dims, Sop.dimsd = (nwins, int(dims[0] // nwins), dims[1]), dimsd - Sop.name = name - return Sop + @reshaped + def _matvec(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + if self.tapertype is not None: + self.taps = to_cupy_conditional(x, self.taps) + y = ncp.zeros(self.dimsd, dtype=self.dtype) + if self.simOp: + x = self.Op @ x + for iwin0 in range(self.dims[0]): + if self.simOp: + xx = x[iwin0].reshape(self.nwin, self.dimsd[-1]) + else: + xx = self.Op.matvec(x[iwin0].ravel()).reshape(self.nwin, self.dimsd[-1]) + if self.tapertype is not None: + xxwin = self.taps[iwin0] * xx + else: + xxwin = xx + y[self.dwin_inends[0][iwin0] : self.dwin_inends[1][iwin0]] += xxwin + return y + + @reshaped + def _rmatvec(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + ncp_sliding_window_view = get_sliding_window_view(x) + if self.tapertype is not None: + self.taps = to_cupy_conditional(x, self.taps) + ywins = ncp_sliding_window_view(x, self.nwin, axis=0)[ + :: self.nwin - self.nover + ].transpose(0, 2, 1) + if self.tapertype is not None: + ywins = ywins * self.taps + if self.simOp: + y = self.Op.H @ ywins + else: + y = ncp.zeros(self.dims, dtype=self.dtype) + for iwin0 in range(self.dims[0]): + y[iwin0] = self.Op.rmatvec(ywins[iwin0].ravel()).reshape( + self.dims[1], self.dims[2] + ) + return y diff --git a/pylops/signalprocessing/sliding3d.py b/pylops/signalprocessing/sliding3d.py index 7d21018c..12bae843 100644 --- a/pylops/signalprocessing/sliding3d.py +++ b/pylops/signalprocessing/sliding3d.py @@ -6,9 +6,17 @@ import logging from typing import Tuple +import numpy as np + from pylops import LinearOperator -from pylops.basicoperators import BlockDiag, Diagonal, HStack, Restriction from pylops.signalprocessing.sliding2d import _slidingsteps +from pylops.utils._internal import _value_or_sized_to_tuple +from pylops.utils.backend import ( + get_array_module, + get_sliding_window_view, + to_cupy_conditional, +) +from pylops.utils.decorators import reshaped from pylops.utils.tapers import taper3d from pylops.utils.typing import InputDimsLike, NDArray @@ -89,17 +97,7 @@ def sliding3d_design( return nwins, dims, mwins_inends, dwins_inends -def Sliding3D( - Op: LinearOperator, - dims: InputDimsLike, - dimsd: InputDimsLike, - nwin: Tuple[int, int], - nover: Tuple[int, int], - nop: Tuple[int, int, int], - tapertype: str = "hanning", - nproc: int = 1, - name: str = "P", -) -> LinearOperator: +class Sliding3D(LinearOperator): """3D Sliding transform operator.w Apply a transform operator ``Op`` repeatedly to patches of the model @@ -121,6 +119,12 @@ def Sliding3D( ``nover``, it is recommended to first run ``sliding3d_design`` to obtain the corresponding ``dims`` and number of windows. + .. note:: Two kind of operators ``Op`` can be provided: the first + applies a single transformation to each window separately; the second + applies the transformation to all of the windows at the same time. This + is directly inferred during initialization when the following condition + holds ``Op.shape[1] == np.prod(dims)``. + .. warning:: Depending on the choice of `nwin` and `nover` as well as the size of the data, sliding windows may not cover the entire data. The start and end indices of each window will be displayed and returned @@ -146,8 +150,8 @@ def Sliding3D( tapertype : :obj:`str`, optional Type of taper (``hanning``, ``cosine``, ``cosinesquare`` or ``None``) nproc : :obj:`int`, optional - Number of processes used to evaluate the N operators in parallel - using ``multiprocessing``. If ``nproc=1``, work in serial mode. + *Deprecated*, will be removed in v3.0.0. Simply kept for + back-compatibility with previous implementation name : :obj:`str`, optional .. versionadded:: 2.0.0 @@ -165,66 +169,167 @@ def Sliding3D( shape (``dims``). """ - # data windows - dwin0_ins, dwin0_ends = _slidingsteps(dimsd[0], nwin[0], nover[0]) - dwin1_ins, dwin1_ends = _slidingsteps(dimsd[1], nwin[1], nover[1]) - nwins0 = len(dwin0_ins) - nwins1 = len(dwin1_ins) - nwins = nwins0 * nwins1 - - # check windows - if nwins * Op.shape[1] // dims[2] != dims[0] * dims[1]: - raise ValueError( - f"Model shape (dims={dims}) is not consistent with chosen " - f"number of windows. Run sliding3d_design to identify the " - f"correct number of windows for the current " - "model size..." - ) - # create tapers - if tapertype is not None: - tap = taper3d(dimsd[2], nwin, nover, tapertype=tapertype).astype(Op.dtype) - - # transform to apply - if tapertype is None: - OOp = BlockDiag([Op for _ in range(nwins)], nproc=nproc) - else: - OOp = BlockDiag( - [Diagonal(tap.ravel(), dtype=Op.dtype) * Op for _ in range(nwins)], - nproc=nproc, + def __init__( + self, + Op: LinearOperator, + dims: InputDimsLike, + dimsd: InputDimsLike, + nwin: Tuple[int, int], + nover: Tuple[int, int], + nop: Tuple[int, int, int], + tapertype: str = "hanning", + nproc: int = 1, + name: str = "P", + ) -> None: + + dims: Tuple[int, ...] = _value_or_sized_to_tuple(dims) + dimsd: Tuple[int, ...] = _value_or_sized_to_tuple(dimsd) + + # data windows + dwin0_ins, dwin0_ends = _slidingsteps(dimsd[0], nwin[0], nover[0]) + dwin1_ins, dwin1_ends = _slidingsteps(dimsd[1], nwin[1], nover[1]) + self.dwins_inends = ((dwin0_ins, dwin0_ends), (dwin1_ins, dwin1_ends)) + nwins0 = len(dwin0_ins) + nwins1 = len(dwin1_ins) + nwins = nwins0 * nwins1 + self.nwin = nwin + self.nover = nover + + # model windows + mwin0_ins, mwin0_ends = _slidingsteps(dims[0], nop[0], 0) + mwin1_ins, mwin1_ends = _slidingsteps(dims[1], nop[1], 0) + self.mwins_inends = ((mwin0_ins, mwin0_ends), (mwin1_ins, mwin1_ends)) + + # check windows + if nwins * Op.shape[1] // dims[2] != dims[0] * dims[1] and Op.shape[ + 1 + ] != np.prod(dims): + raise ValueError( + f"Model shape (dims={dims}) is not consistent with chosen " + f"number of windows. Run sliding3d_design to identify the " + f"correct number of windows for the current " + "model size..." + ) + + # create tapers + self.tapertype = tapertype + if self.tapertype is not None: + tap = taper3d(dimsd[2], nwin, nover, tapertype=tapertype).astype(Op.dtype) + taps = [ + tap, + ] * nwins # {itap: tap for itap in range(nwins)} + + # topmost tapers + taptop = tap.copy() + taptop[: nover[0]] = tap[nwin[0] // 2] + for itap in range(0, nwins1): + taps[itap] = taptop + # bottommost tapers + tapbottom = tap.copy() + tapbottom[-nover[0] :] = tap[nwin[0] // 2] + for itap in range(nwins - nwins1, nwins): + taps[itap] = tapbottom + # leftmost tapers + tapleft = tap.copy() + tapleft[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis] + for itap in range(0, nwins, nwins1): + taps[itap] = tapleft + # rightmost tapers + tapright = tap.copy() + tapright[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis] + for itap in range(nwins1 - 1, nwins, nwins1): + taps[itap] = tapright + # lefttopcorner taper + taplefttop = tap.copy() + taplefttop[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis] + taplefttop[: nover[0]] = taplefttop[nwin[0] // 2] + taps[0] = taplefttop + # righttopcorner taper + taprighttop = tap.copy() + taprighttop[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis] + taprighttop[: nover[0]] = taprighttop[nwin[0] // 2] + taps[nwins1 - 1] = taprighttop + # leftbottomcorner taper + tapleftbottom = tap.copy() + tapleftbottom[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis] + tapleftbottom[-nover[0] :] = tapleftbottom[nwin[0] // 2] + taps[nwins - nwins1] = tapleftbottom + # rightbottomcorner taper + taprightbottom = tap.copy() + taprightbottom[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis] + taprightbottom[-nover[0] :] = taprightbottom[nwin[0] // 2] + taps[nwins - 1] = taprightbottom + self.taps = np.vstack(taps).reshape( + nwins0, nwins1, nwin[0], nwin[1], dimsd[2] + ) + + # check if operator is applied to all windows simultaneously + self.simOp = False + if Op.shape[1] == np.prod(dims): + self.simOp = True + self.Op = Op + + super().__init__( + dtype=Op.dtype, + dims=( + nwins0, + nwins1, + int(dims[0] // nwins0), + int(dims[1] // nwins1), + dims[2], + ), + dimsd=dimsd, + clinear=False, + name=name, ) - hstack = HStack( - [ - Restriction( - (nwin[0], dimsd[1], dimsd[2]), - range(win_in, win_end), - axis=1, - dtype=Op.dtype, - ).H - for win_in, win_end in zip(dwin1_ins, dwin1_ends) - ] - ) + @reshaped + def _matvec(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + if self.tapertype is not None: + self.taps = to_cupy_conditional(x, self.taps) + y = ncp.zeros(self.dimsd, dtype=self.dtype) + if self.simOp: + x = self.Op @ x + for iwin0 in range(self.dims[0]): + for iwin1 in range(self.dims[1]): + if self.simOp: + xx = x[iwin0, iwin1].reshape( + self.nwin[0], self.nwin[1], self.dimsd[-1] + ) + else: + xx = self.Op.matvec(x[iwin0, iwin1].ravel()).reshape( + self.nwin[0], self.nwin[1], self.dimsd[-1] + ) + if self.tapertype is not None: + xxwin = self.taps[iwin0, iwin1] * xx + else: + xxwin = xx + y[ + self.dwins_inends[0][0][iwin0] : self.dwins_inends[0][1][iwin0], + self.dwins_inends[1][0][iwin1] : self.dwins_inends[1][1][iwin1], + ] += xxwin + return y - combining1 = BlockDiag([hstack] * nwins0) - combining0 = HStack( - [ - Restriction( - dimsd, - range(win_in, win_end), - axis=0, - dtype=Op.dtype, - ).H - for win_in, win_end in zip(dwin0_ins, dwin0_ends) - ] - ) - Sop = LinearOperator(combining0 * combining1 * OOp) - Sop.dims, Sop.dimsd = ( - nwins0, - nwins1, - int(dims[0] // nwins0), - int(dims[1] // nwins1), - dims[2], - ), dimsd - Sop.name = name - return Sop + @reshaped + def _rmatvec(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + ncp_sliding_window_view = get_sliding_window_view(x) + if self.tapertype is not None: + self.taps = to_cupy_conditional(x, self.taps) + ywins = ncp_sliding_window_view(x, self.nwin, axis=(0, 1))[ + :: self.nwin[0] - self.nover[0], :: self.nwin[1] - self.nover[1] + ].transpose(0, 1, 3, 4, 2) + if self.tapertype is not None: + ywins = ywins * self.taps + if self.simOp: + y = self.Op.H @ ywins + else: + y = ncp.zeros(self.dims, dtype=self.dtype) + for iwin0 in range(self.dims[0]): + for iwin1 in range(self.dims[1]): + y[iwin0, iwin1] = self.Op.rmatvec( + ywins[iwin0, iwin1].ravel() + ).reshape(self.dims[2], self.dims[3], self.dims[4]) + return y From 9cabdf4d318fde7b0467694b8dba67422e09cf27 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sun, 12 May 2024 23:01:10 +0300 Subject: [PATCH 03/11] feature: improved Patch2D Changed patch2d from using other operators to being directly implemented (including also an option to apply Op simultaneously to all windows). --- pylops/signalprocessing/patch2d.py | 262 ++++++++++++++++----------- pylops/signalprocessing/sliding3d.py | 2 +- 2 files changed, 155 insertions(+), 109 deletions(-) diff --git a/pylops/signalprocessing/patch2d.py b/pylops/signalprocessing/patch2d.py index d95ddeb5..1d54e97a 100644 --- a/pylops/signalprocessing/patch2d.py +++ b/pylops/signalprocessing/patch2d.py @@ -9,8 +9,14 @@ import numpy as np from pylops import LinearOperator -from pylops.basicoperators import BlockDiag, Diagonal, HStack, Restriction from pylops.signalprocessing.sliding2d import _slidingsteps +from pylops.utils._internal import _value_or_sized_to_tuple +from pylops.utils.backend import ( + get_array_module, + get_sliding_window_view, + to_cupy_conditional, +) +from pylops.utils.decorators import reshaped from pylops.utils.tapers import taper2d from pylops.utils.typing import InputDimsLike, NDArray @@ -91,17 +97,7 @@ def patch2d_design( return nwins, dims, mwins_inends, dwins_inends -def Patch2D( - Op: LinearOperator, - dims: InputDimsLike, - dimsd: InputDimsLike, - nwin: Tuple[int, int], - nover: Tuple[int, int], - nop: Tuple[int, int], - tapertype: str = "hanning", - scalings: Optional[Sequence[float]] = None, - name: str = "P", -) -> LinearOperator: +class Patch2D(LinearOperator): """2D Patch transform operator. Apply a transform operator ``Op`` repeatedly to patches of the model @@ -172,104 +168,154 @@ def Patch2D( Patch3D: 3D Patching transform operator. """ - # data windows - dwin0_ins, dwin0_ends = _slidingsteps(dimsd[0], nwin[0], nover[0]) - dwin1_ins, dwin1_ends = _slidingsteps(dimsd[1], nwin[1], nover[1]) - nwins0 = len(dwin0_ins) - nwins1 = len(dwin1_ins) - nwins = nwins0 * nwins1 - - # check patching - if nwins0 * nop[0] != dims[0] or nwins1 * nop[1] != dims[1]: - raise ValueError( - f"Model shape (dims={dims}) is not consistent with chosen " - f"number of windows. Run patch2d_design to identify the " - f"correct number of windows for the current " - "model size..." - ) - # create tapers - if tapertype is not None: - tap = taper2d(nwin[1], nwin[0], nover, tapertype=tapertype).astype(Op.dtype) - taps = {itap: tap for itap in range(nwins)} - # topmost tapers - taptop = tap.copy() - taptop[: nover[0]] = tap[nwin[0] // 2] - for itap in range(0, nwins1): - taps[itap] = taptop - # bottommost tapers - tapbottom = tap.copy() - tapbottom[-nover[0] :] = tap[nwin[0] // 2] - for itap in range(nwins - nwins1, nwins): - taps[itap] = tapbottom - # leftmost tapers - tapleft = tap.copy() - tapleft[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis] - for itap in range(0, nwins, nwins1): - taps[itap] = tapleft - # rightmost tapers - tapright = tap.copy() - tapright[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis] - for itap in range(nwins1 - 1, nwins, nwins1): - taps[itap] = tapright - # lefttopcorner taper - taplefttop = tap.copy() - taplefttop[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis] - taplefttop[: nover[0]] = taplefttop[nwin[0] // 2] - taps[0] = taplefttop - # righttopcorner taper - taprighttop = tap.copy() - taprighttop[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis] - taprighttop[: nover[0]] = taprighttop[nwin[0] // 2] - taps[nwins1 - 1] = taprighttop - # leftbottomcorner taper - tapleftbottom = tap.copy() - tapleftbottom[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis] - tapleftbottom[-nover[0] :] = tapleftbottom[nwin[0] // 2] - taps[nwins - nwins1] = tapleftbottom - # rightbottomcorner taper - taprightbottom = tap.copy() - taprightbottom[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis] - taprightbottom[-nover[0] :] = taprightbottom[nwin[0] // 2] - taps[nwins - 1] = taprightbottom - - # define scalings - if scalings is None: - scalings = [1.0] * nwins - - # transform to apply - if tapertype is None: - OOp = BlockDiag([scalings[itap] * Op for itap in range(nwins)]) - else: - OOp = BlockDiag( - [ - scalings[itap] * Diagonal(taps[itap].ravel(), dtype=Op.dtype) * Op - for itap in range(nwins) - ] + def __init__( + self, + Op: LinearOperator, + dims: InputDimsLike, + dimsd: InputDimsLike, + nwin: Tuple[int, int], + nover: Tuple[int, int], + nop: Tuple[int, int], + tapertype: str = "hanning", + scalings: Optional[Sequence[float]] = None, + name: str = "P", + ) -> None: + + dims: Tuple[int, ...] = _value_or_sized_to_tuple(dims) + dimsd: Tuple[int, ...] = _value_or_sized_to_tuple(dimsd) + + # data windows + dwin0_ins, dwin0_ends = _slidingsteps(dimsd[0], nwin[0], nover[0]) + dwin1_ins, dwin1_ends = _slidingsteps(dimsd[1], nwin[1], nover[1]) + self.dwins_inends = ((dwin0_ins, dwin0_ends), (dwin1_ins, dwin1_ends)) + nwins0 = len(dwin0_ins) + nwins1 = len(dwin1_ins) + nwins = nwins0 * nwins1 + self.nwin = nwin + self.nover = nover + + # check patching + if nwins0 * nop[0] != dims[0] or nwins1 * nop[1] != dims[1]: + raise ValueError( + f"Model shape (dims={dims}) is not consistent with chosen " + f"number of windows. Run patch2d_design to identify the " + f"correct number of windows for the current " + "model size..." + ) + + # create tapers + self.tapertype = tapertype + if self.tapertype is not None: + tap = taper2d(nwin[1], nwin[0], nover, tapertype=tapertype).astype(Op.dtype) + taps = [ + tap, + ] * nwins + # topmost tapers + taptop = tap.copy() + taptop[: nover[0]] = tap[nwin[0] // 2] + for itap in range(0, nwins1): + taps[itap] = taptop + # bottommost tapers + tapbottom = tap.copy() + tapbottom[-nover[0] :] = tap[nwin[0] // 2] + for itap in range(nwins - nwins1, nwins): + taps[itap] = tapbottom + # leftmost tapers + tapleft = tap.copy() + tapleft[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis] + for itap in range(0, nwins, nwins1): + taps[itap] = tapleft + # rightmost tapers + tapright = tap.copy() + tapright[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis] + for itap in range(nwins1 - 1, nwins, nwins1): + taps[itap] = tapright + # lefttopcorner taper + taplefttop = tap.copy() + taplefttop[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis] + taplefttop[: nover[0]] = taplefttop[nwin[0] // 2] + taps[0] = taplefttop + # righttopcorner taper + taprighttop = tap.copy() + taprighttop[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis] + taprighttop[: nover[0]] = taprighttop[nwin[0] // 2] + taps[nwins1 - 1] = taprighttop + # leftbottomcorner taper + tapleftbottom = tap.copy() + tapleftbottom[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis] + tapleftbottom[-nover[0] :] = tapleftbottom[nwin[0] // 2] + taps[nwins - nwins1] = tapleftbottom + # rightbottomcorner taper + taprightbottom = tap.copy() + taprightbottom[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis] + taprightbottom[-nover[0] :] = taprightbottom[nwin[0] // 2] + taps[nwins - 1] = taprightbottom + self.taps = np.vstack(taps).reshape(nwins0, nwins1, nwin[0], nwin[1]) + + # define scalings + if scalings is None: + self.scalings = [1.0] * nwins + else: + self.scalings = scalings + + # check if operator is applied to all windows simultaneously + self.simOp = False + if Op.shape[1] == np.prod(dims): + self.simOp = True + self.Op = Op + + super().__init__( + dtype=Op.dtype, + dims=(nwins0, nwins1, int(dims[0] // nwins0), int(dims[1] // nwins1)), + dimsd=dimsd, + clinear=False, + name=name, ) - hstack = HStack( - [ - Restriction( - (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) - ] - ) - combining1 = BlockDiag([hstack] * nwins0) + @reshaped() + def _matvec(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + if self.tapertype is not None: + self.taps = to_cupy_conditional(x, self.taps) + y = ncp.zeros(self.dimsd, dtype=self.dtype) + if self.simOp: + x = self.Op @ x + for iwin0 in range(self.dims[0]): + for iwin1 in range(self.dims[1]): + if self.simOp: + xx = x[iwin0, iwin1].reshape(self.nwin) + else: + xx = self.Op.matvec(x[iwin0, iwin1].ravel()).reshape(self.nwin) + if self.tapertype is not None: + xxwin = self.taps[iwin0, iwin1] * xx + else: + xxwin = xx - combining0 = HStack( - [ - Restriction(dimsd, range(win_in, win_end), axis=0, dtype=Op.dtype).H - for win_in, win_end in zip(dwin0_ins, dwin0_ends) + y[ + self.dwins_inends[0][0][iwin0] : self.dwins_inends[0][1][iwin0], + self.dwins_inends[1][0][iwin1] : self.dwins_inends[1][1][iwin1], + ] += xxwin + return y + + @reshaped + def _rmatvec(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + ncp_sliding_window_view = get_sliding_window_view(x) + if self.tapertype is not None: + self.taps = to_cupy_conditional(x, self.taps) + ywins = ncp_sliding_window_view(x, self.nwin)[ + :: self.nwin[0] - self.nover[0], :: self.nwin[1] - self.nover[1] ] - ) - Pop = LinearOperator(combining0 * combining1 * OOp) - Pop.dims, Pop.dimsd = ( - nwins0, - nwins1, - int(dims[0] // nwins0), - int(dims[1] // nwins1), - ), dimsd - Pop.name = name - return Pop + if self.tapertype is not None: + ywins = ywins * self.taps + if self.simOp: + y = self.Op.H @ ywins + else: + y = ncp.zeros(self.dims, dtype=self.dtype) + for iwin0 in range(self.dims[0]): + for iwin1 in range(self.dims[1]): + y[iwin0, iwin1] = self.Op.rmatvec( + ywins[iwin0, iwin1].ravel() + ).reshape(self.dims[2], self.dims[3]) + return y diff --git a/pylops/signalprocessing/sliding3d.py b/pylops/signalprocessing/sliding3d.py index 12bae843..e583f171 100644 --- a/pylops/signalprocessing/sliding3d.py +++ b/pylops/signalprocessing/sliding3d.py @@ -218,7 +218,7 @@ def __init__( tap = taper3d(dimsd[2], nwin, nover, tapertype=tapertype).astype(Op.dtype) taps = [ tap, - ] * nwins # {itap: tap for itap in range(nwins)} + ] * nwins # topmost tapers taptop = tap.copy() From 63425e36f46a3b370ca3f7c45ddacbbbacc99b7f Mon Sep 17 00:00:00 2001 From: mrava87 Date: Mon, 13 May 2024 21:13:59 +0300 Subject: [PATCH 04/11] feature: improved Patch3D Changed patch3d from using other operators to being directly implemented (including also an option to apply Op simultaneously to all windows). --- pylops/signalprocessing/patch3d.py | 601 ++++++++++++++++------------- 1 file changed, 329 insertions(+), 272 deletions(-) diff --git a/pylops/signalprocessing/patch3d.py b/pylops/signalprocessing/patch3d.py index 011963cc..d7a87747 100644 --- a/pylops/signalprocessing/patch3d.py +++ b/pylops/signalprocessing/patch3d.py @@ -9,8 +9,14 @@ import numpy as np from pylops import LinearOperator -from pylops.basicoperators import BlockDiag, Diagonal, HStack, Restriction from pylops.signalprocessing.sliding2d import _slidingsteps +from pylops.utils._internal import _value_or_sized_to_tuple +from pylops.utils.backend import ( + get_array_module, + get_sliding_window_view, + to_cupy_conditional, +) +from pylops.utils.decorators import reshaped from pylops.utils.tapers import tapernd from pylops.utils.typing import InputDimsLike, NDArray @@ -106,17 +112,7 @@ def patch3d_design( return nwins, dims, mwins_inends, dwins_inends -def Patch3D( - Op, - dims: InputDimsLike, - dimsd: InputDimsLike, - nwin: Tuple[int, int, int], - nover: Tuple[int, int, int], - nop: Tuple[int, int, int], - tapertype: str = "hanning", - scalings: Optional[Sequence[float]] = None, - name: str = "P", -) -> LinearOperator: +class Patch3D(LinearOperator): """3D Patch transform operator. Apply a transform operator ``Op`` repeatedly to patches of the model @@ -185,272 +181,333 @@ def Patch3D( Patch2D: 2D Patching transform operator. """ - # data windows - dwin0_ins, dwin0_ends = _slidingsteps(dimsd[0], nwin[0], nover[0]) - dwin1_ins, dwin1_ends = _slidingsteps(dimsd[1], nwin[1], nover[1]) - dwin2_ins, dwin2_ends = _slidingsteps(dimsd[2], nwin[2], nover[2]) - nwins0 = len(dwin0_ins) - nwins1 = len(dwin1_ins) - nwins2 = len(dwin2_ins) - nwins = nwins0 * nwins1 * nwins2 - - # check patching - if ( - nwins0 * nop[0] != dims[0] - or nwins1 * nop[1] != dims[1] - or nwins2 * nop[2] != dims[2] - ): - raise ValueError( - f"Model shape (dims={dims}) is not consistent with chosen " - f"number of windows. Run patch3d_design to identify the " - f"correct number of windows for the current " - "model size..." + + def __init__( + self, + Op: LinearOperator, + dims: InputDimsLike, + dimsd: InputDimsLike, + nwin: Tuple[int, int, int], + nover: Tuple[int, int, int], + nop: Tuple[int, int, int], + tapertype: str = "hanning", + scalings: Optional[Sequence[float]] = None, + name: str = "P", + ) -> None: + + dims: Tuple[int, ...] = _value_or_sized_to_tuple(dims) + dimsd: Tuple[int, ...] = _value_or_sized_to_tuple(dimsd) + + # data windows + dwin0_ins, dwin0_ends = _slidingsteps(dimsd[0], nwin[0], nover[0]) + dwin1_ins, dwin1_ends = _slidingsteps(dimsd[1], nwin[1], nover[1]) + dwin2_ins, dwin2_ends = _slidingsteps(dimsd[2], nwin[2], nover[2]) + self.dwins_inends = ( + (dwin0_ins, dwin0_ends), + (dwin1_ins, dwin1_ends), + (dwin2_ins, dwin2_ends), ) + nwins0 = len(dwin0_ins) + nwins1 = len(dwin1_ins) + nwins2 = len(dwin2_ins) + nwins = nwins0 * nwins1 * nwins2 + self.nwin = nwin + self.nover = nover - # create tapers - if tapertype is not None: - tap = tapernd(nwin, nover, tapertype=tapertype).astype(Op.dtype) - taps = {itap: tap for itap in range(nwins)} - # 1, sides - # topmost tapers - taptop = tap.copy() - taptop[: nover[0]] = tap[nwin[0] // 2] - for itap in range(0, nwins1 * nwins2): - taps[itap] = taptop - # bottommost tapers - tapbottom = tap.copy() - tapbottom[-nover[0] :] = tap[nwin[0] // 2] - for itap in range(nwins - nwins1 * nwins2, nwins): - taps[itap] = tapbottom - # frontmost tapers - tapfront = tap.copy() - tapfront[:, :, : nover[2]] = tap[:, :, nwin[2] // 2][:, :, np.newaxis] - for itap in range(0, nwins, nwins2): - taps[itap] = tapfront - # backmost tapers - tapback = tap.copy() - tapback[:, :, -nover[2] :] = tap[:, :, nwin[2] // 2][:, :, np.newaxis] - for itap in range(nwins2 - 1, nwins, nwins2): - taps[itap] = tapback - # leftmost tapers - tapleft = tap.copy() - tapleft[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis, :] - for itap in range(0, nwins, nwins1 * nwins2): - for i in range(nwins2): - taps[itap + i] = tapleft - # rightmost tapers - tapright = tap.copy() - tapright[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis, :] - for itap in range(nwins2 * (nwins1 - 1), nwins, nwins2 * nwins1): - for i in range(nwins2): - taps[itap + i] = tapright - # 2. pillars - # topleftmost tapers - taplefttop = tap.copy() - taplefttop[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis, :] - taplefttop[: nover[0]] = taplefttop[nwin[0] // 2] - for itap in range(nwins2): - taps[itap] = taplefttop - # toprightmost tapers - taprighttop = tap.copy() - taprighttop[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis, :] - taprighttop[: nover[0]] = taprighttop[nwin[0] // 2] - for itap in range(nwins2): - taps[nwins2 * (nwins1 - 1) + itap] = taprighttop - # topfrontmost tapers - tapfronttop = tap.copy() - tapfronttop[:, :, : nover[2]] = tap[:, :, nwin[2] // 2][:, :, np.newaxis] - tapfronttop[: nover[0]] = tapfronttop[nwin[0] // 2] - for itap in range(0, nwins1 * nwins2, nwins2): - taps[itap] = tapfronttop - # topbackmost tapers - tapbacktop = tap.copy() - tapbacktop[:, :, -nover[2] :] = tap[:, :, nwin[2] // 2][:, :, np.newaxis] - tapbacktop[: nover[0]] = tapbacktop[nwin[0] // 2] - for itap in range(nwins2 - 1, nwins1 * nwins2, nwins2): - taps[itap] = tapbacktop - # bottomleftmost tapers - tapleftbottom = tap.copy() - tapleftbottom[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis, :] - tapleftbottom[-nover[0] :] = tapleftbottom[nwin[0] // 2] - for itap in range(nwins2): - taps[(nwins0 - 1) * nwins1 * nwins2 + itap] = tapleftbottom - # bottomrightmost tapers - taprightbottom = tap.copy() - taprightbottom[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis, :] - taprightbottom[-nover[0] :] = taprightbottom[nwin[0] // 2] - for itap in range(nwins2): + # check patching + if ( + nwins0 * nop[0] != dims[0] + or nwins1 * nop[1] != dims[1] + or nwins2 * nop[2] != dims[2] + ): + raise ValueError( + f"Model shape (dims={dims}) is not consistent with chosen " + f"number of windows. Run patch3d_design to identify the " + f"correct number of windows for the current " + "model size..." + ) + + # create tapers + self.tapertype = tapertype + if tapertype is not None: + tap = tapernd(nwin, nover, tapertype=tapertype).astype(Op.dtype) + taps = [ + tap, + ] * nwins + # 1, sides + # topmost tapers + taptop = tap.copy() + taptop[: nover[0]] = tap[nwin[0] // 2] + for itap in range(0, nwins1 * nwins2): + taps[itap] = taptop + # bottommost tapers + tapbottom = tap.copy() + tapbottom[-nover[0] :] = tap[nwin[0] // 2] + for itap in range(nwins - nwins1 * nwins2, nwins): + taps[itap] = tapbottom + # frontmost tapers + tapfront = tap.copy() + tapfront[:, :, : nover[2]] = tap[:, :, nwin[2] // 2][:, :, np.newaxis] + for itap in range(0, nwins, nwins2): + taps[itap] = tapfront + # backmost tapers + tapback = tap.copy() + tapback[:, :, -nover[2] :] = tap[:, :, nwin[2] // 2][:, :, np.newaxis] + for itap in range(nwins2 - 1, nwins, nwins2): + taps[itap] = tapback + # leftmost tapers + tapleft = tap.copy() + tapleft[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis, :] + for itap in range(0, nwins, nwins1 * nwins2): + for i in range(nwins2): + taps[itap + i] = tapleft + # rightmost tapers + tapright = tap.copy() + tapright[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis, :] + for itap in range(nwins2 * (nwins1 - 1), nwins, nwins2 * nwins1): + for i in range(nwins2): + taps[itap + i] = tapright + # 2. pillars + # topleftmost tapers + taplefttop = tap.copy() + taplefttop[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis, :] + taplefttop[: nover[0]] = taplefttop[nwin[0] // 2] + for itap in range(nwins2): + taps[itap] = taplefttop + # toprightmost tapers + taprighttop = tap.copy() + taprighttop[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis, :] + taprighttop[: nover[0]] = taprighttop[nwin[0] // 2] + for itap in range(nwins2): + taps[nwins2 * (nwins1 - 1) + itap] = taprighttop + # topfrontmost tapers + tapfronttop = tap.copy() + tapfronttop[:, :, : nover[2]] = tap[:, :, nwin[2] // 2][:, :, np.newaxis] + tapfronttop[: nover[0]] = tapfronttop[nwin[0] // 2] + for itap in range(0, nwins1 * nwins2, nwins2): + taps[itap] = tapfronttop + # topbackmost tapers + tapbacktop = tap.copy() + tapbacktop[:, :, -nover[2] :] = tap[:, :, nwin[2] // 2][:, :, np.newaxis] + tapbacktop[: nover[0]] = tapbacktop[nwin[0] // 2] + for itap in range(nwins2 - 1, nwins1 * nwins2, nwins2): + taps[itap] = tapbacktop + # bottomleftmost tapers + tapleftbottom = tap.copy() + tapleftbottom[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis, :] + tapleftbottom[-nover[0] :] = tapleftbottom[nwin[0] // 2] + for itap in range(nwins2): + taps[(nwins0 - 1) * nwins1 * nwins2 + itap] = tapleftbottom + # bottomrightmost tapers + taprightbottom = tap.copy() + taprightbottom[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis, :] + taprightbottom[-nover[0] :] = taprightbottom[nwin[0] // 2] + for itap in range(nwins2): + taps[ + (nwins0 - 1) * nwins1 * nwins2 + (nwins1 - 1) * nwins2 + itap + ] = taprightbottom + # bottomfrontmost tapers + tapfrontbottom = tap.copy() + tapfrontbottom[:, :, : nover[2]] = tap[:, :, nwin[2] // 2][:, :, np.newaxis] + tapfrontbottom[-nover[0] :] = tapfrontbottom[nwin[0] // 2] + for itap in range(0, nwins1 * nwins2, nwins2): + taps[(nwins0 - 1) * nwins1 * nwins2 + itap] = tapfrontbottom + # bottombackmost tapers + tapbackbottom = tap.copy() + tapbackbottom[:, :, -nover[2] :] = tap[:, :, nwin[2] // 2][:, :, np.newaxis] + tapbackbottom[-nover[0] :] = tapbackbottom[nwin[0] // 2] + for itap in range(0, nwins1 * nwins2, nwins2): + taps[(nwins0 - 1) * nwins1 * nwins2 + nwins2 + itap - 1] = tapbackbottom + # leftfrontmost tapers + tapleftfront = tap.copy() + tapleftfront[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis, :] + tapleftfront[:, :, : nover[2]] = tapleftfront[:, :, nwin[2] // 2][ + :, :, np.newaxis + ] + for itap in range(0, nwins, nwins1 * nwins2): + taps[itap] = tapleftfront + # rightfrontmost tapers + taprightfront = tap.copy() + taprightfront[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis, :] + taprightfront[:, :, : nover[2]] = taprightfront[:, :, nwin[2] // 2][ + :, :, np.newaxis + ] + for itap in range(0, nwins, nwins1 * nwins2): + taps[(nwins1 - 1) * nwins2 + itap] = taprightfront + # leftbackmost tapers + tapleftback = tap.copy() + tapleftback[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis, :] + tapleftback[:, :, -nover[2] :] = tapleftback[:, :, nwin[2] // 2][ + :, :, np.newaxis + ] + for itap in range(0, nwins, nwins1 * nwins2): + taps[nwins2 + itap - 1] = tapleftback + # rightbackmost tapers + taprightback = tap.copy() + taprightback[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis, :] + taprightback[:, :, -nover[2] :] = taprightback[:, :, nwin[2] // 2][ + :, :, np.newaxis + ] + for itap in range(0, nwins, nwins1 * nwins2): + taps[(nwins1 - 1) * nwins2 + nwins2 + itap - 1] = taprightback + # 3. corners + # lefttopfrontcorner taper + taplefttop = tap.copy() + taplefttop[: nover[0]] = tap[nwin[0] // 2] + taplefttop[:, : nover[1]] = taplefttop[:, nwin[1] // 2][:, np.newaxis, :] + taplefttop[:, :, : nover[2]] = taplefttop[:, :, nwin[2] // 2][ + :, :, np.newaxis + ] + taps[0] = taplefttop + # lefttopbackcorner taper + taplefttop = tap.copy() + taplefttop[: nover[0]] = tap[nwin[0] // 2] + taplefttop[:, : nover[1]] = taplefttop[:, nwin[1] // 2][:, np.newaxis, :] + taplefttop[:, :, -nover[2] :] = taplefttop[:, :, nwin[2] // 2][ + :, :, np.newaxis + ] + taps[nwins2 - 1] = taplefttop + # righttopfrontcorner taper + taprighttop = tap.copy() + taprighttop[: nover[0]] = tap[nwin[0] // 2] + taprighttop[:, -nover[1] :] = taprighttop[:, nwin[1] // 2][:, np.newaxis, :] + taprighttop[:, :, : nover[2]] = taprighttop[:, :, nwin[2] // 2][ + :, :, np.newaxis + ] + taps[(nwins1 - 1) * nwins2] = taprighttop + # righttopbackcorner taper + taprighttop = tap.copy() + taprighttop[: nover[0]] = tap[nwin[0] // 2] + taprighttop[:, -nover[1] :] = taprighttop[:, nwin[1] // 2][:, np.newaxis, :] + taprighttop[:, :, -nover[2] :] = taprighttop[:, :, nwin[2] // 2][ + :, :, np.newaxis + ] + taps[(nwins1 - 1) * nwins2 + nwins2 - 1] = taprighttop + # leftbottomfrontcorner taper + tapleftbottom = tap.copy() + tapleftbottom[-nover[0] :] = tap[nwin[0] // 2] + tapleftbottom[:, : nover[1]] = tapleftbottom[:, nwin[1] // 2][ + :, np.newaxis, : + ] + tapleftbottom[:, :, : nover[2]] = tapleftbottom[:, :, nwin[2] // 2][ + :, :, np.newaxis + ] + taps[(nwins0 - 1) * nwins1 * nwins2] = tapleftbottom + # leftbottombackcorner taper + tapleftbottom = tap.copy() + tapleftbottom[-nover[0] :] = tap[nwin[0] // 2] + tapleftbottom[:, : nover[1]] = tapleftbottom[:, nwin[1] // 2][ + :, np.newaxis, : + ] + tapleftbottom[:, :, -nover[2] :] = tapleftbottom[:, :, nwin[2] // 2][ + :, :, np.newaxis + ] + taps[(nwins0 - 1) * nwins1 * nwins2 + nwins2 - 1] = tapleftbottom + # rightbottomfrontcorner taper + taprightbottom = tap.copy() + taprightbottom[-nover[0] :] = tap[nwin[0] // 2] + taprightbottom[:, -nover[1] :] = taprightbottom[:, nwin[1] // 2][ + :, np.newaxis, : + ] + taprightbottom[:, :, : nover[2]] = taprightbottom[:, :, nwin[2] // 2][ + :, :, np.newaxis + ] taps[ - (nwins0 - 1) * nwins1 * nwins2 + (nwins1 - 1) * nwins2 + itap + (nwins0 - 1) * nwins1 * nwins2 + (nwins1 - 1) * nwins2 ] = taprightbottom - # bottomfrontmost tapers - tapfrontbottom = tap.copy() - tapfrontbottom[:, :, : nover[2]] = tap[:, :, nwin[2] // 2][:, :, np.newaxis] - tapfrontbottom[-nover[0] :] = tapfrontbottom[nwin[0] // 2] - for itap in range(0, nwins1 * nwins2, nwins2): - taps[(nwins0 - 1) * nwins1 * nwins2 + itap] = tapfrontbottom - # bottombackmost tapers - tapbackbottom = tap.copy() - tapbackbottom[:, :, -nover[2] :] = tap[:, :, nwin[2] // 2][:, :, np.newaxis] - tapbackbottom[-nover[0] :] = tapbackbottom[nwin[0] // 2] - for itap in range(0, nwins1 * nwins2, nwins2): - taps[(nwins0 - 1) * nwins1 * nwins2 + nwins2 + itap - 1] = tapbackbottom - # leftfrontmost tapers - tapleftfront = tap.copy() - tapleftfront[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis, :] - tapleftfront[:, :, : nover[2]] = tapleftfront[:, :, nwin[2] // 2][ - :, :, np.newaxis - ] - for itap in range(0, nwins, nwins1 * nwins2): - taps[itap] = tapleftfront - # rightfrontmost tapers - taprightfront = tap.copy() - taprightfront[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis, :] - taprightfront[:, :, : nover[2]] = taprightfront[:, :, nwin[2] // 2][ - :, :, np.newaxis - ] - for itap in range(0, nwins, nwins1 * nwins2): - taps[(nwins1 - 1) * nwins2 + itap] = taprightfront - # leftbackmost tapers - tapleftback = tap.copy() - tapleftback[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis, :] - tapleftback[:, :, -nover[2] :] = tapleftback[:, :, nwin[2] // 2][ - :, :, np.newaxis - ] - for itap in range(0, nwins, nwins1 * nwins2): - taps[nwins2 + itap - 1] = tapleftback - # rightbackmost tapers - taprightback = tap.copy() - taprightback[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis, :] - taprightback[:, :, -nover[2] :] = taprightback[:, :, nwin[2] // 2][ - :, :, np.newaxis - ] - for itap in range(0, nwins, nwins1 * nwins2): - taps[(nwins1 - 1) * nwins2 + nwins2 + itap - 1] = taprightback - # 3. corners - # lefttopfrontcorner taper - taplefttop = tap.copy() - taplefttop[: nover[0]] = tap[nwin[0] // 2] - taplefttop[:, : nover[1]] = taplefttop[:, nwin[1] // 2][:, np.newaxis, :] - taplefttop[:, :, : nover[2]] = taplefttop[:, :, nwin[2] // 2][:, :, np.newaxis] - taps[0] = taplefttop - # lefttopbackcorner taper - taplefttop = tap.copy() - taplefttop[: nover[0]] = tap[nwin[0] // 2] - taplefttop[:, : nover[1]] = taplefttop[:, nwin[1] // 2][:, np.newaxis, :] - taplefttop[:, :, -nover[2] :] = taplefttop[:, :, nwin[2] // 2][:, :, np.newaxis] - taps[nwins2 - 1] = taplefttop - # righttopfrontcorner taper - taprighttop = tap.copy() - taprighttop[: nover[0]] = tap[nwin[0] // 2] - taprighttop[:, -nover[1] :] = taprighttop[:, nwin[1] // 2][:, np.newaxis, :] - taprighttop[:, :, : nover[2]] = taprighttop[:, :, nwin[2] // 2][ - :, :, np.newaxis - ] - taps[(nwins1 - 1) * nwins2] = taprighttop - # righttopbackcorner taper - taprighttop = tap.copy() - taprighttop[: nover[0]] = tap[nwin[0] // 2] - taprighttop[:, -nover[1] :] = taprighttop[:, nwin[1] // 2][:, np.newaxis, :] - taprighttop[:, :, -nover[2] :] = taprighttop[:, :, nwin[2] // 2][ - :, :, np.newaxis - ] - taps[(nwins1 - 1) * nwins2 + nwins2 - 1] = taprighttop - # leftbottomfrontcorner taper - tapleftbottom = tap.copy() - tapleftbottom[-nover[0] :] = tap[nwin[0] // 2] - tapleftbottom[:, : nover[1]] = tapleftbottom[:, nwin[1] // 2][:, np.newaxis, :] - tapleftbottom[:, :, : nover[2]] = tapleftbottom[:, :, nwin[2] // 2][ - :, :, np.newaxis - ] - taps[(nwins0 - 1) * nwins1 * nwins2] = tapleftbottom - # leftbottombackcorner taper - tapleftbottom = tap.copy() - tapleftbottom[-nover[0] :] = tap[nwin[0] // 2] - tapleftbottom[:, : nover[1]] = tapleftbottom[:, nwin[1] // 2][:, np.newaxis, :] - tapleftbottom[:, :, -nover[2] :] = tapleftbottom[:, :, nwin[2] // 2][ - :, :, np.newaxis - ] - taps[(nwins0 - 1) * nwins1 * nwins2 + nwins2 - 1] = tapleftbottom - # rightbottomfrontcorner taper - taprightbottom = tap.copy() - taprightbottom[-nover[0] :] = tap[nwin[0] // 2] - taprightbottom[:, -nover[1] :] = taprightbottom[:, nwin[1] // 2][ - :, np.newaxis, : - ] - taprightbottom[:, :, : nover[2]] = taprightbottom[:, :, nwin[2] // 2][ - :, :, np.newaxis - ] - taps[(nwins0 - 1) * nwins1 * nwins2 + (nwins1 - 1) * nwins2] = taprightbottom - # rightbottombackcorner taper - taprightbottom = tap.copy() - taprightbottom[-nover[0] :] = tap[nwin[0] // 2] - taprightbottom[:, -nover[1] :] = taprightbottom[:, nwin[1] // 2][ - :, np.newaxis, : - ] - taprightbottom[:, :, -nover[2] :] = taprightbottom[:, :, nwin[2] // 2][ - :, :, np.newaxis - ] - taps[ - (nwins0 - 1) * nwins1 * nwins2 + (nwins1 - 1) * nwins2 + nwins2 - 1 - ] = taprightbottom - - # define scalings - if scalings is None: - scalings = [1.0] * nwins - - # transform to apply - if tapertype is None: - OOp = BlockDiag([scalings[itap] * Op for itap in range(nwins)]) - else: - OOp = BlockDiag( - [ - scalings[itap] * Diagonal(taps[itap].ravel(), dtype=Op.dtype) * Op - for itap in range(nwins) + # rightbottombackcorner taper + taprightbottom = tap.copy() + taprightbottom[-nover[0] :] = tap[nwin[0] // 2] + taprightbottom[:, -nover[1] :] = taprightbottom[:, nwin[1] // 2][ + :, np.newaxis, : ] + taprightbottom[:, :, -nover[2] :] = taprightbottom[:, :, nwin[2] // 2][ + :, :, np.newaxis + ] + taps[ + (nwins0 - 1) * nwins1 * nwins2 + (nwins1 - 1) * nwins2 + nwins2 - 1 + ] = taprightbottom + self.taps = np.vstack(taps).reshape( + nwins0, nwins1, nwins2, nwin[0], nwin[1], nwin[2] + ) + + # define scalings + if scalings is None: + self.scalings = [1.0] * nwins + else: + self.scalings = scalings + + # check if operator is applied to all windows simultaneously + self.simOp = False + if Op.shape[1] == np.prod(dims): + self.simOp = True + self.Op = Op + + super().__init__( + dtype=Op.dtype, + dims=( + nwins0, + nwins1, + nwins2, + int(dims[0] // nwins0), + int(dims[1] // nwins1), + int(dims[2] // nwins2), + ), + dimsd=dimsd, + clinear=False, + name=name, ) - hstack2 = HStack( - [ - Restriction( - (nwin[0], nwin[1], dimsd[2]), - range(win_in, win_end), - axis=2, - dtype=Op.dtype, - ).H - for win_in, win_end in zip(dwin2_ins, dwin2_ends) - ] - ) - combining2 = BlockDiag([hstack2] * (nwins1 * nwins0)) - - hstack1 = HStack( - [ - Restriction( - (nwin[0], dimsd[1], dimsd[2]), - range(win_in, win_end), - axis=1, - dtype=Op.dtype, - ).H - for win_in, win_end in zip(dwin1_ins, dwin1_ends) - ] - ) - combining1 = BlockDiag([hstack1] * nwins0) + @reshaped() + def _matvec(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + if self.tapertype is not None: + self.taps = to_cupy_conditional(x, self.taps) + y = ncp.zeros(self.dimsd, dtype=self.dtype) + if self.simOp: + x = self.Op @ x + for iwin0 in range(self.dims[0]): + for iwin1 in range(self.dims[1]): + for iwin2 in range(self.dims[2]): + if self.simOp: + xx = x[iwin0, iwin1, iwin2].reshape(self.nwin) + else: + xx = self.Op.matvec(x[iwin0, iwin1, iwin2].ravel()).reshape( + self.nwin + ) + if self.tapertype is not None: + xxwin = self.taps[iwin0, iwin1, iwin2] * xx + else: + xxwin = xx - combining0 = HStack( - [ - Restriction(dimsd, range(win_in, win_end), axis=0, dtype=Op.dtype).H - for win_in, win_end in zip(dwin0_ins, dwin0_ends) - ] - ) + y[ + self.dwins_inends[0][0][iwin0] : self.dwins_inends[0][1][iwin0], + self.dwins_inends[1][0][iwin1] : self.dwins_inends[1][1][iwin1], + self.dwins_inends[2][0][iwin2] : self.dwins_inends[2][1][iwin2], + ] += xxwin + return y - Pop = LinearOperator(combining0 * combining1 * combining2 * OOp) - Pop.dims, Pop.dimsd = ( - nwins0, - nwins1, - nwins2, - int(dims[0] // nwins0), - int(dims[1] // nwins1), - int(dims[2] // nwins2), - ), dimsd - Pop.name = name - return Pop + @reshaped + def _rmatvec(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + ncp_sliding_window_view = get_sliding_window_view(x) + if self.tapertype is not None: + self.taps = to_cupy_conditional(x, self.taps) + ywins = ncp_sliding_window_view(x, self.nwin)[ + :: self.nwin[0] - self.nover[0], + :: self.nwin[1] - self.nover[1], + :: self.nwin[2] - self.nover[2], + ] + if self.tapertype is not None: + ywins = ywins * self.taps + if self.simOp: + y = self.Op.H @ ywins + else: + y = ncp.zeros(self.dims, dtype=self.dtype) + for iwin0 in range(self.dims[0]): + for iwin1 in range(self.dims[1]): + for iwin2 in range(self.dims[2]): + y[iwin0, iwin1, iwin2] = self.Op.rmatvec( + ywins[iwin0, iwin1, iwin2].ravel() + ).reshape(self.dims[3], self.dims[4], self.dims[5]) + return y From 94115807a4509481afa83bdfd18b228200c8325b Mon Sep 17 00:00:00 2001 From: mrava87 Date: Wed, 15 May 2024 22:40:08 +0300 Subject: [PATCH 05/11] feat: added savetaper to Sliding1d and Sliding2d --- pylops/signalprocessing/sliding1d.py | 105 ++++++++++++++++++++++---- pylops/signalprocessing/sliding2d.py | 107 ++++++++++++++++++++++++--- pytests/test_sliding.py | 42 ++++++++++- 3 files changed, 227 insertions(+), 27 deletions(-) diff --git a/pylops/signalprocessing/sliding1d.py b/pylops/signalprocessing/sliding1d.py index 76ebc92e..8b3bb5d1 100644 --- a/pylops/signalprocessing/sliding1d.py +++ b/pylops/signalprocessing/sliding1d.py @@ -127,6 +127,9 @@ class Sliding1D(LinearOperator): Number of samples of overlapping part of window tapertype : :obj:`str`, optional Type of taper (``hanning``, ``cosine``, ``cosinesquare`` or ``None``) + savetaper: :obj:`bool`, optional + Save all tapers and apply them in one go (``True``) or save unique tapers and apply them one by one (``False``). + The first option is more computationally efficient, whilst the second is more memory efficient. name : :obj:`str`, optional .. versionadded:: 2.0.0 @@ -148,6 +151,7 @@ def __init__( nwin: int, nover: int, tapertype: str = "hanning", + savetaper: bool = True, name: str = "S", ) -> None: @@ -172,19 +176,23 @@ def __init__( # create tapers self.tapertype = tapertype + self.savetaper = savetaper if self.tapertype is not None: tap = taper(nwin, nover, tapertype=self.tapertype) tapin = tap.copy() tapin[:nover] = 1 tapend = tap.copy() tapend[-nover:] = 1 - self.taps = [ - tapin, - ] - for i in range(1, nwins - 1): - self.taps.append(tap) - self.taps.append(tapend) - self.taps = np.vstack(self.taps) + if self.savetaper: + self.taps = [ + tapin, + ] + for i in range(1, nwins - 1): + self.taps.append(tap) + self.taps.append(tapend) + self.taps = np.vstack(self.taps) + else: + self.taps = np.vstack([tapin, tap, tapend]) # check if operator is applied to all windows simultaneously self.simOp = False @@ -200,28 +208,30 @@ def __init__( name=name, ) + self._register_multiplications(self.savetaper) + @reshaped - def _matvec(self, x: NDArray) -> NDArray: + def _matvec_savetaper(self, x: NDArray) -> NDArray: ncp = get_array_module(x) if self.tapertype is not None: self.taps = to_cupy_conditional(x, self.taps) y = ncp.zeros(self.dimsd, dtype=self.dtype) if self.simOp: x = self.Op @ x + if self.tapertype is not None: + x = self.taps * x for iwin0 in range(self.dims[0]): if self.simOp: - xx = x[iwin0] - else: - xx = self.Op.matvec(x[iwin0]) - if self.tapertype is not None: - xxwin = self.taps[iwin0] * xx + xxwin = x[iwin0] else: - xxwin = xx + xxwin = self.Op.matvec(x[iwin0]) + if self.tapertype is not None: + xxwin = self.taps[iwin0] * xxwin y[self.dwin_inends[0][iwin0] : self.dwin_inends[1][iwin0]] += xxwin return y @reshaped - def _rmatvec(self, x: NDArray) -> NDArray: + def _rmatvec_savetaper(self, x: NDArray) -> NDArray: ncp = get_array_module(x) ncp_sliding_window_view = get_sliding_window_view(x) if self.tapertype is not None: @@ -236,3 +246,68 @@ def _rmatvec(self, x: NDArray) -> NDArray: for iwin0 in range(self.dims[0]): y[iwin0] = self.Op.rmatvec(ywins[iwin0]) return y + + @reshaped + def _matvec_nosavetaper(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + if self.tapertype is not None: + self.taps = to_cupy_conditional(x, self.taps) + y = ncp.zeros(self.dimsd, dtype=self.dtype) + if self.simOp: + x = self.Op @ x + for iwin0 in range(self.dims[0]): + if self.simOp: + xxwin = x[iwin0] + else: + xxwin = self.Op.matvec(x[iwin0]) + if self.tapertype is not None: + if iwin0 == 0: + xxwin = self.taps[0] * xxwin + elif iwin0 == self.dims[0] - 1: + xxwin = self.taps[-1] * xxwin + else: + xxwin = self.taps[1] * xxwin + y[self.dwin_inends[0][iwin0] : self.dwin_inends[1][iwin0]] += xxwin + return y + + @reshaped + def _rmatvec_nosavetaper(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + ncp_sliding_window_view = get_sliding_window_view(x) + if self.tapertype is not None: + self.taps = to_cupy_conditional(x, self.taps) + ywins = ncp_sliding_window_view(x, self.nwin)[:: self.nwin - self.nover].copy() + if self.simOp: + if self.tapertype is not None: + for iwin0 in range(self.dims[0]): + if iwin0 == 0: + ywins[0] = ywins[0] * self.taps[0] + elif iwin0 == self.dims[0] - 1: + ywins[-1] = ywins[-1] * self.taps[-1] + else: + ywins[iwin0] = ywins[iwin0] * self.taps[1] + y = self.Op.H @ ywins + else: + y = ncp.zeros(self.dims, dtype=self.dtype) + for iwin0 in range(self.dims[0]): + if iwin0 == 0: + if self.tapertype is not None: + ywins[0] = ywins[0] * self.taps[0] + y[0] = self.Op.rmatvec(ywins[0]) + elif iwin0 == self.dims[0] - 1: + if self.tapertype is not None: + ywins[-1] = ywins[-1] * self.taps[-1] + y[-1] = self.Op.rmatvec(ywins[-1]) + else: + if self.tapertype is not None: + ywins[iwin0] = ywins[iwin0] * self.taps[1] + y[iwin0] = self.Op.rmatvec(ywins[iwin0]) + return y + + def _register_multiplications(self, savetaper: bool) -> None: + if savetaper: + self._matvec = self._matvec_savetaper + self._rmatvec = self._rmatvec_savetaper + else: + self._matvec = self._matvec_nosavetaper + self._rmatvec = self._rmatvec_nosavetaper diff --git a/pylops/signalprocessing/sliding2d.py b/pylops/signalprocessing/sliding2d.py index 71d7ebc9..22268471 100644 --- a/pylops/signalprocessing/sliding2d.py +++ b/pylops/signalprocessing/sliding2d.py @@ -163,6 +163,9 @@ class Sliding2D(LinearOperator): Number of samples of overlapping part of window tapertype : :obj:`str`, optional Type of taper (``hanning``, ``cosine``, ``cosinesquare`` or ``None``) + savetaper: :obj:`bool`, optional + Save all tapers and apply them in one go (``True``) or save unique tapers and apply them one by one (``False``). + The first option is more computationally efficient, whilst the second is more memory efficient. name : :obj:`str`, optional .. versionadded:: 2.0.0 @@ -189,6 +192,7 @@ def __init__( nwin: int, nover: int, tapertype: str = "hanning", + savetaper: bool = True, name: str = "S", ) -> None: @@ -213,19 +217,25 @@ def __init__( # create tapers self.tapertype = tapertype + self.savetaper = savetaper if self.tapertype is not None: tap = taper2d(dimsd[1], nwin, nover, tapertype=self.tapertype) tapin = tap.copy() tapin[:nover] = 1 tapend = tap.copy() tapend[-nover:] = 1 - self.taps = [ - tapin[np.newaxis, :], - ] - for i in range(1, nwins - 1): - self.taps.append(tap[np.newaxis, :]) - self.taps.append(tapend[np.newaxis, :]) - self.taps = np.concatenate(self.taps, axis=0) + if self.savetaper: + self.taps = [ + tapin[np.newaxis, :], + ] + for i in range(1, nwins - 1): + self.taps.append(tap[np.newaxis, :]) + self.taps.append(tapend[np.newaxis, :]) + self.taps = np.concatenate(self.taps, axis=0) + else: + self.taps = np.vstack( + [tapin[np.newaxis, :], tap[np.newaxis, :], tapend[np.newaxis, :]] + ) # check if operator is applied to all windows simultaneously self.simOp = False @@ -241,8 +251,10 @@ def __init__( name=name, ) + self._register_multiplications(self.savetaper) + @reshaped - def _matvec(self, x: NDArray) -> NDArray: + def _matvec_savetaper(self, x: NDArray) -> NDArray: ncp = get_array_module(x) if self.tapertype is not None: self.taps = to_cupy_conditional(x, self.taps) @@ -262,7 +274,7 @@ def _matvec(self, x: NDArray) -> NDArray: return y @reshaped - def _rmatvec(self, x: NDArray) -> NDArray: + def _rmatvec_savetaper(self, x: NDArray) -> NDArray: ncp = get_array_module(x) ncp_sliding_window_view = get_sliding_window_view(x) if self.tapertype is not None: @@ -281,3 +293,80 @@ def _rmatvec(self, x: NDArray) -> NDArray: self.dims[1], self.dims[2] ) return y + + @reshaped + def _matvec_nosavetaper(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + if self.tapertype is not None: + self.taps = to_cupy_conditional(x, self.taps) + y = ncp.zeros(self.dimsd, dtype=self.dtype) + if self.simOp: + x = self.Op @ x + for iwin0 in range(self.dims[0]): + if self.simOp: + xxwin = x[iwin0].reshape(self.nwin, self.dimsd[-1]) + else: + xxwin = self.Op.matvec(x[iwin0].ravel()).reshape( + self.nwin, self.dimsd[-1] + ) + if self.tapertype is not None: + if iwin0 == 0: + xxwin = self.taps[0] * xxwin + elif iwin0 == self.dims[0] - 1: + xxwin = self.taps[-1] * xxwin + else: + xxwin = self.taps[1] * xxwin + y[self.dwin_inends[0][iwin0] : self.dwin_inends[1][iwin0]] += xxwin + return y + + @reshaped + def _rmatvec_nosavetaper(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + ncp_sliding_window_view = get_sliding_window_view(x) + if self.tapertype is not None: + self.taps = to_cupy_conditional(x, self.taps) + ywins = ( + ncp_sliding_window_view(x, self.nwin, axis=0)[:: self.nwin - self.nover] + .transpose(0, 2, 1) + .copy() + ) + if self.simOp: + if self.tapertype is not None: + for iwin0 in range(self.dims[0]): + if iwin0 == 0: + ywins[0] = ywins[0] * self.taps[0] + elif iwin0 == self.dims[0] - 1: + ywins[-1] = ywins[-1] * self.taps[-1] + else: + ywins[iwin0] = ywins[iwin0] * self.taps[1] + y = self.Op.H @ ywins + else: + y = ncp.zeros(self.dims, dtype=self.dtype) + for iwin0 in range(self.dims[0]): + if iwin0 == 0: + if self.tapertype is not None: + ywins[0] = ywins[0] * self.taps[0] + y[0] = self.Op.rmatvec(ywins[0].ravel()).reshape( + self.dims[1], self.dims[2] + ) + elif iwin0 == self.dims[0] - 1: + if self.tapertype is not None: + ywins[-1] = ywins[-1] * self.taps[-1] + y[-1] = self.Op.rmatvec(ywins[-1].ravel()).reshape( + self.dims[1], self.dims[2] + ) + else: + if self.tapertype is not None: + ywins[iwin0] = ywins[iwin0] * self.taps[1] + y[iwin0] = self.Op.rmatvec(ywins[iwin0].ravel()).reshape( + self.dims[1], self.dims[2] + ) + return y + + def _register_multiplications(self, savetaper: bool) -> None: + if savetaper: + self._matvec = self._matvec_savetaper + self._rmatvec = self._rmatvec_savetaper + else: + self._matvec = self._matvec_nosavetaper + self._rmatvec = self._rmatvec_nosavetaper diff --git a/pytests/test_sliding.py b/pytests/test_sliding.py index 6750bbdf..11fb1416 100755 --- a/pytests/test_sliding.py +++ b/pytests/test_sliding.py @@ -22,6 +22,7 @@ "noverx": 0, # "winsx": 2, "tapertype": None, + "savetaper": True, } # no overlap, no taper par2 = { "ny": 6, @@ -36,6 +37,7 @@ "noverx": 0, # "winsx": 2, "tapertype": "hanning", + "savetaper": True, } # no overlap, with taper par3 = { "ny": 6, @@ -50,8 +52,24 @@ "noverx": 2, # "winsx": 4, "tapertype": None, + "savetaper": True, } # overlap, no taper par4 = { + "ny": 6, + "nx": 7, + "nt": 10, + "npy": 15, + "nwiny": 7, + "novery": 3, + # "winsy": 3, + "npx": 10, + "nwinx": 4, + "noverx": 2, + # "winsx": 4, + "tapertype": None, + "savetaper": False, +} # overlap, no taper (non saved) +par5 = { "ny": 6, "nx": 7, "nt": 10, @@ -64,10 +82,26 @@ "noverx": 2, # "winsx": 4, "tapertype": "hanning", + "savetaper": True, } # overlap, with taper +par6 = { + "ny": 6, + "nx": 7, + "nt": 10, + "npy": 15, + "nwiny": 7, + "novery": 3, + # "winsy": 3, + "npx": 10, + "nwinx": 4, + "noverx": 2, + # "winsx": 4, + "tapertype": "hanning", + "savetaper": False, +} # overlap, with taper (non saved) -@pytest.mark.parametrize("par", [(par1), (par2), (par3), (par4)]) +@pytest.mark.parametrize("par", [(par1), (par2), (par3), (par4), (par5), (par6)]) def test_Sliding1D(par): """Dot-test and inverse for Sliding1D operator""" Op = MatrixMult(np.ones((par["nwiny"], par["ny"]))) @@ -83,6 +117,7 @@ def test_Sliding1D(par): nwin=par["nwiny"], nover=par["novery"], tapertype=par["tapertype"], + savetaper=par["savetaper"], ) assert dottest(Slid, par["npy"], par["ny"] * nwins) x = np.ones(par["ny"] * nwins) @@ -92,7 +127,7 @@ def test_Sliding1D(par): assert_array_almost_equal(x.ravel(), xinv) -@pytest.mark.parametrize("par", [(par1), (par2), (par3), (par4)]) +@pytest.mark.parametrize("par", [(par1), (par2), (par3), (par4), (par5), (par6)]) def test_Sliding2D(par): """Dot-test and inverse for Sliding2D operator""" Op = MatrixMult(np.ones((par["nwiny"] * par["nt"], par["ny"] * par["nt"]))) @@ -107,6 +142,7 @@ def test_Sliding2D(par): nwin=par["nwiny"], nover=par["novery"], tapertype=par["tapertype"], + savetaper=par["savetaper"], ) assert dottest(Slid, par["npy"] * par["nt"], par["ny"] * par["nt"] * nwins) x = np.ones((par["ny"] * nwins, par["nt"])) @@ -116,7 +152,7 @@ def test_Sliding2D(par): assert_array_almost_equal(x.ravel(), xinv) -@pytest.mark.parametrize("par", [(par1), (par2), (par3), (par4)]) +@pytest.mark.parametrize("par", [(par1), (par2), (par3), (par5)]) def test_Sliding3D(par): """Dot-test and inverse for Sliding3D operator""" Op = MatrixMult( From c3da1314b298c60634b80b32412bab2f71f4e258 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Thu, 16 May 2024 21:52:07 +0300 Subject: [PATCH 06/11] feat: added savetaper to Sliding3d --- pylops/signalprocessing/sliding1d.py | 14 +-- pylops/signalprocessing/sliding2d.py | 22 ++-- pylops/signalprocessing/sliding3d.py | 178 +++++++++++++++++++++++---- pytests/test_sliding.py | 3 +- 4 files changed, 171 insertions(+), 46 deletions(-) diff --git a/pylops/signalprocessing/sliding1d.py b/pylops/signalprocessing/sliding1d.py index 8b3bb5d1..6e0fb26b 100644 --- a/pylops/signalprocessing/sliding1d.py +++ b/pylops/signalprocessing/sliding1d.py @@ -290,18 +290,14 @@ def _rmatvec_nosavetaper(self, x: NDArray) -> NDArray: else: y = ncp.zeros(self.dims, dtype=self.dtype) for iwin0 in range(self.dims[0]): - if iwin0 == 0: - if self.tapertype is not None: + if self.tapertype is not None: + if iwin0 == 0: ywins[0] = ywins[0] * self.taps[0] - y[0] = self.Op.rmatvec(ywins[0]) - elif iwin0 == self.dims[0] - 1: - if self.tapertype is not None: + elif iwin0 == self.dims[0] - 1: ywins[-1] = ywins[-1] * self.taps[-1] - y[-1] = self.Op.rmatvec(ywins[-1]) - else: - if self.tapertype is not None: + else: ywins[iwin0] = ywins[iwin0] * self.taps[1] - y[iwin0] = self.Op.rmatvec(ywins[iwin0]) + y[iwin0] = self.Op.rmatvec(ywins[iwin0]) return y def _register_multiplications(self, savetaper: bool) -> None: diff --git a/pylops/signalprocessing/sliding2d.py b/pylops/signalprocessing/sliding2d.py index 22268471..fa930662 100644 --- a/pylops/signalprocessing/sliding2d.py +++ b/pylops/signalprocessing/sliding2d.py @@ -343,24 +343,16 @@ def _rmatvec_nosavetaper(self, x: NDArray) -> NDArray: else: y = ncp.zeros(self.dims, dtype=self.dtype) for iwin0 in range(self.dims[0]): - if iwin0 == 0: - if self.tapertype is not None: + if self.tapertype is not None: + if iwin0 == 0: ywins[0] = ywins[0] * self.taps[0] - y[0] = self.Op.rmatvec(ywins[0].ravel()).reshape( - self.dims[1], self.dims[2] - ) - elif iwin0 == self.dims[0] - 1: - if self.tapertype is not None: + elif iwin0 == self.dims[0] - 1: ywins[-1] = ywins[-1] * self.taps[-1] - y[-1] = self.Op.rmatvec(ywins[-1].ravel()).reshape( - self.dims[1], self.dims[2] - ) - else: - if self.tapertype is not None: + else: ywins[iwin0] = ywins[iwin0] * self.taps[1] - y[iwin0] = self.Op.rmatvec(ywins[iwin0].ravel()).reshape( - self.dims[1], self.dims[2] - ) + y[iwin0] = self.Op.rmatvec(ywins[iwin0].ravel()).reshape( + self.dims[1], self.dims[2] + ) return y def _register_multiplications(self, savetaper: bool) -> None: diff --git a/pylops/signalprocessing/sliding3d.py b/pylops/signalprocessing/sliding3d.py index e583f171..846a678e 100644 --- a/pylops/signalprocessing/sliding3d.py +++ b/pylops/signalprocessing/sliding3d.py @@ -149,6 +149,9 @@ class Sliding3D(LinearOperator): to spatial axes in the data tapertype : :obj:`str`, optional Type of taper (``hanning``, ``cosine``, ``cosinesquare`` or ``None``) + savetaper: :obj:`bool`, optional + Save all tapers and apply them in one go (``True``) or save unique tapers and apply them one by one (``False``). + The first option is more computationally efficient, whilst the second is more memory efficient. nproc : :obj:`int`, optional *Deprecated*, will be removed in v3.0.0. Simply kept for back-compatibility with previous implementation @@ -179,6 +182,7 @@ def __init__( nover: Tuple[int, int], nop: Tuple[int, int, int], tapertype: str = "hanning", + savetaper: bool = True, nproc: int = 1, name: str = "P", ) -> None: @@ -214,56 +218,71 @@ def __init__( # create tapers self.tapertype = tapertype + self.savetaper = savetaper if self.tapertype is not None: tap = taper3d(dimsd[2], nwin, nover, tapertype=tapertype).astype(Op.dtype) - taps = [ - tap, - ] * nwins - # topmost tapers taptop = tap.copy() taptop[: nover[0]] = tap[nwin[0] // 2] - for itap in range(0, nwins1): - taps[itap] = taptop # bottommost tapers tapbottom = tap.copy() tapbottom[-nover[0] :] = tap[nwin[0] // 2] - for itap in range(nwins - nwins1, nwins): - taps[itap] = tapbottom # leftmost tapers tapleft = tap.copy() tapleft[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis] - for itap in range(0, nwins, nwins1): - taps[itap] = tapleft # rightmost tapers tapright = tap.copy() tapright[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis] - for itap in range(nwins1 - 1, nwins, nwins1): - taps[itap] = tapright # lefttopcorner taper taplefttop = tap.copy() taplefttop[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis] taplefttop[: nover[0]] = taplefttop[nwin[0] // 2] - taps[0] = taplefttop # righttopcorner taper taprighttop = tap.copy() taprighttop[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis] taprighttop[: nover[0]] = taprighttop[nwin[0] // 2] - taps[nwins1 - 1] = taprighttop # leftbottomcorner taper tapleftbottom = tap.copy() tapleftbottom[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis] tapleftbottom[-nover[0] :] = tapleftbottom[nwin[0] // 2] - taps[nwins - nwins1] = tapleftbottom # rightbottomcorner taper taprightbottom = tap.copy() taprightbottom[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis] taprightbottom[-nover[0] :] = taprightbottom[nwin[0] // 2] - taps[nwins - 1] = taprightbottom - self.taps = np.vstack(taps).reshape( - nwins0, nwins1, nwin[0], nwin[1], dimsd[2] - ) + if self.savetaper: + taps = [ + tap, + ] * nwins + + for itap in range(0, nwins1): + taps[itap] = taptop + for itap in range(nwins - nwins1, nwins): + taps[itap] = tapbottom + for itap in range(0, nwins, nwins1): + taps[itap] = tapleft + for itap in range(nwins1 - 1, nwins, nwins1): + taps[itap] = tapright + taps[0] = taplefttop + taps[nwins1 - 1] = taprighttop + taps[nwins - nwins1] = tapleftbottom + taps[nwins - 1] = taprightbottom + self.taps = np.vstack(taps).reshape( + nwins0, nwins1, nwin[0], nwin[1], dimsd[2] + ) + else: + taps = [ + taplefttop, + taptop, + taprighttop, + tapleft, + tap, + tapright, + tapleftbottom, + tapbottom, + taprightbottom, + ] + self.taps = np.vstack(taps).reshape(3, 3, nwin[0], nwin[1], dimsd[2]) # check if operator is applied to all windows simultaneously self.simOp = False if Op.shape[1] == np.prod(dims): @@ -284,8 +303,10 @@ def __init__( name=name, ) + self._register_multiplications(self.savetaper) + @reshaped - def _matvec(self, x: NDArray) -> NDArray: + def _matvec_savetaper(self, x: NDArray) -> NDArray: ncp = get_array_module(x) if self.tapertype is not None: self.taps = to_cupy_conditional(x, self.taps) @@ -313,7 +334,7 @@ def _matvec(self, x: NDArray) -> NDArray: return y @reshaped - def _rmatvec(self, x: NDArray) -> NDArray: + def _rmatvec_savetaper(self, x: NDArray) -> NDArray: ncp = get_array_module(x) ncp_sliding_window_view = get_sliding_window_view(x) if self.tapertype is not None: @@ -333,3 +354,118 @@ def _rmatvec(self, x: NDArray) -> NDArray: ywins[iwin0, iwin1].ravel() ).reshape(self.dims[2], self.dims[3], self.dims[4]) return y + + @reshaped + def _matvec_nosavetaper(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + if self.tapertype is not None: + self.taps = to_cupy_conditional(x, self.taps) + y = ncp.zeros(self.dimsd, dtype=self.dtype) + if self.simOp: + x = self.Op @ x + for iwin0 in range(self.dims[0]): + for iwin1 in range(self.dims[1]): + if self.simOp: + xxwin = x[iwin0, iwin1].reshape( + self.nwin[0], self.nwin[1], self.dimsd[-1] + ) + else: + xxwin = self.Op.matvec(x[iwin0, iwin1].ravel()).reshape( + self.nwin[0], self.nwin[1], self.dimsd[-1] + ) + if self.tapertype is not None: + if iwin0 == 0 and iwin1 == 0: + xxwin = self.taps[0, 0] * xxwin + elif iwin0 == 0 and iwin1 == self.dims[1] - 1: + xxwin = self.taps[0, -1] * xxwin + elif iwin0 == 0: + xxwin = self.taps[0, 1] * xxwin + elif iwin0 == self.dims[0] - 1 and iwin1 == 0: + xxwin = self.taps[-1, 0] * xxwin + elif iwin0 == self.dims[0] - 1 and iwin1 == self.dims[1] - 1: + xxwin = self.taps[-1, -1] * xxwin + elif iwin0 == self.dims[0] - 1: + xxwin = self.taps[-1, 1] * xxwin + elif iwin1 == 0: + xxwin = self.taps[1, 0] * xxwin + elif iwin1 == self.dims[1] - 1: + xxwin = self.taps[1, -1] * xxwin + else: + xxwin = self.taps[1, 1] * xxwin + y[ + self.dwins_inends[0][0][iwin0] : self.dwins_inends[0][1][iwin0], + self.dwins_inends[1][0][iwin1] : self.dwins_inends[1][1][iwin1], + ] += xxwin + return y + + @reshaped + def _rmatvec_nosavetaper(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + ncp_sliding_window_view = get_sliding_window_view(x) + if self.tapertype is not None: + self.taps = to_cupy_conditional(x, self.taps) + ywins = ( + ncp_sliding_window_view(x, self.nwin, axis=(0, 1))[ + :: self.nwin[0] - self.nover[0], :: self.nwin[1] - self.nover[1] + ] + .transpose(0, 1, 3, 4, 2) + .copy() + ) + if self.simOp: + if self.tapertype is not None: + for iwin0 in range(self.dims[0]): + for iwin1 in range(self.dims[1]): + if iwin0 == 0 and iwin1 == 0: + ywins[0, 0] = self.taps[0, 0] * ywins[0, 0] + elif iwin0 == 0 and iwin1 == self.dims[1] - 1: + ywins[0, -1] = self.taps[0, -1] * ywins[0, -1] + elif iwin0 == 0: + ywins[0, iwin1] = self.taps[0, 1] * ywins[0, iwin1] + elif iwin0 == self.dims[0] - 1 and iwin1 == 0: + ywins[-1, 0] = self.taps[-1, 0] * ywins[-1, 0] + elif iwin0 == self.dims[0] - 1 and iwin1 == self.dims[1] - 1: + ywins[-1, -1] = self.taps[-1, -1] * ywins[-1, -1] + elif iwin0 == self.dims[0] - 1: + ywins[-1, iwin1] = self.taps[-1, 1] * ywins[-1, iwin1] + elif iwin1 == 0: + ywins[iwin0, 0] = self.taps[1, 0] * ywins[iwin0, 0] + elif iwin1 == self.dims[1] - 1: + ywins[iwin0, -1] = self.taps[1, -1] * ywins[iwin0, -1] + else: + ywins[iwin0, iwin1] = self.taps[1, 1] * ywins[iwin0, iwin1] + y = self.Op.H @ ywins + else: + y = ncp.zeros(self.dims, dtype=self.dtype) + for iwin0 in range(self.dims[0]): + for iwin1 in range(self.dims[1]): + if self.tapertype is not None: + if iwin0 == 0 and iwin1 == 0: + ywins[0, 0] = self.taps[0, 0] * ywins[0, 0] + elif iwin0 == 0 and iwin1 == self.dims[1] - 1: + ywins[0, -1] = self.taps[0, -1] * ywins[0, -1] + elif iwin0 == 0: + ywins[0, iwin1] = self.taps[0, 1] * ywins[0, iwin1] + elif iwin0 == self.dims[0] - 1 and iwin1 == 0: + ywins[-1, 0] = self.taps[-1, 0] * ywins[-1, 0] + elif iwin0 == self.dims[0] - 1 and iwin1 == self.dims[1] - 1: + ywins[-1, -1] = self.taps[-1, -1] * ywins[-1, -1] + elif iwin0 == self.dims[0] - 1: + ywins[-1, iwin1] = self.taps[-1, 1] * ywins[-1, iwin1] + elif iwin1 == 0: + ywins[iwin0, 0] = self.taps[1, 0] * ywins[iwin0, 0] + elif iwin1 == self.dims[1] - 1: + ywins[iwin0, -1] = self.taps[1, -1] * ywins[iwin0, -1] + else: + ywins[iwin0, iwin1] = self.taps[1, 1] * ywins[iwin0, iwin1] + y[iwin0, iwin1] = self.Op.rmatvec( + ywins[iwin0, iwin1].ravel() + ).reshape(self.dims[2], self.dims[3], self.dims[4]) + return y + + def _register_multiplications(self, savetaper: bool) -> None: + if savetaper: + self._matvec = self._matvec_savetaper + self._rmatvec = self._rmatvec_savetaper + else: + self._matvec = self._matvec_nosavetaper + self._rmatvec = self._rmatvec_nosavetaper diff --git a/pytests/test_sliding.py b/pytests/test_sliding.py index 11fb1416..55a62199 100755 --- a/pytests/test_sliding.py +++ b/pytests/test_sliding.py @@ -152,7 +152,7 @@ def test_Sliding2D(par): assert_array_almost_equal(x.ravel(), xinv) -@pytest.mark.parametrize("par", [(par1), (par2), (par3), (par5)]) +@pytest.mark.parametrize("par", [(par1), (par2), (par3), (par4), (par5), (par6)]) def test_Sliding3D(par): """Dot-test and inverse for Sliding3D operator""" Op = MatrixMult( @@ -176,6 +176,7 @@ def test_Sliding3D(par): nover=(par["novery"], par["noverx"]), nop=(par["ny"], par["nx"]), tapertype=par["tapertype"], + savetaper=par["savetaper"], ) assert dottest( Slid, From 3066a4fadb9a3fd83339b8f2e933a6ffb75e0c5d Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sat, 18 May 2024 15:54:24 +0300 Subject: [PATCH 07/11] feat: added savetaper to Patch2D and Patch3D --- pylops/signalprocessing/patch2d.py | 166 +++++++- pylops/signalprocessing/patch3d.py | 646 ++++++++++++++++++++++++----- pytests/test_patching.py | 49 ++- 3 files changed, 741 insertions(+), 120 deletions(-) diff --git a/pylops/signalprocessing/patch2d.py b/pylops/signalprocessing/patch2d.py index 1d54e97a..e60fd807 100644 --- a/pylops/signalprocessing/patch2d.py +++ b/pylops/signalprocessing/patch2d.py @@ -141,6 +141,9 @@ class Patch2D(LinearOperator): Size of model in the transformed domain tapertype : :obj:`str`, optional Type of taper (``hanning``, ``cosine``, ``cosinesquare`` or ``None``) + savetaper: :obj:`bool`, optional + Save all tapers and apply them in one go (``True``) or save unique tapers and apply them one by one (``False``). + The first option is more computationally efficient, whilst the second is more memory efficient. scalings : :obj:`tuple` or :obj:`list`, optional Set of scalings to apply to each patch. If ``None``, no scale will be applied @@ -178,6 +181,7 @@ def __init__( nover: Tuple[int, int], nop: Tuple[int, int], tapertype: str = "hanning", + savetaper: bool = True, scalings: Optional[Sequence[float]] = None, name: str = "P", ) -> None: @@ -206,52 +210,68 @@ def __init__( # create tapers self.tapertype = tapertype + self.savetaper = savetaper if self.tapertype is not None: tap = taper2d(nwin[1], nwin[0], nover, tapertype=tapertype).astype(Op.dtype) - taps = [ - tap, - ] * nwins # topmost tapers taptop = tap.copy() taptop[: nover[0]] = tap[nwin[0] // 2] - for itap in range(0, nwins1): - taps[itap] = taptop # bottommost tapers tapbottom = tap.copy() tapbottom[-nover[0] :] = tap[nwin[0] // 2] - for itap in range(nwins - nwins1, nwins): - taps[itap] = tapbottom # leftmost tapers tapleft = tap.copy() tapleft[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis] - for itap in range(0, nwins, nwins1): - taps[itap] = tapleft # rightmost tapers tapright = tap.copy() tapright[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis] - for itap in range(nwins1 - 1, nwins, nwins1): - taps[itap] = tapright # lefttopcorner taper taplefttop = tap.copy() taplefttop[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis] taplefttop[: nover[0]] = taplefttop[nwin[0] // 2] - taps[0] = taplefttop # righttopcorner taper taprighttop = tap.copy() taprighttop[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis] taprighttop[: nover[0]] = taprighttop[nwin[0] // 2] - taps[nwins1 - 1] = taprighttop # leftbottomcorner taper tapleftbottom = tap.copy() tapleftbottom[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis] tapleftbottom[-nover[0] :] = tapleftbottom[nwin[0] // 2] - taps[nwins - nwins1] = tapleftbottom # rightbottomcorner taper taprightbottom = tap.copy() taprightbottom[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis] taprightbottom[-nover[0] :] = taprightbottom[nwin[0] // 2] - taps[nwins - 1] = taprightbottom - self.taps = np.vstack(taps).reshape(nwins0, nwins1, nwin[0], nwin[1]) + + if self.savetaper: + taps = [ + tap, + ] * nwins + for itap in range(0, nwins1): + taps[itap] = taptop + for itap in range(nwins - nwins1, nwins): + taps[itap] = tapbottom + for itap in range(0, nwins, nwins1): + taps[itap] = tapleft + for itap in range(nwins1 - 1, nwins, nwins1): + taps[itap] = tapright + taps[0] = taplefttop + taps[nwins1 - 1] = taprighttop + taps[nwins - nwins1] = tapleftbottom + taps[nwins - 1] = taprightbottom + self.taps = np.vstack(taps).reshape(nwins0, nwins1, nwin[0], nwin[1]) + else: + taps = [ + taplefttop, + taptop, + taprighttop, + tapleft, + tap, + tapright, + tapleftbottom, + tapbottom, + taprightbottom, + ] + self.taps = np.vstack(taps).reshape(3, 3, nwin[0], nwin[1]) # define scalings if scalings is None: @@ -273,8 +293,10 @@ def __init__( name=name, ) + self._register_multiplications(self.savetaper) + @reshaped() - def _matvec(self, x: NDArray) -> NDArray: + def _matvec_savetaper(self, x: NDArray) -> NDArray: ncp = get_array_module(x) if self.tapertype is not None: self.taps = to_cupy_conditional(x, self.taps) @@ -299,7 +321,7 @@ def _matvec(self, x: NDArray) -> NDArray: return y @reshaped - def _rmatvec(self, x: NDArray) -> NDArray: + def _rmatvec_savetaper(self, x: NDArray) -> NDArray: ncp = get_array_module(x) ncp_sliding_window_view = get_sliding_window_view(x) if self.tapertype is not None: @@ -319,3 +341,111 @@ def _rmatvec(self, x: NDArray) -> NDArray: ywins[iwin0, iwin1].ravel() ).reshape(self.dims[2], self.dims[3]) return y + + @reshaped() + def _matvec_nosavetaper(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + if self.tapertype is not None: + self.taps = to_cupy_conditional(x, self.taps) + y = ncp.zeros(self.dimsd, dtype=self.dtype) + if self.simOp: + x = self.Op @ x + for iwin0 in range(self.dims[0]): + for iwin1 in range(self.dims[1]): + if self.simOp: + xxwin = x[iwin0, iwin1].reshape(self.nwin) + else: + xxwin = self.Op.matvec(x[iwin0, iwin1].ravel()).reshape(self.nwin) + if self.tapertype is not None: + if iwin0 == 0 and iwin1 == 0: + xxwin = self.taps[0, 0] * xxwin + elif iwin0 == 0 and iwin1 == self.dims[1] - 1: + xxwin = self.taps[0, -1] * xxwin + elif iwin0 == 0: + xxwin = self.taps[0, 1] * xxwin + elif iwin0 == self.dims[0] - 1 and iwin1 == 0: + xxwin = self.taps[-1, 0] * xxwin + elif iwin0 == self.dims[0] - 1 and iwin1 == self.dims[1] - 1: + xxwin = self.taps[-1, -1] * xxwin + elif iwin0 == self.dims[0] - 1: + xxwin = self.taps[-1, 1] * xxwin + elif iwin1 == 0: + xxwin = self.taps[1, 0] * xxwin + elif iwin1 == self.dims[1] - 1: + xxwin = self.taps[1, -1] * xxwin + else: + xxwin = self.taps[1, 1] * xxwin + + y[ + self.dwins_inends[0][0][iwin0] : self.dwins_inends[0][1][iwin0], + self.dwins_inends[1][0][iwin1] : self.dwins_inends[1][1][iwin1], + ] += xxwin + return y + + @reshaped + def _rmatvec_nosavetaper(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + ncp_sliding_window_view = get_sliding_window_view(x) + if self.tapertype is not None: + self.taps = to_cupy_conditional(x, self.taps) + ywins = ncp_sliding_window_view(x, self.nwin)[ + :: self.nwin[0] - self.nover[0], :: self.nwin[1] - self.nover[1] + ].copy() + if self.simOp: + if self.tapertype is not None: + for iwin0 in range(self.dims[0]): + for iwin1 in range(self.dims[1]): + if iwin0 == 0 and iwin1 == 0: + ywins[0, 0] = self.taps[0, 0] * ywins[0, 0] + elif iwin0 == 0 and iwin1 == self.dims[1] - 1: + ywins[0, -1] = self.taps[0, -1] * ywins[0, -1] + elif iwin0 == 0: + ywins[0, iwin1] = self.taps[0, 1] * ywins[0, iwin1] + elif iwin0 == self.dims[0] - 1 and iwin1 == 0: + ywins[-1, 0] = self.taps[-1, 0] * ywins[-1, 0] + elif iwin0 == self.dims[0] - 1 and iwin1 == self.dims[1] - 1: + ywins[-1, -1] = self.taps[-1, -1] * ywins[-1, -1] + elif iwin0 == self.dims[0] - 1: + ywins[-1, iwin1] = self.taps[-1, 1] * ywins[-1, iwin1] + elif iwin1 == 0: + ywins[iwin0, 0] = self.taps[1, 0] * ywins[iwin0, 0] + elif iwin1 == self.dims[1] - 1: + ywins[iwin0, -1] = self.taps[1, -1] * ywins[iwin0, -1] + else: + ywins[iwin0, iwin1] = self.taps[1, 1] * ywins[iwin0, iwin1] + y = self.Op.H @ ywins + else: + y = ncp.zeros(self.dims, dtype=self.dtype) + for iwin0 in range(self.dims[0]): + for iwin1 in range(self.dims[1]): + if self.tapertype is not None: + if iwin0 == 0 and iwin1 == 0: + ywins[0, 0] = self.taps[0, 0] * ywins[0, 0] + elif iwin0 == 0 and iwin1 == self.dims[1] - 1: + ywins[0, -1] = self.taps[0, -1] * ywins[0, -1] + elif iwin0 == 0: + ywins[0, iwin1] = self.taps[0, 1] * ywins[0, iwin1] + elif iwin0 == self.dims[0] - 1 and iwin1 == 0: + ywins[-1, 0] = self.taps[-1, 0] * ywins[-1, 0] + elif iwin0 == self.dims[0] - 1 and iwin1 == self.dims[1] - 1: + ywins[-1, -1] = self.taps[-1, -1] * ywins[-1, -1] + elif iwin0 == self.dims[0] - 1: + ywins[-1, iwin1] = self.taps[-1, 1] * ywins[-1, iwin1] + elif iwin1 == 0: + ywins[iwin0, 0] = self.taps[1, 0] * ywins[iwin0, 0] + elif iwin1 == self.dims[1] - 1: + ywins[iwin0, -1] = self.taps[1, -1] * ywins[iwin0, -1] + else: + ywins[iwin0, iwin1] = self.taps[1, 1] * ywins[iwin0, iwin1] + y[iwin0, iwin1] = self.Op.rmatvec( + ywins[iwin0, iwin1].ravel() + ).reshape(self.dims[2], self.dims[3]) + return y + + def _register_multiplications(self, savetaper: bool) -> None: + if savetaper: + self._matvec = self._matvec_savetaper + self._rmatvec = self._rmatvec_savetaper + else: + self._matvec = self._matvec_nosavetaper + self._rmatvec = self._rmatvec_nosavetaper diff --git a/pylops/signalprocessing/patch3d.py b/pylops/signalprocessing/patch3d.py index d7a87747..7a893b54 100644 --- a/pylops/signalprocessing/patch3d.py +++ b/pylops/signalprocessing/patch3d.py @@ -156,6 +156,9 @@ class Patch3D(LinearOperator): Size of model in the transformed domain tapertype : :obj:`str`, optional Type of taper (``hanning``, ``cosine``, ``cosinesquare`` or ``None``) + savetaper: :obj:`bool`, optional + Save all tapers and apply them in one go (``True``) or save unique tapers and apply them one by one (``False``). + The first option is more computationally efficient, whilst the second is more memory efficient. scalings : :obj:`tuple` or :obj:`list`, optional Set of scalings to apply to each patch. If ``None``, no scale will be applied @@ -191,6 +194,7 @@ def __init__( nover: Tuple[int, int, int], nop: Tuple[int, int, int], tapertype: str = "hanning", + savetaper: bool = True, scalings: Optional[Sequence[float]] = None, name: str = "P", ) -> None: @@ -229,207 +233,254 @@ def __init__( # create tapers self.tapertype = tapertype + self.savetaper = savetaper if tapertype is not None: tap = tapernd(nwin, nover, tapertype=tapertype).astype(Op.dtype) - taps = [ - tap, - ] * nwins # 1, sides # topmost tapers taptop = tap.copy() taptop[: nover[0]] = tap[nwin[0] // 2] - for itap in range(0, nwins1 * nwins2): - taps[itap] = taptop # bottommost tapers tapbottom = tap.copy() tapbottom[-nover[0] :] = tap[nwin[0] // 2] - for itap in range(nwins - nwins1 * nwins2, nwins): - taps[itap] = tapbottom # frontmost tapers tapfront = tap.copy() tapfront[:, :, : nover[2]] = tap[:, :, nwin[2] // 2][:, :, np.newaxis] - for itap in range(0, nwins, nwins2): - taps[itap] = tapfront # backmost tapers tapback = tap.copy() tapback[:, :, -nover[2] :] = tap[:, :, nwin[2] // 2][:, :, np.newaxis] - for itap in range(nwins2 - 1, nwins, nwins2): - taps[itap] = tapback # leftmost tapers tapleft = tap.copy() tapleft[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis, :] - for itap in range(0, nwins, nwins1 * nwins2): - for i in range(nwins2): - taps[itap + i] = tapleft # rightmost tapers tapright = tap.copy() tapright[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis, :] - for itap in range(nwins2 * (nwins1 - 1), nwins, nwins2 * nwins1): - for i in range(nwins2): - taps[itap + i] = tapright + # 2. pillars # topleftmost tapers taplefttop = tap.copy() taplefttop[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis, :] taplefttop[: nover[0]] = taplefttop[nwin[0] // 2] - for itap in range(nwins2): - taps[itap] = taplefttop # toprightmost tapers taprighttop = tap.copy() taprighttop[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis, :] taprighttop[: nover[0]] = taprighttop[nwin[0] // 2] - for itap in range(nwins2): - taps[nwins2 * (nwins1 - 1) + itap] = taprighttop # topfrontmost tapers tapfronttop = tap.copy() tapfronttop[:, :, : nover[2]] = tap[:, :, nwin[2] // 2][:, :, np.newaxis] tapfronttop[: nover[0]] = tapfronttop[nwin[0] // 2] - for itap in range(0, nwins1 * nwins2, nwins2): - taps[itap] = tapfronttop # topbackmost tapers tapbacktop = tap.copy() tapbacktop[:, :, -nover[2] :] = tap[:, :, nwin[2] // 2][:, :, np.newaxis] tapbacktop[: nover[0]] = tapbacktop[nwin[0] // 2] - for itap in range(nwins2 - 1, nwins1 * nwins2, nwins2): - taps[itap] = tapbacktop # bottomleftmost tapers tapleftbottom = tap.copy() tapleftbottom[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis, :] tapleftbottom[-nover[0] :] = tapleftbottom[nwin[0] // 2] - for itap in range(nwins2): - taps[(nwins0 - 1) * nwins1 * nwins2 + itap] = tapleftbottom # bottomrightmost tapers taprightbottom = tap.copy() taprightbottom[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis, :] taprightbottom[-nover[0] :] = taprightbottom[nwin[0] // 2] - for itap in range(nwins2): - taps[ - (nwins0 - 1) * nwins1 * nwins2 + (nwins1 - 1) * nwins2 + itap - ] = taprightbottom # bottomfrontmost tapers tapfrontbottom = tap.copy() tapfrontbottom[:, :, : nover[2]] = tap[:, :, nwin[2] // 2][:, :, np.newaxis] tapfrontbottom[-nover[0] :] = tapfrontbottom[nwin[0] // 2] - for itap in range(0, nwins1 * nwins2, nwins2): - taps[(nwins0 - 1) * nwins1 * nwins2 + itap] = tapfrontbottom # bottombackmost tapers tapbackbottom = tap.copy() tapbackbottom[:, :, -nover[2] :] = tap[:, :, nwin[2] // 2][:, :, np.newaxis] tapbackbottom[-nover[0] :] = tapbackbottom[nwin[0] // 2] - for itap in range(0, nwins1 * nwins2, nwins2): - taps[(nwins0 - 1) * nwins1 * nwins2 + nwins2 + itap - 1] = tapbackbottom # leftfrontmost tapers tapleftfront = tap.copy() tapleftfront[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis, :] tapleftfront[:, :, : nover[2]] = tapleftfront[:, :, nwin[2] // 2][ :, :, np.newaxis ] - for itap in range(0, nwins, nwins1 * nwins2): - taps[itap] = tapleftfront # rightfrontmost tapers taprightfront = tap.copy() taprightfront[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis, :] taprightfront[:, :, : nover[2]] = taprightfront[:, :, nwin[2] // 2][ :, :, np.newaxis ] - for itap in range(0, nwins, nwins1 * nwins2): - taps[(nwins1 - 1) * nwins2 + itap] = taprightfront # leftbackmost tapers tapleftback = tap.copy() tapleftback[:, : nover[1]] = tap[:, nwin[1] // 2][:, np.newaxis, :] tapleftback[:, :, -nover[2] :] = tapleftback[:, :, nwin[2] // 2][ :, :, np.newaxis ] - for itap in range(0, nwins, nwins1 * nwins2): - taps[nwins2 + itap - 1] = tapleftback # rightbackmost tapers taprightback = tap.copy() taprightback[:, -nover[1] :] = tap[:, nwin[1] // 2][:, np.newaxis, :] taprightback[:, :, -nover[2] :] = taprightback[:, :, nwin[2] // 2][ :, :, np.newaxis ] - for itap in range(0, nwins, nwins1 * nwins2): - taps[(nwins1 - 1) * nwins2 + nwins2 + itap - 1] = taprightback + # 3. corners # lefttopfrontcorner taper - taplefttop = tap.copy() - taplefttop[: nover[0]] = tap[nwin[0] // 2] - taplefttop[:, : nover[1]] = taplefttop[:, nwin[1] // 2][:, np.newaxis, :] - taplefttop[:, :, : nover[2]] = taplefttop[:, :, nwin[2] // 2][ + taplefttopfront = tap.copy() + taplefttopfront[: nover[0]] = tap[nwin[0] // 2] + taplefttopfront[:, : nover[1]] = taplefttopfront[:, nwin[1] // 2][ + :, np.newaxis, : + ] + taplefttopfront[:, :, : nover[2]] = taplefttopfront[:, :, nwin[2] // 2][ :, :, np.newaxis ] - taps[0] = taplefttop # lefttopbackcorner taper - taplefttop = tap.copy() - taplefttop[: nover[0]] = tap[nwin[0] // 2] - taplefttop[:, : nover[1]] = taplefttop[:, nwin[1] // 2][:, np.newaxis, :] - taplefttop[:, :, -nover[2] :] = taplefttop[:, :, nwin[2] // 2][ + taplefttopback = tap.copy() + taplefttopback[: nover[0]] = tap[nwin[0] // 2] + taplefttopback[:, : nover[1]] = taplefttopback[:, nwin[1] // 2][ + :, np.newaxis, : + ] + taplefttopback[:, :, -nover[2] :] = taplefttopback[:, :, nwin[2] // 2][ :, :, np.newaxis ] - taps[nwins2 - 1] = taplefttop # righttopfrontcorner taper - taprighttop = tap.copy() - taprighttop[: nover[0]] = tap[nwin[0] // 2] - taprighttop[:, -nover[1] :] = taprighttop[:, nwin[1] // 2][:, np.newaxis, :] - taprighttop[:, :, : nover[2]] = taprighttop[:, :, nwin[2] // 2][ + taprighttopfront = tap.copy() + taprighttopfront[: nover[0]] = tap[nwin[0] // 2] + taprighttopfront[:, -nover[1] :] = taprighttopfront[:, nwin[1] // 2][ + :, np.newaxis, : + ] + taprighttopfront[:, :, : nover[2]] = taprighttopfront[:, :, nwin[2] // 2][ :, :, np.newaxis ] - taps[(nwins1 - 1) * nwins2] = taprighttop # righttopbackcorner taper - taprighttop = tap.copy() - taprighttop[: nover[0]] = tap[nwin[0] // 2] - taprighttop[:, -nover[1] :] = taprighttop[:, nwin[1] // 2][:, np.newaxis, :] - taprighttop[:, :, -nover[2] :] = taprighttop[:, :, nwin[2] // 2][ + taprighttopback = tap.copy() + taprighttopback[: nover[0]] = tap[nwin[0] // 2] + taprighttopback[:, -nover[1] :] = taprighttopback[:, nwin[1] // 2][ + :, np.newaxis, : + ] + taprighttopback[:, :, -nover[2] :] = taprighttopback[:, :, nwin[2] // 2][ :, :, np.newaxis ] - taps[(nwins1 - 1) * nwins2 + nwins2 - 1] = taprighttop # leftbottomfrontcorner taper - tapleftbottom = tap.copy() - tapleftbottom[-nover[0] :] = tap[nwin[0] // 2] - tapleftbottom[:, : nover[1]] = tapleftbottom[:, nwin[1] // 2][ + tapleftbottomfront = tap.copy() + tapleftbottomfront[-nover[0] :] = tap[nwin[0] // 2] + tapleftbottomfront[:, : nover[1]] = tapleftbottomfront[:, nwin[1] // 2][ :, np.newaxis, : ] - tapleftbottom[:, :, : nover[2]] = tapleftbottom[:, :, nwin[2] // 2][ - :, :, np.newaxis - ] - taps[(nwins0 - 1) * nwins1 * nwins2] = tapleftbottom + tapleftbottomfront[:, :, : nover[2]] = tapleftbottomfront[ + :, :, nwin[2] // 2 + ][:, :, np.newaxis] # leftbottombackcorner taper - tapleftbottom = tap.copy() - tapleftbottom[-nover[0] :] = tap[nwin[0] // 2] - tapleftbottom[:, : nover[1]] = tapleftbottom[:, nwin[1] // 2][ + tapleftbottomback = tap.copy() + tapleftbottomback[-nover[0] :] = tap[nwin[0] // 2] + tapleftbottomback[:, : nover[1]] = tapleftbottomback[:, nwin[1] // 2][ :, np.newaxis, : ] - tapleftbottom[:, :, -nover[2] :] = tapleftbottom[:, :, nwin[2] // 2][ - :, :, np.newaxis - ] - taps[(nwins0 - 1) * nwins1 * nwins2 + nwins2 - 1] = tapleftbottom + tapleftbottomback[:, :, -nover[2] :] = tapleftbottomback[ + :, :, nwin[2] // 2 + ][:, :, np.newaxis] # rightbottomfrontcorner taper - taprightbottom = tap.copy() - taprightbottom[-nover[0] :] = tap[nwin[0] // 2] - taprightbottom[:, -nover[1] :] = taprightbottom[:, nwin[1] // 2][ + taprightbottomfront = tap.copy() + taprightbottomfront[-nover[0] :] = tap[nwin[0] // 2] + taprightbottomfront[:, -nover[1] :] = taprightbottomfront[:, nwin[1] // 2][ :, np.newaxis, : ] - taprightbottom[:, :, : nover[2]] = taprightbottom[:, :, nwin[2] // 2][ - :, :, np.newaxis - ] - taps[ - (nwins0 - 1) * nwins1 * nwins2 + (nwins1 - 1) * nwins2 - ] = taprightbottom + taprightbottomfront[:, :, : nover[2]] = taprightbottomfront[ + :, :, nwin[2] // 2 + ][:, :, np.newaxis] # rightbottombackcorner taper - taprightbottom = tap.copy() - taprightbottom[-nover[0] :] = tap[nwin[0] // 2] - taprightbottom[:, -nover[1] :] = taprightbottom[:, nwin[1] // 2][ + taprightbottomback = tap.copy() + taprightbottomback[-nover[0] :] = tap[nwin[0] // 2] + taprightbottomback[:, -nover[1] :] = taprightbottomback[:, nwin[1] // 2][ :, np.newaxis, : ] - taprightbottom[:, :, -nover[2] :] = taprightbottom[:, :, nwin[2] // 2][ - :, :, np.newaxis - ] - taps[ - (nwins0 - 1) * nwins1 * nwins2 + (nwins1 - 1) * nwins2 + nwins2 - 1 - ] = taprightbottom - self.taps = np.vstack(taps).reshape( - nwins0, nwins1, nwins2, nwin[0], nwin[1], nwin[2] - ) + taprightbottomback[:, :, -nover[2] :] = taprightbottomback[ + :, :, nwin[2] // 2 + ][:, :, np.newaxis] + if self.savetaper: + taps = [ + tap, + ] * nwins + for itap in range(0, nwins1 * nwins2): + taps[itap] = taptop + for itap in range(nwins - nwins1 * nwins2, nwins): + taps[itap] = tapbottom + for itap in range(0, nwins, nwins2): + taps[itap] = tapfront + for itap in range(nwins2 - 1, nwins, nwins2): + taps[itap] = tapback + for itap in range(0, nwins, nwins1 * nwins2): + for i in range(nwins2): + taps[itap + i] = tapleft + for itap in range(nwins2 * (nwins1 - 1), nwins, nwins2 * nwins1): + for i in range(nwins2): + taps[itap + i] = tapright + for itap in range(nwins2): + taps[itap] = taplefttop + for itap in range(nwins2): + taps[nwins2 * (nwins1 - 1) + itap] = taprighttop + for itap in range(0, nwins1 * nwins2, nwins2): + taps[itap] = tapfronttop + for itap in range(nwins2 - 1, nwins1 * nwins2, nwins2): + taps[itap] = tapbacktop + for itap in range(nwins2): + taps[(nwins0 - 1) * nwins1 * nwins2 + itap] = tapleftbottom + for itap in range(nwins2): + taps[ + (nwins0 - 1) * nwins1 * nwins2 + (nwins1 - 1) * nwins2 + itap + ] = taprightbottom + for itap in range(0, nwins1 * nwins2, nwins2): + taps[(nwins0 - 1) * nwins1 * nwins2 + itap] = tapfrontbottom + for itap in range(0, nwins1 * nwins2, nwins2): + taps[ + (nwins0 - 1) * nwins1 * nwins2 + nwins2 + itap - 1 + ] = tapbackbottom + for itap in range(0, nwins, nwins1 * nwins2): + taps[itap] = tapleftfront + for itap in range(0, nwins, nwins1 * nwins2): + taps[(nwins1 - 1) * nwins2 + itap] = taprightfront + for itap in range(0, nwins, nwins1 * nwins2): + taps[nwins2 + itap - 1] = tapleftback + for itap in range(0, nwins, nwins1 * nwins2): + taps[(nwins1 - 1) * nwins2 + nwins2 + itap - 1] = taprightback + taps[0] = taplefttopfront + taps[nwins2 - 1] = taplefttopback + taps[(nwins1 - 1) * nwins2] = taprighttopfront + taps[(nwins1 - 1) * nwins2 + nwins2 - 1] = taprighttopback + taps[(nwins0 - 1) * nwins1 * nwins2] = tapleftbottomfront + taps[(nwins0 - 1) * nwins1 * nwins2 + nwins2 - 1] = tapleftbottomback + taps[ + (nwins0 - 1) * nwins1 * nwins2 + (nwins1 - 1) * nwins2 + ] = taprightbottomfront + taps[ + (nwins0 - 1) * nwins1 * nwins2 + (nwins1 - 1) * nwins2 + nwins2 - 1 + ] = taprightbottomback + self.taps = np.vstack(taps).reshape( + nwins0, nwins1, nwins2, nwin[0], nwin[1], nwin[2] + ) + else: + self.taps = np.zeros( + (3, 3, 3, nwin[0], nwin[1], nwin[2]), dtype=Op.dtype + ) + self.taps[0, 0, 0] = taplefttopfront + self.taps[0, 0, 1] = taplefttop + self.taps[0, 0, -1] = taplefttopback + self.taps[0, 1, 0] = tapfronttop + self.taps[0, 1, 1] = taptop + self.taps[0, 1, -1] = tapbacktop + self.taps[0, -1, 0] = taprighttopfront + self.taps[0, -1, 1] = taprighttop + self.taps[0, -1, -1] = taprighttopback + + self.taps[1, 0, 0] = tapleftfront + self.taps[1, 0, 1] = tapleft + self.taps[1, 0, -1] = tapleftback + self.taps[1, 1, 0] = tapfront + self.taps[1, 1, 1] = tap + self.taps[1, 1, -1] = tapback + self.taps[1, -1, 0] = taprightfront + self.taps[1, -1, 1] = tapright + self.taps[1, -1, -1] = taprightback + + self.taps[-1, 0, 0] = tapleftbottomfront + self.taps[-1, 0, 1] = tapleftbottom + self.taps[-1, 0, -1] = tapleftbottomback + self.taps[-1, 1, 0] = tapfrontbottom + self.taps[-1, 1, 1] = tapbottom + self.taps[-1, 1, -1] = tapbackbottom + self.taps[-1, -1, 0] = taprightbottomfront + self.taps[-1, -1, 1] = taprightbottom + self.taps[-1, -1, -1] = taprightbottomback # define scalings if scalings is None: @@ -458,8 +509,10 @@ def __init__( name=name, ) + self._register_multiplications(self.savetaper) + @reshaped() - def _matvec(self, x: NDArray) -> NDArray: + def _matvec_savetaper(self, x: NDArray) -> NDArray: ncp = get_array_module(x) if self.tapertype is not None: self.taps = to_cupy_conditional(x, self.taps) @@ -488,7 +541,7 @@ def _matvec(self, x: NDArray) -> NDArray: return y @reshaped - def _rmatvec(self, x: NDArray) -> NDArray: + def _rmatvec_savetaper(self, x: NDArray) -> NDArray: ncp = get_array_module(x) ncp_sliding_window_view = get_sliding_window_view(x) if self.tapertype is not None: @@ -511,3 +564,396 @@ def _rmatvec(self, x: NDArray) -> NDArray: ywins[iwin0, iwin1, iwin2].ravel() ).reshape(self.dims[3], self.dims[4], self.dims[5]) return y + + @reshaped() + def _matvec_nosavetaper(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + if self.tapertype is not None: + self.taps = to_cupy_conditional(x, self.taps) + y = ncp.zeros(self.dimsd, dtype=self.dtype) + if self.simOp: + x = self.Op @ x + for iwin0 in range(self.dims[0]): + for iwin1 in range(self.dims[1]): + for iwin2 in range(self.dims[2]): + if self.simOp: + xxwin = x[iwin0, iwin1, iwin2].reshape(self.nwin) + else: + xxwin = self.Op.matvec(x[iwin0, iwin1, iwin2].ravel()).reshape( + self.nwin + ) + if self.tapertype is not None: + if iwin0 == 0 and iwin1 == 0 and iwin2 == 0: + xxwin = self.taps[0, 0, 0] * xxwin + elif iwin0 == 0 and iwin1 == 0 and iwin2 == self.dims[2] - 1: + xxwin = self.taps[0, 0, -1] * xxwin + elif ( + iwin0 == 0 + and iwin1 == self.dims[1] - 1 + and iwin2 == self.dims[2] - 1 + ): + xxwin = self.taps[0, -1, -1] * xxwin + elif iwin0 == 0 and iwin1 == self.dims[1] - 1 and iwin2 == 0: + xxwin = self.taps[0, -1, 0] * xxwin + elif iwin0 == 0 and iwin1 == 0: + xxwin = self.taps[0, 0, 1] * xxwin + elif iwin0 == 0 and iwin1 == self.dims[1] - 1: + xxwin = self.taps[0, -1, 1] * xxwin + elif iwin0 == 0 and iwin2 == 0: + xxwin = self.taps[0, 1, 0] * xxwin + elif iwin0 == 0 and iwin2 == self.dims[2] - 1: + xxwin = self.taps[0, 1, -1] * xxwin + elif iwin0 == 0: + xxwin = self.taps[0, 1, 1] * xxwin + + elif iwin0 == self.dims[0] - 1 and iwin1 == 0 and iwin2 == 0: + xxwin = self.taps[-1, 0, 0] * xxwin + elif ( + iwin0 == self.dims[0] - 1 + and iwin1 == 0 + and iwin2 == self.dims[2] - 1 + ): + xxwin = self.taps[-1, 0, -1] * xxwin + elif ( + iwin0 == self.dims[0] - 1 + and iwin1 == self.dims[1] - 1 + and iwin2 == self.dims[2] - 1 + ): + xxwin = self.taps[-1, -1, -1] * xxwin + elif ( + iwin0 == self.dims[0] - 1 + and iwin1 == self.dims[1] - 1 + and iwin2 == 0 + ): + xxwin = self.taps[-1, -1, 0] * xxwin + elif iwin0 == self.dims[0] - 1 and iwin1 == 0: + xxwin = self.taps[-1, 0, 1] * xxwin + elif iwin0 == self.dims[0] - 1 and iwin1 == self.dims[1] - 1: + xxwin = self.taps[-1, -1, 1] * xxwin + elif iwin0 == self.dims[0] - 1 and iwin2 == 0: + xxwin = self.taps[-1, 1, 0] * xxwin + elif iwin0 == self.dims[0] - 1 and iwin2 == self.dims[2] - 1: + xxwin = self.taps[-1, 1, -1] * xxwin + elif iwin0 == self.dims[0] - 1: + xxwin = self.taps[-1, 1, 1] * xxwin + + elif iwin1 == 0 and iwin2 == 0: + xxwin = self.taps[1, 0, 0] * xxwin + elif iwin1 == 0 and iwin2 == self.dims[2] - 1: + xxwin = self.taps[1, 0, -1] * xxwin + elif iwin1 == self.dims[1] - 1 and iwin2 == self.dims[2] - 1: + xxwin = self.taps[1, -1, -1] * xxwin + elif iwin1 == self.dims[1] - 1 and iwin2 == 0: + xxwin = self.taps[1, -1, 0] * xxwin + elif iwin1 == 0: + xxwin = self.taps[1, 0, 1] * xxwin + elif iwin1 == self.dims[1] - 1: + xxwin = self.taps[1, -1, 1] * xxwin + elif iwin2 == 0: + xxwin = self.taps[1, 1, 0] * xxwin + elif iwin2 == self.dims[2] - 1: + xxwin = self.taps[1, 1, -1] * xxwin + else: + xxwin = self.taps[1, 1, 1] * xxwin + y[ + self.dwins_inends[0][0][iwin0] : self.dwins_inends[0][1][iwin0], + self.dwins_inends[1][0][iwin1] : self.dwins_inends[1][1][iwin1], + self.dwins_inends[2][0][iwin2] : self.dwins_inends[2][1][iwin2], + ] += xxwin + return y + + @reshaped() + def _rmatvec_nosavetaper(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + ncp_sliding_window_view = get_sliding_window_view(x) + if self.tapertype is not None: + self.taps = to_cupy_conditional(x, self.taps) + ywins = ncp_sliding_window_view(x, self.nwin)[ + :: self.nwin[0] - self.nover[0], + :: self.nwin[1] - self.nover[1], + :: self.nwin[2] - self.nover[2], + ].copy() + if self.simOp: + if self.tapertype is not None: + for iwin0 in range(self.dims[0]): + for iwin1 in range(self.dims[1]): + for iwin2 in range(self.dims[2]): + if iwin0 == 0 and iwin1 == 0 and iwin2 == 0: + ywins[0, 0, 0] = self.taps[0, 0, 0] * ywins[0, 0, 0] + elif ( + iwin0 == 0 and iwin1 == 0 and iwin2 == self.dims[2] - 1 + ): + ywins[0, 0, -1] = self.taps[0, 0, -1] * ywins[0, 0, -1] + elif ( + iwin0 == 0 + and iwin1 == self.dims[1] - 1 + and iwin2 == self.dims[2] - 1 + ): + ywins[0, -1, -1] = ( + self.taps[0, -1, -1] * ywins[0, -1, -1] + ) + elif ( + iwin0 == 0 and iwin1 == self.dims[1] - 1 and iwin2 == 0 + ): + ywins[0, -1, 0] = self.taps[0, -1, 0] * ywins[0, -1, 0] + elif iwin0 == 0 and iwin1 == 0: + ywins[0, 0, iwin2] = ( + self.taps[0, 0, 1] * ywins[0, 0, iwin2] + ) + elif iwin0 == 0 and iwin1 == self.dims[1] - 1: + ywins[0, -1, iwin2] = ( + self.taps[0, -1, 1] * ywins[0, -1, iwin2] + ) + elif iwin0 == 0 and iwin2 == 0: + ywins[0, iwin1, 0] = ( + self.taps[0, 1, 0] * ywins[0, iwin1, 0] + ) + elif iwin0 == 0 and iwin2 == self.dims[2] - 1: + ywins[0, iwin1, -1] = ( + self.taps[0, 1, -1] * ywins[0, iwin1, -1] + ) + elif iwin0 == 0: + ywins[0, iwin1, iwin2] = ( + self.taps[0, 1, 1] * ywins[0, iwin1, iwin2] + ) + + elif ( + iwin0 == self.dims[0] - 1 and iwin1 == 0 and iwin2 == 0 + ): + ywins[-1, 0, 0] = self.taps[-1, 0, 0] * ywins[-1, 0, 0] + elif ( + iwin0 == self.dims[0] - 1 + and iwin1 == 0 + and iwin2 == self.dims[2] - 1 + ): + ywins[-1, 0, -1] = ( + self.taps[-1, 0, -1] * ywins[-1, 0, -1] + ) + elif ( + iwin0 == self.dims[0] - 1 + and iwin1 == self.dims[1] - 1 + and iwin2 == self.dims[2] - 1 + ): + ywins[-1, -1, -1] = ( + self.taps[-1, -1, -1] * ywins[-1, -1, -1] + ) + elif ( + iwin0 == self.dims[0] - 1 + and iwin1 == self.dims[1] - 1 + and iwin2 == 0 + ): + ywins[-1, -1, 0] = ( + self.taps[-1, -1, 0] * ywins[-1, -1, 0] + ) + elif iwin0 == self.dims[0] - 1 and iwin1 == 0: + ywins[-1, 0, iwin2] = ( + self.taps[-1, 0, 1] * ywins[-1, 0, iwin2] + ) + elif ( + iwin0 == self.dims[0] - 1 and iwin1 == self.dims[1] - 1 + ): + ywins[-1, -1, iwin2] = ( + self.taps[-1, -1, 1] * ywins[-1, -1, iwin2] + ) + elif iwin0 == self.dims[0] - 1 and iwin2 == 0: + ywins[-1, iwin1, 0] = ( + self.taps[-1, 1, 0] * ywins[-1, iwin1, 0] + ) + elif ( + iwin0 == self.dims[0] - 1 and iwin2 == self.dims[2] - 1 + ): + ywins[-1, iwin1, -1] = ( + self.taps[-1, 1, -1] * ywins[-1, iwin1, -1] + ) + elif iwin0 == self.dims[0] - 1: + ywins[-1, iwin1, iwin2] = ( + self.taps[-1, 1, 1] * ywins[-1, iwin1, iwin2] + ) + + elif iwin1 == 0 and iwin2 == 0: + ywins[iwin0, 0, 0] = ( + self.taps[1, 0, 0] * ywins[iwin0, 0, 0] + ) + elif iwin1 == 0 and iwin2 == self.dims[2] - 1: + ywins[iwin0, 0, -1] = ( + self.taps[1, 0, -1] * ywins[iwin0, 0, -1] + ) + elif ( + iwin1 == self.dims[1] - 1 and iwin2 == self.dims[2] - 1 + ): + ywins[iwin0, -1, -1] = ( + self.taps[1, -1, -1] * ywins[iwin0, -1, -1] + ) + elif iwin1 == self.dims[1] - 1 and iwin2 == 0: + ywins[iwin0, -1, 0] = ( + self.taps[1, -1, 0] * ywins[iwin0, -1, 0] + ) + elif iwin1 == 0: + ywins[iwin0, 0, iwin2] = ( + self.taps[1, 0, 1] * ywins[iwin0, 0, iwin2] + ) + elif iwin1 == self.dims[1] - 1: + ywins[iwin0, -1, iwin2] = ( + self.taps[1, -1, 1] * ywins[iwin0, -1, iwin2] + ) + elif iwin2 == 0: + ywins[iwin0, iwin1, 0] = ( + self.taps[1, 1, 0] * ywins[iwin0, iwin1, 0] + ) + elif iwin2 == self.dims[2] - 1: + ywins[iwin0, iwin1, -1] = ( + self.taps[1, 1, -1] * ywins[iwin0, iwin1, -1] + ) + else: + ywins[iwin0, iwin1, iwin2] = ( + self.taps[1, 1, 1] * ywins[iwin0, iwin1, iwin2] + ) + y = self.Op.H @ ywins + else: + y = ncp.zeros(self.dims, dtype=self.dtype) + for iwin0 in range(self.dims[0]): + for iwin1 in range(self.dims[1]): + for iwin2 in range(self.dims[2]): + if self.tapertype is not None: + if iwin0 == 0 and iwin1 == 0 and iwin2 == 0: + ywins[0, 0, 0] = self.taps[0, 0, 0] * ywins[0, 0, 0] + elif ( + iwin0 == 0 and iwin1 == 0 and iwin2 == self.dims[2] - 1 + ): + ywins[0, 0, -1] = self.taps[0, 0, -1] * ywins[0, 0, -1] + elif ( + iwin0 == 0 + and iwin1 == self.dims[1] - 1 + and iwin2 == self.dims[2] - 1 + ): + ywins[0, -1, -1] = ( + self.taps[0, -1, -1] * ywins[0, -1, -1] + ) + elif ( + iwin0 == 0 and iwin1 == self.dims[1] - 1 and iwin2 == 0 + ): + ywins[0, -1, 0] = self.taps[0, -1, 0] * ywins[0, -1, 0] + elif iwin0 == 0 and iwin1 == 0: + ywins[0, 0, iwin2] = ( + self.taps[0, 0, 1] * ywins[0, 0, iwin2] + ) + elif iwin0 == 0 and iwin1 == self.dims[1] - 1: + ywins[0, -1, iwin2] = ( + self.taps[0, -1, 1] * ywins[0, -1, iwin2] + ) + elif iwin0 == 0 and iwin2 == 0: + ywins[0, iwin1, 0] = ( + self.taps[0, 1, 0] * ywins[0, iwin1, 0] + ) + elif iwin0 == 0 and iwin2 == self.dims[2] - 1: + ywins[0, iwin1, -1] = ( + self.taps[0, 1, -1] * ywins[0, iwin1, -1] + ) + elif iwin0 == 0: + ywins[0, iwin1, iwin2] = ( + self.taps[0, 1, 1] * ywins[0, iwin1, iwin2] + ) + + elif ( + iwin0 == self.dims[0] - 1 and iwin1 == 0 and iwin2 == 0 + ): + ywins[-1, 0, 0] = self.taps[-1, 0, 0] * ywins[-1, 0, 0] + elif ( + iwin0 == self.dims[0] - 1 + and iwin1 == 0 + and iwin2 == self.dims[2] - 1 + ): + ywins[-1, 0, -1] = ( + self.taps[-1, 0, -1] * ywins[-1, 0, -1] + ) + elif ( + iwin0 == self.dims[0] - 1 + and iwin1 == self.dims[1] - 1 + and iwin2 == self.dims[2] - 1 + ): + ywins[-1, -1, -1] = ( + self.taps[-1, -1, -1] * ywins[-1, -1, -1] + ) + elif ( + iwin0 == self.dims[0] - 1 + and iwin1 == self.dims[1] - 1 + and iwin2 == 0 + ): + ywins[-1, -1, 0] = ( + self.taps[-1, -1, 0] * ywins[-1, -1, 0] + ) + elif iwin0 == self.dims[0] - 1 and iwin1 == 0: + ywins[-1, 0, iwin2] = ( + self.taps[-1, 0, 1] * ywins[-1, 0, iwin2] + ) + elif ( + iwin0 == self.dims[0] - 1 and iwin1 == self.dims[1] - 1 + ): + ywins[-1, -1, iwin2] = ( + self.taps[-1, -1, 1] * ywins[-1, -1, iwin2] + ) + elif iwin0 == self.dims[0] - 1 and iwin2 == 0: + ywins[-1, iwin1, 0] = ( + self.taps[-1, 1, 0] * ywins[-1, iwin1, 0] + ) + elif ( + iwin0 == self.dims[0] - 1 and iwin2 == self.dims[2] - 1 + ): + ywins[-1, iwin1, -1] = ( + self.taps[-1, 1, -1] * ywins[-1, iwin1, -1] + ) + elif iwin0 == self.dims[0] - 1: + ywins[-1, iwin1, iwin2] = ( + self.taps[-1, 1, 1] * ywins[-1, iwin1, iwin2] + ) + + elif iwin1 == 0 and iwin2 == 0: + ywins[iwin0, 0, 0] = ( + self.taps[1, 0, 0] * ywins[iwin0, 0, 0] + ) + elif iwin1 == 0 and iwin2 == self.dims[2] - 1: + ywins[iwin0, 0, -1] = ( + self.taps[1, 0, -1] * ywins[iwin0, 0, -1] + ) + elif ( + iwin1 == self.dims[1] - 1 and iwin2 == self.dims[2] - 1 + ): + ywins[iwin0, -1, -1] = ( + self.taps[1, -1, -1] * ywins[iwin0, -1, -1] + ) + elif iwin1 == self.dims[1] - 1 and iwin2 == 0: + ywins[iwin0, -1, 0] = ( + self.taps[1, -1, 0] * ywins[iwin0, -1, 0] + ) + elif iwin1 == 0: + ywins[iwin0, 0, iwin2] = ( + self.taps[1, 0, 1] * ywins[iwin0, 0, iwin2] + ) + elif iwin1 == self.dims[1] - 1: + ywins[iwin0, -1, iwin2] = ( + self.taps[1, -1, 1] * ywins[iwin0, -1, iwin2] + ) + elif iwin2 == 0: + ywins[iwin0, iwin1, 0] = ( + self.taps[1, 1, 0] * ywins[iwin0, iwin1, 0] + ) + elif iwin2 == self.dims[2] - 1: + ywins[iwin0, iwin1, -1] = ( + self.taps[1, 1, -1] * ywins[iwin0, iwin1, -1] + ) + else: + ywins[iwin0, iwin1, iwin2] = ( + self.taps[1, 1, 1] * ywins[iwin0, iwin1, iwin2] + ) + y[iwin0, iwin1, iwin2] = self.Op.rmatvec( + ywins[iwin0, iwin1, iwin2].ravel() + ).reshape(self.dims[3], self.dims[4], self.dims[5]) + return y + + def _register_multiplications(self, savetaper: bool) -> None: + if savetaper: + self._matvec = self._matvec_savetaper + self._rmatvec = self._rmatvec_savetaper + else: + self._matvec = self._matvec_nosavetaper + self._rmatvec = self._rmatvec_nosavetaper diff --git a/pytests/test_patching.py b/pytests/test_patching.py index 25656724..de61243e 100755 --- a/pytests/test_patching.py +++ b/pytests/test_patching.py @@ -25,6 +25,7 @@ "novert": 0, # "winst": 2, "tapertype": None, + "savetaper": True, } # no overlap, no taper par2 = { "ny": 6, @@ -43,6 +44,7 @@ "novert": 0, # "winst": 2, "tapertype": "hanning", + "savetaper": True, } # no overlap, with taper par3 = { "ny": 6, @@ -61,8 +63,28 @@ "novert": 2, # "winst": 4, "tapertype": None, + "savetaper": True, } # overlap, no taper par4 = { + "ny": 6, + "nx": 7, + "nt": 10, + "npy": 15, + "nwiny": 7, + "novery": 3, + # "winsy": 3, + "npx": 13, + "nwinx": 5, + "noverx": 2, + # "winsx": 3, + "npt": 10, + "nwint": 4, + "novert": 2, + # "winst": 4, + "tapertype": None, + "savetaper": False, +} # overlap, no taper (non saved +par5 = { "ny": 6, "nx": 7, "nt": 10, @@ -79,10 +101,30 @@ "novert": 2, # "winst": 4, "tapertype": "hanning", + "savetaper": True, } # overlap, with taper +par6 = { + "ny": 6, + "nx": 7, + "nt": 10, + "npy": 15, + "nwiny": 7, + "novery": 3, + # "winsy": 3, + "npx": 13, + "nwinx": 5, + "noverx": 2, + # "winsx": 3, + "npt": 10, + "nwint": 4, + "novert": 2, + # "winst": 4, + "tapertype": "hanning", + "savetaper": False, +} # overlap, with taper (non saved) -@pytest.mark.parametrize("par", [(par1), (par2), (par3), (par4)]) +@pytest.mark.parametrize("par", [(par1), (par2), (par3), (par4), (par5), (par6)]) def test_Patch2D(par): """Dot-test and inverse for Patch2D operator""" Op = MatrixMult(np.ones((par["nwiny"] * par["nwint"], par["ny"] * par["nt"]))) @@ -101,6 +143,7 @@ def test_Patch2D(par): nover=(par["novery"], par["novert"]), nop=(par["ny"], par["nt"]), tapertype=par["tapertype"], + savetaper=par["savetaper"], ) assert dottest( Pop, @@ -134,6 +177,7 @@ def test_Patch2D_scalings(par): nover=(par["novery"], par["novert"]), nop=(par["ny"], par["nt"]), tapertype=par["tapertype"], + savetaper=par["savetaper"], scalings=scalings, ) assert dottest( @@ -148,7 +192,7 @@ def test_Patch2D_scalings(par): assert_array_almost_equal(x.ravel(), xinv) -@pytest.mark.parametrize("par", [(par1), (par2), (par3), (par4)]) +@pytest.mark.parametrize("par", [(par1), (par2), (par3), (par4), (par5), (par6)]) def test_Patch3D(par): """Dot-test and inverse for Patch3D operator""" Op = MatrixMult( @@ -179,6 +223,7 @@ def test_Patch3D(par): nover=(par["novery"], par["noverx"], par["novert"]), nop=(par["ny"], par["nx"], par["nt"]), tapertype=par["tapertype"], + savetaper=par["savetaper"], ) assert dottest( Pop, From 6c57fc33ad8fe18bbb342eabfa99baa4f20ed51a Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sat, 18 May 2024 20:59:39 +0300 Subject: [PATCH 08/11] minor: added versionadded to docstrings --- pylops/signalprocessing/patch2d.py | 2 ++ pylops/signalprocessing/patch3d.py | 2 ++ pylops/signalprocessing/sliding1d.py | 2 ++ pylops/signalprocessing/sliding2d.py | 2 ++ pylops/signalprocessing/sliding3d.py | 2 ++ 5 files changed, 10 insertions(+) diff --git a/pylops/signalprocessing/patch2d.py b/pylops/signalprocessing/patch2d.py index e60fd807..1a352303 100644 --- a/pylops/signalprocessing/patch2d.py +++ b/pylops/signalprocessing/patch2d.py @@ -142,6 +142,8 @@ class Patch2D(LinearOperator): tapertype : :obj:`str`, optional Type of taper (``hanning``, ``cosine``, ``cosinesquare`` or ``None``) savetaper: :obj:`bool`, optional + .. versionadded:: 2.3.0 + Save all tapers and apply them in one go (``True``) or save unique tapers and apply them one by one (``False``). The first option is more computationally efficient, whilst the second is more memory efficient. scalings : :obj:`tuple` or :obj:`list`, optional diff --git a/pylops/signalprocessing/patch3d.py b/pylops/signalprocessing/patch3d.py index 7a893b54..7b969984 100644 --- a/pylops/signalprocessing/patch3d.py +++ b/pylops/signalprocessing/patch3d.py @@ -157,6 +157,8 @@ class Patch3D(LinearOperator): tapertype : :obj:`str`, optional Type of taper (``hanning``, ``cosine``, ``cosinesquare`` or ``None``) savetaper: :obj:`bool`, optional + .. versionadded:: 2.3.0 + Save all tapers and apply them in one go (``True``) or save unique tapers and apply them one by one (``False``). The first option is more computationally efficient, whilst the second is more memory efficient. scalings : :obj:`tuple` or :obj:`list`, optional diff --git a/pylops/signalprocessing/sliding1d.py b/pylops/signalprocessing/sliding1d.py index 6e0fb26b..767045d1 100644 --- a/pylops/signalprocessing/sliding1d.py +++ b/pylops/signalprocessing/sliding1d.py @@ -128,6 +128,8 @@ class Sliding1D(LinearOperator): tapertype : :obj:`str`, optional Type of taper (``hanning``, ``cosine``, ``cosinesquare`` or ``None``) savetaper: :obj:`bool`, optional + .. versionadded:: 2.3.0 + Save all tapers and apply them in one go (``True``) or save unique tapers and apply them one by one (``False``). The first option is more computationally efficient, whilst the second is more memory efficient. name : :obj:`str`, optional diff --git a/pylops/signalprocessing/sliding2d.py b/pylops/signalprocessing/sliding2d.py index fa930662..f3b3cfb3 100644 --- a/pylops/signalprocessing/sliding2d.py +++ b/pylops/signalprocessing/sliding2d.py @@ -164,6 +164,8 @@ class Sliding2D(LinearOperator): tapertype : :obj:`str`, optional Type of taper (``hanning``, ``cosine``, ``cosinesquare`` or ``None``) savetaper: :obj:`bool`, optional + .. versionadded:: 2.3.0 + Save all tapers and apply them in one go (``True``) or save unique tapers and apply them one by one (``False``). The first option is more computationally efficient, whilst the second is more memory efficient. name : :obj:`str`, optional diff --git a/pylops/signalprocessing/sliding3d.py b/pylops/signalprocessing/sliding3d.py index 846a678e..9dc7e255 100644 --- a/pylops/signalprocessing/sliding3d.py +++ b/pylops/signalprocessing/sliding3d.py @@ -150,6 +150,8 @@ class Sliding3D(LinearOperator): tapertype : :obj:`str`, optional Type of taper (``hanning``, ``cosine``, ``cosinesquare`` or ``None``) savetaper: :obj:`bool`, optional + .. versionadded:: 2.3.0 + Save all tapers and apply them in one go (``True``) or save unique tapers and apply them one by one (``False``). The first option is more computationally efficient, whilst the second is more memory efficient. nproc : :obj:`int`, optional From cd3aadfc54696487810260a0f0b764510bbb9a63 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sat, 18 May 2024 21:30:18 +0300 Subject: [PATCH 09/11] minor: moved repeated code into private method --- pylops/signalprocessing/patch2d.py | 59 ++--- pylops/signalprocessing/patch3d.py | 325 ++++++--------------------- pylops/signalprocessing/sliding1d.py | 23 +- pylops/signalprocessing/sliding2d.py | 23 +- pylops/signalprocessing/sliding3d.py | 59 ++--- 5 files changed, 133 insertions(+), 356 deletions(-) diff --git a/pylops/signalprocessing/patch2d.py b/pylops/signalprocessing/patch2d.py index 1a352303..3aa157a5 100644 --- a/pylops/signalprocessing/patch2d.py +++ b/pylops/signalprocessing/patch2d.py @@ -297,6 +297,27 @@ def __init__( self._register_multiplications(self.savetaper) + def _apply_taper(self, ywins, iwin0, iwin1): + if iwin0 == 0 and iwin1 == 0: + ywins[0, 0] = self.taps[0, 0] * ywins[0, 0] + elif iwin0 == 0 and iwin1 == self.dims[1] - 1: + ywins[0, -1] = self.taps[0, -1] * ywins[0, -1] + elif iwin0 == 0: + ywins[0, iwin1] = self.taps[0, 1] * ywins[0, iwin1] + elif iwin0 == self.dims[0] - 1 and iwin1 == 0: + ywins[-1, 0] = self.taps[-1, 0] * ywins[-1, 0] + elif iwin0 == self.dims[0] - 1 and iwin1 == self.dims[1] - 1: + ywins[-1, -1] = self.taps[-1, -1] * ywins[-1, -1] + elif iwin0 == self.dims[0] - 1: + ywins[-1, iwin1] = self.taps[-1, 1] * ywins[-1, iwin1] + elif iwin1 == 0: + ywins[iwin0, 0] = self.taps[1, 0] * ywins[iwin0, 0] + elif iwin1 == self.dims[1] - 1: + ywins[iwin0, -1] = self.taps[1, -1] * ywins[iwin0, -1] + else: + ywins[iwin0, iwin1] = self.taps[1, 1] * ywins[iwin0, iwin1] + return ywins + @reshaped() def _matvec_savetaper(self, x: NDArray) -> NDArray: ncp = get_array_module(x) @@ -397,48 +418,14 @@ def _rmatvec_nosavetaper(self, x: NDArray) -> NDArray: if self.tapertype is not None: for iwin0 in range(self.dims[0]): for iwin1 in range(self.dims[1]): - if iwin0 == 0 and iwin1 == 0: - ywins[0, 0] = self.taps[0, 0] * ywins[0, 0] - elif iwin0 == 0 and iwin1 == self.dims[1] - 1: - ywins[0, -1] = self.taps[0, -1] * ywins[0, -1] - elif iwin0 == 0: - ywins[0, iwin1] = self.taps[0, 1] * ywins[0, iwin1] - elif iwin0 == self.dims[0] - 1 and iwin1 == 0: - ywins[-1, 0] = self.taps[-1, 0] * ywins[-1, 0] - elif iwin0 == self.dims[0] - 1 and iwin1 == self.dims[1] - 1: - ywins[-1, -1] = self.taps[-1, -1] * ywins[-1, -1] - elif iwin0 == self.dims[0] - 1: - ywins[-1, iwin1] = self.taps[-1, 1] * ywins[-1, iwin1] - elif iwin1 == 0: - ywins[iwin0, 0] = self.taps[1, 0] * ywins[iwin0, 0] - elif iwin1 == self.dims[1] - 1: - ywins[iwin0, -1] = self.taps[1, -1] * ywins[iwin0, -1] - else: - ywins[iwin0, iwin1] = self.taps[1, 1] * ywins[iwin0, iwin1] + ywins = self._apply_taper(ywins, iwin0, iwin1) y = self.Op.H @ ywins else: y = ncp.zeros(self.dims, dtype=self.dtype) for iwin0 in range(self.dims[0]): for iwin1 in range(self.dims[1]): if self.tapertype is not None: - if iwin0 == 0 and iwin1 == 0: - ywins[0, 0] = self.taps[0, 0] * ywins[0, 0] - elif iwin0 == 0 and iwin1 == self.dims[1] - 1: - ywins[0, -1] = self.taps[0, -1] * ywins[0, -1] - elif iwin0 == 0: - ywins[0, iwin1] = self.taps[0, 1] * ywins[0, iwin1] - elif iwin0 == self.dims[0] - 1 and iwin1 == 0: - ywins[-1, 0] = self.taps[-1, 0] * ywins[-1, 0] - elif iwin0 == self.dims[0] - 1 and iwin1 == self.dims[1] - 1: - ywins[-1, -1] = self.taps[-1, -1] * ywins[-1, -1] - elif iwin0 == self.dims[0] - 1: - ywins[-1, iwin1] = self.taps[-1, 1] * ywins[-1, iwin1] - elif iwin1 == 0: - ywins[iwin0, 0] = self.taps[1, 0] * ywins[iwin0, 0] - elif iwin1 == self.dims[1] - 1: - ywins[iwin0, -1] = self.taps[1, -1] * ywins[iwin0, -1] - else: - ywins[iwin0, iwin1] = self.taps[1, 1] * ywins[iwin0, iwin1] + ywins = self._apply_taper(ywins, iwin0, iwin1) y[iwin0, iwin1] = self.Op.rmatvec( ywins[iwin0, iwin1].ravel() ).reshape(self.dims[2], self.dims[3]) diff --git a/pylops/signalprocessing/patch3d.py b/pylops/signalprocessing/patch3d.py index 7b969984..71ffe341 100644 --- a/pylops/signalprocessing/patch3d.py +++ b/pylops/signalprocessing/patch3d.py @@ -513,6 +513,69 @@ def __init__( self._register_multiplications(self.savetaper) + def _apply_taper(self, ywins, iwin0, iwin1, iwin2): + if iwin0 == 0 and iwin1 == 0 and iwin2 == 0: + ywins[0, 0, 0] = self.taps[0, 0, 0] * ywins[0, 0, 0] + elif iwin0 == 0 and iwin1 == 0 and iwin2 == self.dims[2] - 1: + ywins[0, 0, -1] = self.taps[0, 0, -1] * ywins[0, 0, -1] + elif iwin0 == 0 and iwin1 == self.dims[1] - 1 and iwin2 == self.dims[2] - 1: + ywins[0, -1, -1] = self.taps[0, -1, -1] * ywins[0, -1, -1] + elif iwin0 == 0 and iwin1 == self.dims[1] - 1 and iwin2 == 0: + ywins[0, -1, 0] = self.taps[0, -1, 0] * ywins[0, -1, 0] + elif iwin0 == 0 and iwin1 == 0: + ywins[0, 0, iwin2] = self.taps[0, 0, 1] * ywins[0, 0, iwin2] + elif iwin0 == 0 and iwin1 == self.dims[1] - 1: + ywins[0, -1, iwin2] = self.taps[0, -1, 1] * ywins[0, -1, iwin2] + elif iwin0 == 0 and iwin2 == 0: + ywins[0, iwin1, 0] = self.taps[0, 1, 0] * ywins[0, iwin1, 0] + elif iwin0 == 0 and iwin2 == self.dims[2] - 1: + ywins[0, iwin1, -1] = self.taps[0, 1, -1] * ywins[0, iwin1, -1] + elif iwin0 == 0: + ywins[0, iwin1, iwin2] = self.taps[0, 1, 1] * ywins[0, iwin1, iwin2] + + elif iwin0 == self.dims[0] - 1 and iwin1 == 0 and iwin2 == 0: + ywins[-1, 0, 0] = self.taps[-1, 0, 0] * ywins[-1, 0, 0] + elif iwin0 == self.dims[0] - 1 and iwin1 == 0 and iwin2 == self.dims[2] - 1: + ywins[-1, 0, -1] = self.taps[-1, 0, -1] * ywins[-1, 0, -1] + elif ( + iwin0 == self.dims[0] - 1 + and iwin1 == self.dims[1] - 1 + and iwin2 == self.dims[2] - 1 + ): + ywins[-1, -1, -1] = self.taps[-1, -1, -1] * ywins[-1, -1, -1] + elif iwin0 == self.dims[0] - 1 and iwin1 == self.dims[1] - 1 and iwin2 == 0: + ywins[-1, -1, 0] = self.taps[-1, -1, 0] * ywins[-1, -1, 0] + elif iwin0 == self.dims[0] - 1 and iwin1 == 0: + ywins[-1, 0, iwin2] = self.taps[-1, 0, 1] * ywins[-1, 0, iwin2] + elif iwin0 == self.dims[0] - 1 and iwin1 == self.dims[1] - 1: + ywins[-1, -1, iwin2] = self.taps[-1, -1, 1] * ywins[-1, -1, iwin2] + elif iwin0 == self.dims[0] - 1 and iwin2 == 0: + ywins[-1, iwin1, 0] = self.taps[-1, 1, 0] * ywins[-1, iwin1, 0] + elif iwin0 == self.dims[0] - 1 and iwin2 == self.dims[2] - 1: + ywins[-1, iwin1, -1] = self.taps[-1, 1, -1] * ywins[-1, iwin1, -1] + elif iwin0 == self.dims[0] - 1: + ywins[-1, iwin1, iwin2] = self.taps[-1, 1, 1] * ywins[-1, iwin1, iwin2] + + elif iwin1 == 0 and iwin2 == 0: + ywins[iwin0, 0, 0] = self.taps[1, 0, 0] * ywins[iwin0, 0, 0] + elif iwin1 == 0 and iwin2 == self.dims[2] - 1: + ywins[iwin0, 0, -1] = self.taps[1, 0, -1] * ywins[iwin0, 0, -1] + elif iwin1 == self.dims[1] - 1 and iwin2 == self.dims[2] - 1: + ywins[iwin0, -1, -1] = self.taps[1, -1, -1] * ywins[iwin0, -1, -1] + elif iwin1 == self.dims[1] - 1 and iwin2 == 0: + ywins[iwin0, -1, 0] = self.taps[1, -1, 0] * ywins[iwin0, -1, 0] + elif iwin1 == 0: + ywins[iwin0, 0, iwin2] = self.taps[1, 0, 1] * ywins[iwin0, 0, iwin2] + elif iwin1 == self.dims[1] - 1: + ywins[iwin0, -1, iwin2] = self.taps[1, -1, 1] * ywins[iwin0, -1, iwin2] + elif iwin2 == 0: + ywins[iwin0, iwin1, 0] = self.taps[1, 1, 0] * ywins[iwin0, iwin1, 0] + elif iwin2 == self.dims[2] - 1: + ywins[iwin0, iwin1, -1] = self.taps[1, 1, -1] * ywins[iwin0, iwin1, -1] + else: + ywins[iwin0, iwin1, iwin2] = self.taps[1, 1, 1] * ywins[iwin0, iwin1, iwin2] + return ywins + @reshaped() def _matvec_savetaper(self, x: NDArray) -> NDArray: ncp = get_array_module(x) @@ -680,136 +743,7 @@ def _rmatvec_nosavetaper(self, x: NDArray) -> NDArray: for iwin0 in range(self.dims[0]): for iwin1 in range(self.dims[1]): for iwin2 in range(self.dims[2]): - if iwin0 == 0 and iwin1 == 0 and iwin2 == 0: - ywins[0, 0, 0] = self.taps[0, 0, 0] * ywins[0, 0, 0] - elif ( - iwin0 == 0 and iwin1 == 0 and iwin2 == self.dims[2] - 1 - ): - ywins[0, 0, -1] = self.taps[0, 0, -1] * ywins[0, 0, -1] - elif ( - iwin0 == 0 - and iwin1 == self.dims[1] - 1 - and iwin2 == self.dims[2] - 1 - ): - ywins[0, -1, -1] = ( - self.taps[0, -1, -1] * ywins[0, -1, -1] - ) - elif ( - iwin0 == 0 and iwin1 == self.dims[1] - 1 and iwin2 == 0 - ): - ywins[0, -1, 0] = self.taps[0, -1, 0] * ywins[0, -1, 0] - elif iwin0 == 0 and iwin1 == 0: - ywins[0, 0, iwin2] = ( - self.taps[0, 0, 1] * ywins[0, 0, iwin2] - ) - elif iwin0 == 0 and iwin1 == self.dims[1] - 1: - ywins[0, -1, iwin2] = ( - self.taps[0, -1, 1] * ywins[0, -1, iwin2] - ) - elif iwin0 == 0 and iwin2 == 0: - ywins[0, iwin1, 0] = ( - self.taps[0, 1, 0] * ywins[0, iwin1, 0] - ) - elif iwin0 == 0 and iwin2 == self.dims[2] - 1: - ywins[0, iwin1, -1] = ( - self.taps[0, 1, -1] * ywins[0, iwin1, -1] - ) - elif iwin0 == 0: - ywins[0, iwin1, iwin2] = ( - self.taps[0, 1, 1] * ywins[0, iwin1, iwin2] - ) - - elif ( - iwin0 == self.dims[0] - 1 and iwin1 == 0 and iwin2 == 0 - ): - ywins[-1, 0, 0] = self.taps[-1, 0, 0] * ywins[-1, 0, 0] - elif ( - iwin0 == self.dims[0] - 1 - and iwin1 == 0 - and iwin2 == self.dims[2] - 1 - ): - ywins[-1, 0, -1] = ( - self.taps[-1, 0, -1] * ywins[-1, 0, -1] - ) - elif ( - iwin0 == self.dims[0] - 1 - and iwin1 == self.dims[1] - 1 - and iwin2 == self.dims[2] - 1 - ): - ywins[-1, -1, -1] = ( - self.taps[-1, -1, -1] * ywins[-1, -1, -1] - ) - elif ( - iwin0 == self.dims[0] - 1 - and iwin1 == self.dims[1] - 1 - and iwin2 == 0 - ): - ywins[-1, -1, 0] = ( - self.taps[-1, -1, 0] * ywins[-1, -1, 0] - ) - elif iwin0 == self.dims[0] - 1 and iwin1 == 0: - ywins[-1, 0, iwin2] = ( - self.taps[-1, 0, 1] * ywins[-1, 0, iwin2] - ) - elif ( - iwin0 == self.dims[0] - 1 and iwin1 == self.dims[1] - 1 - ): - ywins[-1, -1, iwin2] = ( - self.taps[-1, -1, 1] * ywins[-1, -1, iwin2] - ) - elif iwin0 == self.dims[0] - 1 and iwin2 == 0: - ywins[-1, iwin1, 0] = ( - self.taps[-1, 1, 0] * ywins[-1, iwin1, 0] - ) - elif ( - iwin0 == self.dims[0] - 1 and iwin2 == self.dims[2] - 1 - ): - ywins[-1, iwin1, -1] = ( - self.taps[-1, 1, -1] * ywins[-1, iwin1, -1] - ) - elif iwin0 == self.dims[0] - 1: - ywins[-1, iwin1, iwin2] = ( - self.taps[-1, 1, 1] * ywins[-1, iwin1, iwin2] - ) - - elif iwin1 == 0 and iwin2 == 0: - ywins[iwin0, 0, 0] = ( - self.taps[1, 0, 0] * ywins[iwin0, 0, 0] - ) - elif iwin1 == 0 and iwin2 == self.dims[2] - 1: - ywins[iwin0, 0, -1] = ( - self.taps[1, 0, -1] * ywins[iwin0, 0, -1] - ) - elif ( - iwin1 == self.dims[1] - 1 and iwin2 == self.dims[2] - 1 - ): - ywins[iwin0, -1, -1] = ( - self.taps[1, -1, -1] * ywins[iwin0, -1, -1] - ) - elif iwin1 == self.dims[1] - 1 and iwin2 == 0: - ywins[iwin0, -1, 0] = ( - self.taps[1, -1, 0] * ywins[iwin0, -1, 0] - ) - elif iwin1 == 0: - ywins[iwin0, 0, iwin2] = ( - self.taps[1, 0, 1] * ywins[iwin0, 0, iwin2] - ) - elif iwin1 == self.dims[1] - 1: - ywins[iwin0, -1, iwin2] = ( - self.taps[1, -1, 1] * ywins[iwin0, -1, iwin2] - ) - elif iwin2 == 0: - ywins[iwin0, iwin1, 0] = ( - self.taps[1, 1, 0] * ywins[iwin0, iwin1, 0] - ) - elif iwin2 == self.dims[2] - 1: - ywins[iwin0, iwin1, -1] = ( - self.taps[1, 1, -1] * ywins[iwin0, iwin1, -1] - ) - else: - ywins[iwin0, iwin1, iwin2] = ( - self.taps[1, 1, 1] * ywins[iwin0, iwin1, iwin2] - ) + ywins = self._apply_taper(ywins, iwin0, iwin1, iwin2) y = self.Op.H @ ywins else: y = ncp.zeros(self.dims, dtype=self.dtype) @@ -817,136 +751,7 @@ def _rmatvec_nosavetaper(self, x: NDArray) -> NDArray: for iwin1 in range(self.dims[1]): for iwin2 in range(self.dims[2]): if self.tapertype is not None: - if iwin0 == 0 and iwin1 == 0 and iwin2 == 0: - ywins[0, 0, 0] = self.taps[0, 0, 0] * ywins[0, 0, 0] - elif ( - iwin0 == 0 and iwin1 == 0 and iwin2 == self.dims[2] - 1 - ): - ywins[0, 0, -1] = self.taps[0, 0, -1] * ywins[0, 0, -1] - elif ( - iwin0 == 0 - and iwin1 == self.dims[1] - 1 - and iwin2 == self.dims[2] - 1 - ): - ywins[0, -1, -1] = ( - self.taps[0, -1, -1] * ywins[0, -1, -1] - ) - elif ( - iwin0 == 0 and iwin1 == self.dims[1] - 1 and iwin2 == 0 - ): - ywins[0, -1, 0] = self.taps[0, -1, 0] * ywins[0, -1, 0] - elif iwin0 == 0 and iwin1 == 0: - ywins[0, 0, iwin2] = ( - self.taps[0, 0, 1] * ywins[0, 0, iwin2] - ) - elif iwin0 == 0 and iwin1 == self.dims[1] - 1: - ywins[0, -1, iwin2] = ( - self.taps[0, -1, 1] * ywins[0, -1, iwin2] - ) - elif iwin0 == 0 and iwin2 == 0: - ywins[0, iwin1, 0] = ( - self.taps[0, 1, 0] * ywins[0, iwin1, 0] - ) - elif iwin0 == 0 and iwin2 == self.dims[2] - 1: - ywins[0, iwin1, -1] = ( - self.taps[0, 1, -1] * ywins[0, iwin1, -1] - ) - elif iwin0 == 0: - ywins[0, iwin1, iwin2] = ( - self.taps[0, 1, 1] * ywins[0, iwin1, iwin2] - ) - - elif ( - iwin0 == self.dims[0] - 1 and iwin1 == 0 and iwin2 == 0 - ): - ywins[-1, 0, 0] = self.taps[-1, 0, 0] * ywins[-1, 0, 0] - elif ( - iwin0 == self.dims[0] - 1 - and iwin1 == 0 - and iwin2 == self.dims[2] - 1 - ): - ywins[-1, 0, -1] = ( - self.taps[-1, 0, -1] * ywins[-1, 0, -1] - ) - elif ( - iwin0 == self.dims[0] - 1 - and iwin1 == self.dims[1] - 1 - and iwin2 == self.dims[2] - 1 - ): - ywins[-1, -1, -1] = ( - self.taps[-1, -1, -1] * ywins[-1, -1, -1] - ) - elif ( - iwin0 == self.dims[0] - 1 - and iwin1 == self.dims[1] - 1 - and iwin2 == 0 - ): - ywins[-1, -1, 0] = ( - self.taps[-1, -1, 0] * ywins[-1, -1, 0] - ) - elif iwin0 == self.dims[0] - 1 and iwin1 == 0: - ywins[-1, 0, iwin2] = ( - self.taps[-1, 0, 1] * ywins[-1, 0, iwin2] - ) - elif ( - iwin0 == self.dims[0] - 1 and iwin1 == self.dims[1] - 1 - ): - ywins[-1, -1, iwin2] = ( - self.taps[-1, -1, 1] * ywins[-1, -1, iwin2] - ) - elif iwin0 == self.dims[0] - 1 and iwin2 == 0: - ywins[-1, iwin1, 0] = ( - self.taps[-1, 1, 0] * ywins[-1, iwin1, 0] - ) - elif ( - iwin0 == self.dims[0] - 1 and iwin2 == self.dims[2] - 1 - ): - ywins[-1, iwin1, -1] = ( - self.taps[-1, 1, -1] * ywins[-1, iwin1, -1] - ) - elif iwin0 == self.dims[0] - 1: - ywins[-1, iwin1, iwin2] = ( - self.taps[-1, 1, 1] * ywins[-1, iwin1, iwin2] - ) - - elif iwin1 == 0 and iwin2 == 0: - ywins[iwin0, 0, 0] = ( - self.taps[1, 0, 0] * ywins[iwin0, 0, 0] - ) - elif iwin1 == 0 and iwin2 == self.dims[2] - 1: - ywins[iwin0, 0, -1] = ( - self.taps[1, 0, -1] * ywins[iwin0, 0, -1] - ) - elif ( - iwin1 == self.dims[1] - 1 and iwin2 == self.dims[2] - 1 - ): - ywins[iwin0, -1, -1] = ( - self.taps[1, -1, -1] * ywins[iwin0, -1, -1] - ) - elif iwin1 == self.dims[1] - 1 and iwin2 == 0: - ywins[iwin0, -1, 0] = ( - self.taps[1, -1, 0] * ywins[iwin0, -1, 0] - ) - elif iwin1 == 0: - ywins[iwin0, 0, iwin2] = ( - self.taps[1, 0, 1] * ywins[iwin0, 0, iwin2] - ) - elif iwin1 == self.dims[1] - 1: - ywins[iwin0, -1, iwin2] = ( - self.taps[1, -1, 1] * ywins[iwin0, -1, iwin2] - ) - elif iwin2 == 0: - ywins[iwin0, iwin1, 0] = ( - self.taps[1, 1, 0] * ywins[iwin0, iwin1, 0] - ) - elif iwin2 == self.dims[2] - 1: - ywins[iwin0, iwin1, -1] = ( - self.taps[1, 1, -1] * ywins[iwin0, iwin1, -1] - ) - else: - ywins[iwin0, iwin1, iwin2] = ( - self.taps[1, 1, 1] * ywins[iwin0, iwin1, iwin2] - ) + ywins = self._apply_taper(ywins, iwin0, iwin1, iwin2) y[iwin0, iwin1, iwin2] = self.Op.rmatvec( ywins[iwin0, iwin1, iwin2].ravel() ).reshape(self.dims[3], self.dims[4], self.dims[5]) diff --git a/pylops/signalprocessing/sliding1d.py b/pylops/signalprocessing/sliding1d.py index 767045d1..9bb331bb 100644 --- a/pylops/signalprocessing/sliding1d.py +++ b/pylops/signalprocessing/sliding1d.py @@ -212,6 +212,15 @@ def __init__( self._register_multiplications(self.savetaper) + def _apply_taper(self, ywins, iwin0): + if iwin0 == 0: + ywins[0] = ywins[0] * self.taps[0] + elif iwin0 == self.dims[0] - 1: + ywins[-1] = ywins[-1] * self.taps[-1] + else: + ywins[iwin0] = ywins[iwin0] * self.taps[1] + return ywins + @reshaped def _matvec_savetaper(self, x: NDArray) -> NDArray: ncp = get_array_module(x) @@ -282,23 +291,13 @@ def _rmatvec_nosavetaper(self, x: NDArray) -> NDArray: if self.simOp: if self.tapertype is not None: for iwin0 in range(self.dims[0]): - if iwin0 == 0: - ywins[0] = ywins[0] * self.taps[0] - elif iwin0 == self.dims[0] - 1: - ywins[-1] = ywins[-1] * self.taps[-1] - else: - ywins[iwin0] = ywins[iwin0] * self.taps[1] + ywins = self._apply_taper(ywins, iwin0) y = self.Op.H @ ywins else: y = ncp.zeros(self.dims, dtype=self.dtype) for iwin0 in range(self.dims[0]): if self.tapertype is not None: - if iwin0 == 0: - ywins[0] = ywins[0] * self.taps[0] - elif iwin0 == self.dims[0] - 1: - ywins[-1] = ywins[-1] * self.taps[-1] - else: - ywins[iwin0] = ywins[iwin0] * self.taps[1] + ywins = self._apply_taper(ywins, iwin0) y[iwin0] = self.Op.rmatvec(ywins[iwin0]) return y diff --git a/pylops/signalprocessing/sliding2d.py b/pylops/signalprocessing/sliding2d.py index f3b3cfb3..4a76f5e6 100644 --- a/pylops/signalprocessing/sliding2d.py +++ b/pylops/signalprocessing/sliding2d.py @@ -255,6 +255,15 @@ def __init__( self._register_multiplications(self.savetaper) + def _apply_taper(self, ywins, iwin0): + if iwin0 == 0: + ywins[0] = ywins[0] * self.taps[0] + elif iwin0 == self.dims[0] - 1: + ywins[-1] = ywins[-1] * self.taps[-1] + else: + ywins[iwin0] = ywins[iwin0] * self.taps[1] + return ywins + @reshaped def _matvec_savetaper(self, x: NDArray) -> NDArray: ncp = get_array_module(x) @@ -335,23 +344,13 @@ def _rmatvec_nosavetaper(self, x: NDArray) -> NDArray: if self.simOp: if self.tapertype is not None: for iwin0 in range(self.dims[0]): - if iwin0 == 0: - ywins[0] = ywins[0] * self.taps[0] - elif iwin0 == self.dims[0] - 1: - ywins[-1] = ywins[-1] * self.taps[-1] - else: - ywins[iwin0] = ywins[iwin0] * self.taps[1] + ywins = self._apply_taper(ywins, iwin0) y = self.Op.H @ ywins else: y = ncp.zeros(self.dims, dtype=self.dtype) for iwin0 in range(self.dims[0]): if self.tapertype is not None: - if iwin0 == 0: - ywins[0] = ywins[0] * self.taps[0] - elif iwin0 == self.dims[0] - 1: - ywins[-1] = ywins[-1] * self.taps[-1] - else: - ywins[iwin0] = ywins[iwin0] * self.taps[1] + ywins = self._apply_taper(ywins, iwin0) y[iwin0] = self.Op.rmatvec(ywins[iwin0].ravel()).reshape( self.dims[1], self.dims[2] ) diff --git a/pylops/signalprocessing/sliding3d.py b/pylops/signalprocessing/sliding3d.py index 9dc7e255..49324230 100644 --- a/pylops/signalprocessing/sliding3d.py +++ b/pylops/signalprocessing/sliding3d.py @@ -307,6 +307,27 @@ def __init__( self._register_multiplications(self.savetaper) + def _apply_taper(self, ywins, iwin0, iwin1): + if iwin0 == 0 and iwin1 == 0: + ywins[0, 0] = self.taps[0, 0] * ywins[0, 0] + elif iwin0 == 0 and iwin1 == self.dims[1] - 1: + ywins[0, -1] = self.taps[0, -1] * ywins[0, -1] + elif iwin0 == 0: + ywins[0, iwin1] = self.taps[0, 1] * ywins[0, iwin1] + elif iwin0 == self.dims[0] - 1 and iwin1 == 0: + ywins[-1, 0] = self.taps[-1, 0] * ywins[-1, 0] + elif iwin0 == self.dims[0] - 1 and iwin1 == self.dims[1] - 1: + ywins[-1, -1] = self.taps[-1, -1] * ywins[-1, -1] + elif iwin0 == self.dims[0] - 1: + ywins[-1, iwin1] = self.taps[-1, 1] * ywins[-1, iwin1] + elif iwin1 == 0: + ywins[iwin0, 0] = self.taps[1, 0] * ywins[iwin0, 0] + elif iwin1 == self.dims[1] - 1: + ywins[iwin0, -1] = self.taps[1, -1] * ywins[iwin0, -1] + else: + ywins[iwin0, iwin1] = self.taps[1, 1] * ywins[iwin0, iwin1] + return ywins + @reshaped def _matvec_savetaper(self, x: NDArray) -> NDArray: ncp = get_array_module(x) @@ -417,48 +438,14 @@ def _rmatvec_nosavetaper(self, x: NDArray) -> NDArray: if self.tapertype is not None: for iwin0 in range(self.dims[0]): for iwin1 in range(self.dims[1]): - if iwin0 == 0 and iwin1 == 0: - ywins[0, 0] = self.taps[0, 0] * ywins[0, 0] - elif iwin0 == 0 and iwin1 == self.dims[1] - 1: - ywins[0, -1] = self.taps[0, -1] * ywins[0, -1] - elif iwin0 == 0: - ywins[0, iwin1] = self.taps[0, 1] * ywins[0, iwin1] - elif iwin0 == self.dims[0] - 1 and iwin1 == 0: - ywins[-1, 0] = self.taps[-1, 0] * ywins[-1, 0] - elif iwin0 == self.dims[0] - 1 and iwin1 == self.dims[1] - 1: - ywins[-1, -1] = self.taps[-1, -1] * ywins[-1, -1] - elif iwin0 == self.dims[0] - 1: - ywins[-1, iwin1] = self.taps[-1, 1] * ywins[-1, iwin1] - elif iwin1 == 0: - ywins[iwin0, 0] = self.taps[1, 0] * ywins[iwin0, 0] - elif iwin1 == self.dims[1] - 1: - ywins[iwin0, -1] = self.taps[1, -1] * ywins[iwin0, -1] - else: - ywins[iwin0, iwin1] = self.taps[1, 1] * ywins[iwin0, iwin1] + ywins = self._apply_taper(ywins, iwin0, iwin1) y = self.Op.H @ ywins else: y = ncp.zeros(self.dims, dtype=self.dtype) for iwin0 in range(self.dims[0]): for iwin1 in range(self.dims[1]): if self.tapertype is not None: - if iwin0 == 0 and iwin1 == 0: - ywins[0, 0] = self.taps[0, 0] * ywins[0, 0] - elif iwin0 == 0 and iwin1 == self.dims[1] - 1: - ywins[0, -1] = self.taps[0, -1] * ywins[0, -1] - elif iwin0 == 0: - ywins[0, iwin1] = self.taps[0, 1] * ywins[0, iwin1] - elif iwin0 == self.dims[0] - 1 and iwin1 == 0: - ywins[-1, 0] = self.taps[-1, 0] * ywins[-1, 0] - elif iwin0 == self.dims[0] - 1 and iwin1 == self.dims[1] - 1: - ywins[-1, -1] = self.taps[-1, -1] * ywins[-1, -1] - elif iwin0 == self.dims[0] - 1: - ywins[-1, iwin1] = self.taps[-1, 1] * ywins[-1, iwin1] - elif iwin1 == 0: - ywins[iwin0, 0] = self.taps[1, 0] * ywins[iwin0, 0] - elif iwin1 == self.dims[1] - 1: - ywins[iwin0, -1] = self.taps[1, -1] * ywins[iwin0, -1] - else: - ywins[iwin0, iwin1] = self.taps[1, 1] * ywins[iwin0, iwin1] + ywins = self._apply_taper(ywins, iwin0, iwin1) y[iwin0, iwin1] = self.Op.rmatvec( ywins[iwin0, iwin1].ravel() ).reshape(self.dims[2], self.dims[3], self.dims[4]) From af5d8ae328d1e419ce41bbfcc43a5fbebe426053 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sat, 18 May 2024 21:57:57 +0300 Subject: [PATCH 10/11] minor: remove unused i --- pylops/signalprocessing/sliding1d.py | 2 +- pylops/signalprocessing/sliding2d.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/pylops/signalprocessing/sliding1d.py b/pylops/signalprocessing/sliding1d.py index 9bb331bb..9ddc4842 100644 --- a/pylops/signalprocessing/sliding1d.py +++ b/pylops/signalprocessing/sliding1d.py @@ -189,7 +189,7 @@ def __init__( self.taps = [ tapin, ] - for i in range(1, nwins - 1): + for _ in range(1, nwins - 1): self.taps.append(tap) self.taps.append(tapend) self.taps = np.vstack(self.taps) diff --git a/pylops/signalprocessing/sliding2d.py b/pylops/signalprocessing/sliding2d.py index 4a76f5e6..ad0be2f5 100644 --- a/pylops/signalprocessing/sliding2d.py +++ b/pylops/signalprocessing/sliding2d.py @@ -230,7 +230,7 @@ def __init__( self.taps = [ tapin[np.newaxis, :], ] - for i in range(1, nwins - 1): + for _ in range(1, nwins - 1): self.taps.append(tap[np.newaxis, :]) self.taps.append(tapend[np.newaxis, :]) self.taps = np.concatenate(self.taps, axis=0) From 4c82d2166c7e5b3bda093d2027f767bb49349c51 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Wed, 22 May 2024 21:00:40 +0300 Subject: [PATCH 11/11] minor: added verb to *_design methods and few other small fixes --- pylops/signalprocessing/patch2d.py | 44 ++++++++++++----------- pylops/signalprocessing/patch3d.py | 54 ++++++++++++++-------------- pylops/signalprocessing/sliding1d.py | 27 ++++++++------ pylops/signalprocessing/sliding2d.py | 27 ++++++++------ pylops/signalprocessing/sliding3d.py | 35 ++++++++++-------- 5 files changed, 103 insertions(+), 84 deletions(-) diff --git a/pylops/signalprocessing/patch2d.py b/pylops/signalprocessing/patch2d.py index 3aa157a5..11b93fca 100644 --- a/pylops/signalprocessing/patch2d.py +++ b/pylops/signalprocessing/patch2d.py @@ -28,6 +28,7 @@ def patch2d_design( nwin: Tuple[int, int], nover: Tuple[int, int], nop: Tuple[int, int], + verb: bool = True, ) -> Tuple[ Tuple[int, int], Tuple[int, int], @@ -51,6 +52,9 @@ def patch2d_design( Number of samples of overlapping part of window. nop : :obj:`tuple` Size of model in the transformed domain. + verb : :obj:`bool`, optional + Verbosity flag. If ``verb==True``, print the data + and model windows start-end indices Returns ------- @@ -79,21 +83,22 @@ def patch2d_design( mwins_inends = ((mwin0_ins, mwin0_ends), (mwin1_ins, mwin1_ends)) # print information about patching - logging.warning("%d-%d windows required...", nwins0, nwins1) - logging.warning( - "data wins - start:%s, end:%s / start:%s, end:%s", - dwin0_ins, - dwin0_ends, - dwin1_ins, - dwin1_ends, - ) - logging.warning( - "model wins - start:%s, end:%s / start:%s, end:%s", - mwin0_ins, - mwin0_ends, - mwin1_ins, - mwin1_ends, - ) + if verb: + logging.warning("%d-%d windows required...", nwins0, nwins1) + logging.warning( + "data wins - start:%s, end:%s / start:%s, end:%s", + dwin0_ins, + dwin0_ends, + dwin1_ins, + dwin1_ends, + ) + logging.warning( + "model wins - start:%s, end:%s / start:%s, end:%s", + mwin0_ins, + mwin0_ends, + mwin1_ins, + mwin1_ends, + ) return nwins, dims, mwins_inends, dwins_inends @@ -276,10 +281,7 @@ def __init__( self.taps = np.vstack(taps).reshape(3, 3, nwin[0], nwin[1]) # define scalings - if scalings is None: - self.scalings = [1.0] * nwins - else: - self.scalings = scalings + self.scalings = [1.0] * nwins if scalings is None else scalings # check if operator is applied to all windows simultaneously self.simOp = False @@ -318,7 +320,7 @@ def _apply_taper(self, ywins, iwin0, iwin1): ywins[iwin0, iwin1] = self.taps[1, 1] * ywins[iwin0, iwin1] return ywins - @reshaped() + @reshaped def _matvec_savetaper(self, x: NDArray) -> NDArray: ncp = get_array_module(x) if self.tapertype is not None: @@ -365,7 +367,7 @@ def _rmatvec_savetaper(self, x: NDArray) -> NDArray: ).reshape(self.dims[2], self.dims[3]) return y - @reshaped() + @reshaped def _matvec_nosavetaper(self, x: NDArray) -> NDArray: ncp = get_array_module(x) if self.tapertype is not None: diff --git a/pylops/signalprocessing/patch3d.py b/pylops/signalprocessing/patch3d.py index 71ffe341..26cd9086 100644 --- a/pylops/signalprocessing/patch3d.py +++ b/pylops/signalprocessing/patch3d.py @@ -28,6 +28,7 @@ def patch3d_design( nwin: Tuple[int, int, int], nover: Tuple[int, int, int], nop: Tuple[int, int, int], + verb: bool = True, ) -> Tuple[ Tuple[int, int, int], Tuple[int, int, int], @@ -51,6 +52,9 @@ def patch3d_design( Number of samples of overlapping part of window. nop : :obj:`tuple` Size of model in the transformed domain. + verb : :obj:`bool`, optional + Verbosity flag. If ``verb==True``, print the data + and model windows start-end indices Returns ------- @@ -90,25 +94,26 @@ def patch3d_design( ) # print information about patching - logging.warning("%d-%d-%d windows required...", nwins0, nwins1, nwins2) - logging.warning( - "data wins - start:%s, end:%s / start:%s, end:%s / start:%s, end:%s", - dwin0_ins, - dwin0_ends, - dwin1_ins, - dwin1_ends, - dwin2_ins, - dwin2_ends, - ) - logging.warning( - "model wins - start:%s, end:%s / start:%s, end:%s / start:%s, end:%s", - mwin0_ins, - mwin0_ends, - mwin1_ins, - mwin1_ends, - mwin2_ins, - mwin2_ends, - ) + if verb: + logging.warning("%d-%d-%d windows required...", nwins0, nwins1, nwins2) + logging.warning( + "data wins - start:%s, end:%s / start:%s, end:%s / start:%s, end:%s", + dwin0_ins, + dwin0_ends, + dwin1_ins, + dwin1_ends, + dwin2_ins, + dwin2_ends, + ) + logging.warning( + "model wins - start:%s, end:%s / start:%s, end:%s / start:%s, end:%s", + mwin0_ins, + mwin0_ends, + mwin1_ins, + mwin1_ends, + mwin2_ins, + mwin2_ends, + ) return nwins, dims, mwins_inends, dwins_inends @@ -485,10 +490,7 @@ def __init__( self.taps[-1, -1, -1] = taprightbottomback # define scalings - if scalings is None: - self.scalings = [1.0] * nwins - else: - self.scalings = scalings + self.scalings = [1.0] * nwins if scalings is None else scalings # check if operator is applied to all windows simultaneously self.simOp = False @@ -576,7 +578,7 @@ def _apply_taper(self, ywins, iwin0, iwin1, iwin2): ywins[iwin0, iwin1, iwin2] = self.taps[1, 1, 1] * ywins[iwin0, iwin1, iwin2] return ywins - @reshaped() + @reshaped def _matvec_savetaper(self, x: NDArray) -> NDArray: ncp = get_array_module(x) if self.tapertype is not None: @@ -630,7 +632,7 @@ def _rmatvec_savetaper(self, x: NDArray) -> NDArray: ).reshape(self.dims[3], self.dims[4], self.dims[5]) return y - @reshaped() + @reshaped def _matvec_nosavetaper(self, x: NDArray) -> NDArray: ncp = get_array_module(x) if self.tapertype is not None: @@ -727,7 +729,7 @@ def _matvec_nosavetaper(self, x: NDArray) -> NDArray: ] += xxwin return y - @reshaped() + @reshaped def _rmatvec_nosavetaper(self, x: NDArray) -> NDArray: ncp = get_array_module(x) ncp_sliding_window_view = get_sliding_window_view(x) diff --git a/pylops/signalprocessing/sliding1d.py b/pylops/signalprocessing/sliding1d.py index 9ddc4842..b203c1a1 100644 --- a/pylops/signalprocessing/sliding1d.py +++ b/pylops/signalprocessing/sliding1d.py @@ -28,6 +28,7 @@ def sliding1d_design( nwin: int, nover: int, nop: int, + verb: bool = True, ) -> Tuple[int, int, Tuple[NDArray, NDArray], Tuple[NDArray, NDArray]]: """Design Sliding1D operator @@ -46,6 +47,9 @@ def sliding1d_design( Number of samples of overlapping part of window. nop : :obj:`tuple` Size of model in the transformed domain. + verb : :obj:`bool`, optional + Verbosity flag. If ``verb==True``, print the data + and model windows start-end indices Returns ------- @@ -70,17 +74,18 @@ def sliding1d_design( mwins_inends = (mwin_ins, mwin_ends) # print information about patching - logging.warning("%d windows required...", nwins) - logging.warning( - "data wins - start:%s, end:%s", - dwin_ins, - dwin_ends, - ) - logging.warning( - "model wins - start:%s, end:%s", - mwin_ins, - mwin_ends, - ) + if verb: + logging.warning("%d windows required...", nwins) + logging.warning( + "data wins - start:%s, end:%s", + dwin_ins, + dwin_ends, + ) + logging.warning( + "model wins - start:%s, end:%s", + mwin_ins, + mwin_ends, + ) return nwins, dim, mwins_inends, dwins_inends diff --git a/pylops/signalprocessing/sliding2d.py b/pylops/signalprocessing/sliding2d.py index ad0be2f5..27c40089 100644 --- a/pylops/signalprocessing/sliding2d.py +++ b/pylops/signalprocessing/sliding2d.py @@ -60,6 +60,7 @@ def sliding2d_design( nwin: int, nover: int, nop: Tuple[int, int], + verb: bool = True, ) -> Tuple[int, Tuple[int, int], Tuple[NDArray, NDArray], Tuple[NDArray, NDArray]]: """Design Sliding2D operator @@ -78,6 +79,9 @@ def sliding2d_design( Number of samples of overlapping part of window. nop : :obj:`tuple` Size of model in the transformed domain. + verb : :obj:`bool`, optional + Verbosity flag. If ``verb==True``, print the data + and model windows start-end indices Returns ------- @@ -102,17 +106,18 @@ def sliding2d_design( mwins_inends = (mwin_ins, mwin_ends) # print information about patching - logging.warning("%d windows required...", nwins) - logging.warning( - "data wins - start:%s, end:%s", - dwin_ins, - dwin_ends, - ) - logging.warning( - "model wins - start:%s, end:%s", - mwin_ins, - mwin_ends, - ) + if verb: + logging.warning("%d windows required...", nwins) + logging.warning( + "data wins - start:%s, end:%s", + dwin_ins, + dwin_ends, + ) + logging.warning( + "model wins - start:%s, end:%s", + mwin_ins, + mwin_ends, + ) return nwins, dims, mwins_inends, dwins_inends diff --git a/pylops/signalprocessing/sliding3d.py b/pylops/signalprocessing/sliding3d.py index 49324230..db052180 100644 --- a/pylops/signalprocessing/sliding3d.py +++ b/pylops/signalprocessing/sliding3d.py @@ -28,6 +28,7 @@ def sliding3d_design( nwin: Tuple[int, int], nover: Tuple[int, int], nop: Tuple[int, int, int], + verb: bool = True, ) -> Tuple[ Tuple[int, int], Tuple[int, int, int], @@ -51,6 +52,9 @@ def sliding3d_design( Number of samples of overlapping part of window. nop : :obj:`tuple` Size of model in the transformed domain. + verb : :obj:`bool`, optional + Verbosity flag. If ``verb==True``, print the data + and model windows start-end indices Returns ------- @@ -79,21 +83,22 @@ def sliding3d_design( mwins_inends = ((mwin0_ins, mwin0_ends), (mwin1_ins, mwin1_ends)) # print information about patching - logging.warning("%d-%d windows required...", nwins0, nwins1) - logging.warning( - "data wins - start:%s, end:%s / start:%s, end:%s", - dwin0_ins, - dwin0_ends, - dwin1_ins, - dwin1_ends, - ) - logging.warning( - "model wins - start:%s, end:%s / start:%s, end:%s", - mwin0_ins, - mwin0_ends, - mwin1_ins, - mwin1_ends, - ) + if verb: + logging.warning("%d-%d windows required...", nwins0, nwins1) + logging.warning( + "data wins - start:%s, end:%s / start:%s, end:%s", + dwin0_ins, + dwin0_ends, + dwin1_ins, + dwin1_ends, + ) + logging.warning( + "model wins - start:%s, end:%s / start:%s, end:%s", + mwin0_ins, + mwin0_ends, + mwin1_ins, + mwin1_ends, + ) return nwins, dims, mwins_inends, dwins_inends