From e1e367509812f065cf8a4b8bdeeae93cb0d37187 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sun, 26 Mar 2023 21:57:22 +0300 Subject: [PATCH 01/69] fix: modify solvers to use matvec/rmatvec instead of @/.H @ This commits modify all solvers to use matvec/rmatvec as they are more performant than their equivalent @/.H @. Note that only in the case of ISTA and FISTA the @/.H @ are still used to allow working with either single or multiple right-hand-sides. --- MIGRATION_V1_V2.md | 2 +- pylops/optimization/cls_leastsquares.py | 36 ++++++-------- pylops/optimization/cls_sparsity.py | 62 ++++++++++++------------- 3 files changed, 44 insertions(+), 56 deletions(-) diff --git a/MIGRATION_V1_V2.md b/MIGRATION_V1_V2.md index e198965d..d26271b8 100644 --- a/MIGRATION_V1_V2.md +++ b/MIGRATION_V1_V2.md @@ -10,7 +10,7 @@ should be used as a checklist when converting a piece of code using PyLops from for every operator by default. While the change is mostly backwards compatible, there are some operators (e.g. the ``Bilinear`` transpose/conjugate) which can output reshaped arrays instead of 1d-arrays. To ensure no breakage, you can entirely disable this feature either globally by setting ``pylops.set_ndarray_multiplication(False)``, or locally with the context manager - ``pylops.disabled_ndarray_multiplication()``. Both will revert to v1.x behavior. At this time, PyLops sparse solvers do + ``pylops.disabled_ndarray_multiplication()``. Both will revert to v1.x behavior. At this time, PyLops solvers do *not* support N-D array multiplication. See the table at the end of this document for support ndarray operations. diff --git a/pylops/optimization/cls_leastsquares.py b/pylops/optimization/cls_leastsquares.py index 8382a2fc..a9e20fec 100644 --- a/pylops/optimization/cls_leastsquares.py +++ b/pylops/optimization/cls_leastsquares.py @@ -16,7 +16,6 @@ from pylops.optimization.basesolver import Solver from pylops.optimization.basic import cg, cgls from pylops.utils.backend import get_array_module -from pylops.utils.decorators import disable_ndarray_multiplication from pylops.utils.typing import NDArray if TYPE_CHECKING: @@ -89,7 +88,6 @@ def _print_finalize(self) -> None: print(f"\nTotal time (s) = {self.telapsed:.2f}") print("-" * 55 + "\n") - @disable_ndarray_multiplication def setup( self, y: NDArray, @@ -157,13 +155,13 @@ def setup( # normal equations if Weight is not None: - self.y_normal = self.OpH * Weight * y + self.y_normal = self.Op.rmatvec(Weight.matvec(y)) else: - self.y_normal = self.OpH * y + self.y_normal = self.Op.rmatvec(y) if Weight is not None: - self.Op_normal = self.OpH * Weight * self.Op + self.Op_normal = self.OpH @ Weight @ self.Op else: - self.Op_normal = self.OpH * self.Op + self.Op_normal = self.OpH @ self.Op # add regularization terms if epsI > 0: @@ -178,9 +176,8 @@ def setup( and self.dataregs is not None ): for epsR, Reg, datareg in zip(self.epsRs, self.Regs, self.dataregs): - self.RegH = Reg.H - self.y_normal += epsR**2 * self.RegH * datareg - self.Op_normal += epsR**2 * self.RegH * Reg + self.y_normal += epsR**2 * Reg.rmatvec(datareg) + self.Op_normal += epsR**2 * Reg.H @ Reg if epsNRs is not None and NRegs is not None: for epsNR, NReg in zip(epsNRs, NRegs): @@ -197,7 +194,6 @@ def step(self) -> None: "step method is not implemented. Use directly run or solve." ) - @disable_ndarray_multiplication def run( self, x: NDArray, @@ -240,7 +236,7 @@ def run( """ if x is not None: - self.y_normal = self.y_normal - self.Op_normal * x + self.y_normal = self.y_normal - self.Op_normal.matvec(x) if engine == "scipy" and self.ncp == np: if "atol" not in kwargs_solver: kwargs_solver["atol"] = "legacy" @@ -451,7 +447,6 @@ def _print_finalize(self) -> None: print(f"\nTotal time (s) = {self.telapsed:.2f}") print("-" * 65 + "\n") - @disable_ndarray_multiplication def setup( self, y: NDArray, @@ -504,10 +499,10 @@ def setup( self.RegOp: LinearOperator if Weight is not None: if Regs is None: - self.RegOp = Weight * self.Op + self.RegOp = Weight @ self.Op else: self.RegOp = RegularizedOperator( - Weight * self.Op, Regs, epsRs=self.epsRs + Weight @ self.Op, Regs, epsRs=self.epsRs ) else: if Regs is None: @@ -517,7 +512,7 @@ def setup( # augumented data if Weight is not None: - self.datatot: NDArray = Weight * self.y.copy() + self.datatot: NDArray = Weight.matvec(self.y.copy()) else: self.datatot = self.y.copy() @@ -537,7 +532,6 @@ def step(self) -> None: "step method is not implemented. Use directly run or solve." ) - @disable_ndarray_multiplication def run( self, x: NDArray, @@ -586,7 +580,7 @@ def run( """ if x is not None: - self.datatot = self.datatot - self.RegOp * x + self.datatot = self.datatot - self.RegOp.matvec(x) if engine == "scipy" and self.ncp == np: if show: kwargs_solver["show"] = 1 @@ -722,7 +716,6 @@ def _print_finalize(self) -> None: print(f"\nTotal time (s) = {self.telapsed:.2f}") print("-" * 65 + "\n") - @disable_ndarray_multiplication def setup( self, y: NDArray, @@ -746,7 +739,7 @@ def setup( self.ncp = get_array_module(y) # preconditioned operator - self.POp = self.Op * P + self.POp = self.Op @ P # print setup if show: @@ -759,7 +752,6 @@ def step(self) -> None: "step method is not implemented. Use directly run or solve." ) - @disable_ndarray_multiplication def run( self, x: NDArray, @@ -808,7 +800,7 @@ def run( """ if x is not None: - self.y = self.y - self.Op * x + self.y = self.y - self.Op.matvec(x) if engine == "scipy" and self.ncp == np: if show: kwargs_solver["show"] = 1 @@ -830,7 +822,7 @@ def run( pinv = pinv.ravel() else: raise NotImplementedError("Engine must be scipy or pylops") - xinv = self.P * pinv + xinv = self.P.matvec(pinv) if x is not None: xinv = x + xinv return xinv, istop, itn, r1norm, r2norm diff --git a/pylops/optimization/cls_sparsity.py b/pylops/optimization/cls_sparsity.py index 11f38aee..80aae105 100644 --- a/pylops/optimization/cls_sparsity.py +++ b/pylops/optimization/cls_sparsity.py @@ -14,7 +14,6 @@ from pylops.optimization.leastsquares import regularized_inversion from pylops.utils import deps from pylops.utils.backend import get_array_module, get_module_name -from pylops.utils.decorators import disable_ndarray_multiplication from pylops.utils.typing import InputDimsLike, NDArray, SamplingLike spgl1_message = deps.spgl1_import("the spgl1 solver") @@ -422,18 +421,16 @@ def _step_model(self, x: NDArray, **kwargs_solver) -> NDArray: if self.iiter == 0: # first iteration (unweighted least-squares) if self.ncp == np: - x = ( - self.Op.H - @ lsqr( + x = self.Op.rmatvec( + lsqr( self.Op @ self.Op.H + (self.epsI**2) * self.Iop, self.y, **kwargs_solver, )[0] ) else: - x = ( - self.Op.H - @ cgls( + x = self.Op.rmatvec( + cgls( self.Op @ self.Op.H + (self.epsI**2) * self.Iop, self.y, self.ncp.zeros(int(self.Op.shape[0]), dtype=self.Op.dtype), @@ -446,25 +443,25 @@ def _step_model(self, x: NDArray, **kwargs_solver) -> NDArray: self.rw = self.rw / self.rw.max() R = Diagonal(self.rw, dtype=self.rw.dtype) if self.ncp == np: - x = ( - R - @ self.Op.H - @ lsqr( - self.Op @ R @ self.Op.H + self.epsI**2 * self.Iop, - self.y, - **kwargs_solver, - )[0] + x = R.matvec( + self.Op.rmatvec( + lsqr( + self.Op @ R @ self.Op.H + self.epsI**2 * self.Iop, + self.y, + **kwargs_solver, + )[0] + ) ) else: - x = ( - R - @ self.Op.H - @ cgls( - self.Op @ R @ self.Op.H + self.epsI**2 * self.Iop, - self.y, - self.ncp.zeros(int(self.Op.shape[0]), dtype=self.Op.dtype), - **kwargs_solver, - )[0] + x = R.matvec( + self.Op.rmatvec( + cgls( + self.Op @ R @ self.Op.H + self.epsI**2 * self.Iop, + self.y, + self.ncp.zeros(int(self.Op.shape[0]), dtype=self.Op.dtype), + **kwargs_solver, + )[0] + ) ) return x @@ -1140,7 +1137,8 @@ def setup( Parameters ---------- y : :obj:`np.ndarray` - Data of size :math:`[N \times 1]` + Data of size :math:`[N \times 1]` or :math:`[N \times R]` where for a solution for multiple right-hand-side + is found when ``R>1``. x0: :obj:`numpy.ndarray`, optional Initial guess niter : :obj:`int` @@ -1330,10 +1328,10 @@ def step(self, x: NDArray, show: bool = False) -> Tuple[NDArray, float]: else: self.normresold = self.normres - # compute gradient + # compute gradient (must use .H @ for the solver to work with either single or multiple rhs) grad: NDArray = self.alpha * (self.Op.H @ res) - # update inverted model + # update inverted model (must use @ / .H @ for the solver to work with either single or multiple rhs) x_unthesh: NDArray = x + grad if self.SOp is not None: x_unthesh = self.SOp.H @ x_unthesh @@ -1578,7 +1576,7 @@ def step(self, x: NDArray, z: NDArray, show: bool = False) -> NDArray: ---------- x : :obj:`np.ndarray` Current model vector to be updated by a step of ISTA - x : :obj:`np.ndarray` + z : :obj:`np.ndarray` Current auxiliary model vector to be updated by a step of ISTA show : :obj:`bool`, optional Display iteration log @@ -1608,10 +1606,10 @@ def step(self, x: NDArray, z: NDArray, show: bool = False) -> NDArray: else: self.normresold = self.normres - # compute gradient + # compute gradient (must use .H @ for the solver to work with either single or multiple rhs grad: NDArray = self.alpha * (self.Op.H @ resz) - # update inverted model + # update inverted model (must use @ / .H @ for the solver to work with either single or multiple rhs) x_unthesh: NDArray = z + grad if self.SOp is not None: x_unthesh = self.SOp.H @ x_unthesh @@ -1753,7 +1751,6 @@ def _print_finalize(self) -> None: print(f"\nTotal time (s) = {self.telapsed:.2f}") print("-" * 80 + "\n") - @disable_ndarray_multiplication def setup( self, y: NDArray, @@ -1800,7 +1797,6 @@ def step(self) -> None: "step method is not implemented. Use directly run or solve." ) - @disable_ndarray_multiplication def run( self, x: NDArray, @@ -1879,7 +1875,7 @@ def run( **kwargs_spgl1, ) - xinv = pinv.copy() if self.SOp is None else self.SOp.H * pinv + xinv = pinv.copy() if self.SOp is None else self.SOp.rmatvec(pinv) return xinv, pinv, info def solve( From 777eac2e4663f39278870ccd9afbbbf366365a9d Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sun, 26 Mar 2023 21:59:05 +0300 Subject: [PATCH 02/69] feature: added forceflat to operators with ambiguous rmatvec Some operators produce 1D arrays in the forward pass (matvec). This makes it impossible to know whether a 1D array or a ND array should be returned in the adjoint pass (rmatvec). By introducing the optional input parameter forceflat a user can specify the preferred behaviour (default: False - rmatvec returns ND array). The class LinearOperator is also modified such that forceflat drives the dot and __add__ methods. --- pylops/basicoperators/sum.py | 13 ++++- pylops/linearoperator.py | 77 ++++++++++++++++++++++++++--- pylops/signalprocessing/bilinear.py | 13 ++++- tutorials/linearoperator.py | 55 +++++++++++++++------ 4 files changed, 134 insertions(+), 24 deletions(-) diff --git a/pylops/basicoperators/sum.py b/pylops/basicoperators/sum.py index ec42cabe..751bbba1 100644 --- a/pylops/basicoperators/sum.py +++ b/pylops/basicoperators/sum.py @@ -24,6 +24,10 @@ class Sum(LinearOperator): .. versionadded:: 2.0.0 Axis along which model is summed. + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after rmatvec. dtype : :obj:`str`, optional Type of elements in input array. name : :obj:`str`, optional @@ -61,6 +65,7 @@ def __init__( self, dims: InputDimsLike, axis: int = -1, + forceflat: bool = False, dtype: DTypeLike = "float64", name: str = "S", ) -> None: @@ -71,7 +76,13 @@ def __init__( # data dimensions dimsd = list(dims).copy() dimsd.pop(self.axis) - super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dimsd, name=name) + super().__init__( + dtype=np.dtype(dtype), + dims=dims, + dimsd=dimsd, + forceflat=forceflat, + name=name, + ) # array of ones with dims of model in self.axis for np.tile in adjoint mode self.tile = np.ones(len(self.dims), dtype=int) diff --git a/pylops/linearoperator.py b/pylops/linearoperator.py index b16358cb..1997dd03 100644 --- a/pylops/linearoperator.py +++ b/pylops/linearoperator.py @@ -85,13 +85,19 @@ class LinearOperator(_LinearOperator): .. versionadded:: 2.0.0 Dimensions of data. If not provided, ``(self.shape[0],)`` is used. - explicit : :obj:`bool`, optional - Operator contains a matrix that can be solved explicitly - (``True``) or not (``False``). Defaults to ``False``. clinear : :obj:`bool`, optional .. versionadded:: 1.17.0 Operator is complex-linear. Defaults to ``True``. + explicit : :obj:`bool`, optional + Operator contains a matrix that can be solved explicitly + (``True``) or not (``False``). Defaults to ``False``. + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after matvec/rmatvec if the input is ambiguous (i.e., is a 1D array both when + operating with ND arrays and with 1D arrays. Defaults to ``None`` for operators that have no ambiguity (and + to ``False`` for those with ambiguity) name : :obj:`str`, optional .. versionadded:: 2.0.0 @@ -108,6 +114,7 @@ def __init__( dimsd: Optional[ShapeLike] = None, clinear: Optional[bool] = None, explicit: Optional[bool] = None, + forceflat: Optional[bool] = None, name: Optional[str] = None, ) -> None: if Op is not None: @@ -124,6 +131,9 @@ def __init__( ) if explicit and hasattr(Op, "A"): self.A = Op.A + forceflat = ( + getattr(self.Op, "forceflat", None) if forceflat is None else forceflat + ) name = getattr(Op, "name", None) if name is None else name if dtype is not None: @@ -138,6 +148,8 @@ def __init__( self.clinear = clinear if explicit is not None: self.explicit = explicit + if forceflat is not None: + self.forceflat = forceflat self.name = name # counters @@ -263,6 +275,18 @@ def explicit(self, new_explicit: bool) -> None: def explicit(self): del self._explicit + @property + def forceflat(self): + return getattr(self, "_forceflat", None) + + @forceflat.setter + def forceflat(self, new_forceflat: bool) -> None: + self._forceflat = bool(new_forceflat) + + @forceflat.deleter + def forceflat(self): + del self._forceflat + @property def name(self): return getattr(self, "_name", None) @@ -328,11 +352,26 @@ def __add__(self, x: LinearOperator) -> LinearOperator: Op, exclude=[ "explicit", + "forceflat", "name", ], ) Op.clinear = Op.clinear and Opx.clinear Op.explicit = False + # Define forceflat (if differing, raise error) + if isinstance(self.forceflat, bool) and isinstance(Opx.forceflat, bool): + if self.forceflat != Opx.forceflat: + raise ValueError( + f"the two operators have contrasting forceflat {Op.forceflat}-{Opx.forceflatx}" + ) + else: + Op.forceflat = self.forceflat + if isinstance(self.forceflat, bool) or isinstance(Opx.forceflat, bool): + Op.forceflat = ( + self.forceflat if self.forceflat is not None else Opx.forceflat + ) + else: + Op.forceflat = None # Replace if shape-like if len(self.dims) == 1: Op.dims = Opx.dims @@ -374,7 +413,7 @@ def _copy_attributes( """Copy attributes from one LinearOperator to another""" if exclude is None: exclude = ["name"] - attrs = ["dims", "dimsd", "clinear", "explicit", "name"] + attrs = ["dims", "dimsd", "clinear", "explicit", "forceflat", "name"] if exclude is not None: for item in exclude: attrs.remove(item) @@ -582,9 +621,23 @@ def dot(self, x: NDArray) -> NDArray: # to allow mixing pylops and scipy operators) Opx = aslinearoperator(x) Op = _ProductLinearOperator(self, Opx) - self._copy_attributes(Op, exclude=["dims", "explicit", "name"]) + self._copy_attributes(Op, exclude=["dims", "explicit", "forceflat", "name"]) Op.clinear = Op.clinear and Opx.clinear Op.explicit = False + # Define forceflat (if differing, raise error) + if isinstance(self.forceflat, bool) and isinstance(Opx.forceflat, bool): + if self.forceflat != Opx.forceflat: + raise ValueError( + f"the two operators have contrasting forceflat {Op.forceflat}-{Opx.forceflatx}" + ) + else: + Op.forceflat = self.forceflat + if isinstance(self.forceflat, bool) or isinstance(Opx.forceflat, bool): + Op.forceflat = ( + self.forceflat if self.forceflat is not None else Opx.forceflat + ) + else: + Op.forceflat = None Op.dims = Opx.dims return Op elif np.isscalar(x): @@ -608,17 +661,25 @@ def dot(self, x: NDArray) -> NDArray: if is_dims_shaped: # (dims1, ..., dimsK) => (dims1 * ... * dimsK,) == self.shape x = x.ravel() - if is_dims_shaped_matrix: + if is_dims_shaped_matrix and not self.forceflat: # (dims1, ..., dimsK, P) => (dims1 * ... * dimsK, P) x = x.reshape((-1, x.shape[-1])) if x.ndim == 1: y = self.matvec(x) - if is_dims_shaped and get_ndarray_multiplication(): + if ( + is_dims_shaped + and not self.forceflat + and get_ndarray_multiplication() + ): y = y.reshape(self.dimsd) return y elif x.ndim == 2: y = self.matmat(x) - if is_dims_shaped_matrix and get_ndarray_multiplication(): + if ( + is_dims_shaped_matrix + and not self.forceflat + and get_ndarray_multiplication() + ): y = y.reshape((*self.dimsd, -1)) return y else: diff --git a/pylops/signalprocessing/bilinear.py b/pylops/signalprocessing/bilinear.py index 7924e602..395683d9 100644 --- a/pylops/signalprocessing/bilinear.py +++ b/pylops/signalprocessing/bilinear.py @@ -37,6 +37,10 @@ class Bilinear(LinearOperator): for interpolation. dims : :obj:`list` Number of samples for each dimension + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after rmatvec. dtype : :obj:`str`, optional Type of elements in input array. name : :obj:`str`, optional @@ -90,13 +94,20 @@ def __init__( self, iava: IntNDArray, dims: InputDimsLike, + forceflat: bool = False, dtype: DTypeLike = "float64", name: str = "B", ) -> None: # define dimension of data ndims = len(dims) dimsd = [len(iava[1])] + list(dims[2:]) - super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dimsd, name=name) + super().__init__( + dtype=np.dtype(dtype), + dims=dims, + dimsd=dimsd, + forceflat=forceflat, + name=name, + ) ncp = get_array_module(iava) # check non-unique pairs (works only with numpy arrays) diff --git a/tutorials/linearoperator.py b/tutorials/linearoperator.py index e1b86669..6c85f9c4 100755 --- a/tutorials/linearoperator.py +++ b/tutorials/linearoperator.py @@ -5,11 +5,17 @@ library for both new users and developers. We will start by looking at how to initialize a linear operator as well as -different ways to apply the forward and adjoint operations. Finally we will +different ways to apply the forward and adjoint operations. We will then investigate various *special methods*, also called *magic methods* (i.e., methods with the double underscores at the beginning and the end) that -have been implemented for such a class and will allow summing, subtractring, +have been implemented for such a class and will allow summing, subtracting, chaining, etc. multiple operators in very easy and expressive way. + +Finally, we will consider a scenario where the input and output vectors of an +operator are naturally represented by nd-arrays. Since pylops V2, we can both +pass to the operator a flattened array (original behaviour) as well as the +nd-array with shape `dims` (new behaviour); according to input, the output +will be either a vector or a nd-array with shape `dimsd`. """ ############################################################################### @@ -112,10 +118,10 @@ cmd2 = "Dop.rmatvec(x)" # .H* (pre-computed H) -cmd3 = "DopH*x" +cmd3 = "DopH@x" # .H* -cmd4 = "Dop.H*x" +cmd4 = "Dop.H@x" # timing t1 = 1.0e3 * np.array(timeit.repeat(cmd1, setup=cmd_setup, number=500, repeat=5)) @@ -171,7 +177,7 @@ print(Dop**3) # * and / -y = Dop * x +y = Dop @ x print(Dop / y) # eigs @@ -194,14 +200,14 @@ x = np.ones(n) Dop = pylops.Diagonal(d) -print(f"y = Dx = {Dop * x}") -print(f"y = conj(D)x = {Dop.conj() * x}") +print(f"y = Dx = {Dop @ x}") +print(f"y = conj(D)x = {Dop.conj() @ x}") ############################################################################### -# At this point, the concept of linear operator may sound abstract. -# The convinience method :func:`pylops.LinearOperator.todense` can be used to -# create the equivalent dense matrix of any operator. In this case for example -# we expect to see a diagonal matrix with ``d`` values along the main diagonal +# The concept of linear operator may sound abstract. The convenience method +# :func:`pylops.LinearOperator.todense` can be used to create the equivalent +# dense matrix of any operator. In this case for example we expect to see a +# diagonal matrix with ``d`` values along the main diagonal D = Dop.todense() plt.figure(figsize=(5, 5)) @@ -212,14 +218,14 @@ plt.tight_layout() ############################################################################### -# At this point it is worth reiterating that if two linear operators are +# At this point, it is worth reiterating that if two linear operators are # combined by means of the algebraical operations shown above, the resulting # operator is still a :py:class:`pylops.LinearOperator` operator. This means # that we can still apply any of the methods implemented in our class # like ``*`` or ``/``. Dop1 = Dop - Dop.conj() -y = Dop1 * x +y = Dop1 @ x print(f"x = (Dop - conj(Dop))/y = {Dop1 / y}") D1 = Dop1.todense() @@ -232,7 +238,7 @@ plt.tight_layout() ############################################################################### -# Finally, another important feature of PyLops linear operators is that we can +# Another important feature of PyLops linear operators is that we can # always keep track of how many times the forward and adjoint passes have been # applied (and reset when needed). This is particularly useful when running a # third party solver to see how many evaluations of our operator are performed @@ -252,6 +258,27 @@ print(f"Forward evaluations: {Dop.matvec_count}") print(f"Adjoint evaluations: {Dop.rmatvec_count}") + +######################################################################## +# Finally, let's consider a case where the :py:class:`pylops.basicoperators.Diagonal` +# operator acts on one dimension of a 2d-array. We will see how in this case the forward +# and adjoint passes can be applied on the flattened array as well as on the 2d-array directly + +# +m, n = 10, 5 +d = np.arange(n) + 1.0 +x = np.ones((m, n)) +Dop = pylops.Diagonal(d, dims=(m, n), axis=-1) + +yflat = Dop @ x.ravel() +y = Dop @ x + +xadjflat = Dop.H @ yflat +xadj = Dop.H @ y + +print(f"yflat shape = {yflat.shape}, y shape = {y.shape}") +print(f"xadjflat shape = {xadjflat.shape}, xadj shape = {xadj.shape}") + ############################################################################### # This first tutorial is completed. You have seen the basic operations that # can be performed using :py:class:`pylops.LinearOperator` and you From acba1c56d40f0792d18d93078f92733e91a1cdff Mon Sep 17 00:00:00 2001 From: mrava87 Date: Mon, 27 Mar 2023 10:02:41 +0300 Subject: [PATCH 03/69] minor: add check for forceflat based on len(dims) A check on the user provided value of forceflat is added based on len(dims), such that if the input array has a shape larger than 2, forceflat is defaulted back to None as pylops can handle internally the operator without users having to choose how to return the outputs of rmatvec. --- pylops/basicoperators/sum.py | 20 ++++++++++++++++++-- pylops/signalprocessing/bilinear.py | 16 ++++++++++++++-- 2 files changed, 32 insertions(+), 4 deletions(-) diff --git a/pylops/basicoperators/sum.py b/pylops/basicoperators/sum.py index 751bbba1..c3066514 100644 --- a/pylops/basicoperators/sum.py +++ b/pylops/basicoperators/sum.py @@ -1,5 +1,7 @@ __all__ = ["Sum"] +import logging + import numpy as np from pylops import LinearOperator @@ -8,6 +10,8 @@ from pylops.utils.decorators import reshaped from pylops.utils.typing import DTypeLike, InputDimsLike, NDArray +logging.basicConfig(format="%(levelname)s: %(message)s", level=logging.WARNING) + class Sum(LinearOperator): r"""Sum operator. @@ -27,7 +31,9 @@ class Sum(LinearOperator): forceflat : :obj:`bool`, optional .. versionadded:: 2.2.0 - Force an array to be flattened after rmatvec. + Force an array to be flattened after rmatvec. Note that this is only + required when `len(dims)=2`, otherwise pylops will detect whether to + return a 1d or nd array. dtype : :obj:`str`, optional Type of elements in input array. name : :obj:`str`, optional @@ -65,7 +71,7 @@ def __init__( self, dims: InputDimsLike, axis: int = -1, - forceflat: bool = False, + forceflat: bool = None, dtype: DTypeLike = "float64", name: str = "S", ) -> None: @@ -76,6 +82,16 @@ def __init__( # data dimensions dimsd = list(dims).copy() dimsd.pop(self.axis) + # check if forceflat is needed and set it back to None otherwise + if len(dims) > 2: + if forceflat is not None: + logging.warning( + f"setting forceflat=None since len(dims)={len(dims)}>2. " + f"PyLops will automatically detect whether to return " + f"a 1d or nd array based on the shape of the input" + f"array." + ) + forceflat = None super().__init__( dtype=np.dtype(dtype), dims=dims, diff --git a/pylops/signalprocessing/bilinear.py b/pylops/signalprocessing/bilinear.py index 395683d9..6bc62e58 100644 --- a/pylops/signalprocessing/bilinear.py +++ b/pylops/signalprocessing/bilinear.py @@ -40,7 +40,9 @@ class Bilinear(LinearOperator): forceflat : :obj:`bool`, optional .. versionadded:: 2.2.0 - Force an array to be flattened after rmatvec. + Force an array to be flattened after rmatvec. Note that this is only + required when `len(dims)=2`, otherwise pylops will detect whether to + return a 1d or nd array. dtype : :obj:`str`, optional Type of elements in input array. name : :obj:`str`, optional @@ -94,13 +96,23 @@ def __init__( self, iava: IntNDArray, dims: InputDimsLike, - forceflat: bool = False, + forceflat: bool = None, dtype: DTypeLike = "float64", name: str = "B", ) -> None: # define dimension of data ndims = len(dims) dimsd = [len(iava[1])] + list(dims[2:]) + # check if forceflat is needed and set it back to None otherwise + if ndims > 2: + if forceflat is not None: + logging.warning( + f"setting forceflat=None since len(dims)={len(dims)}>2. " + f"PyLops will automatically detect whether to return " + f"a 1d or nd array based on the shape of the input" + f"array." + ) + forceflat = None super().__init__( dtype=np.dtype(dtype), dims=dims, From e7045d06e2f2df9104d73bc4162191d94a3b7e86 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Mon, 27 Mar 2023 11:36:26 +0300 Subject: [PATCH 04/69] minor: fix typo in LinearOperator --- pylops/linearoperator.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pylops/linearoperator.py b/pylops/linearoperator.py index 1997dd03..c3d793a1 100644 --- a/pylops/linearoperator.py +++ b/pylops/linearoperator.py @@ -362,7 +362,7 @@ def __add__(self, x: LinearOperator) -> LinearOperator: if isinstance(self.forceflat, bool) and isinstance(Opx.forceflat, bool): if self.forceflat != Opx.forceflat: raise ValueError( - f"the two operators have contrasting forceflat {Op.forceflat}-{Opx.forceflatx}" + f"the two operators have contrasting forceflat {Op.forceflat}-{Opx.forceflat}" ) else: Op.forceflat = self.forceflat @@ -628,7 +628,7 @@ def dot(self, x: NDArray) -> NDArray: if isinstance(self.forceflat, bool) and isinstance(Opx.forceflat, bool): if self.forceflat != Opx.forceflat: raise ValueError( - f"the two operators have contrasting forceflat {Op.forceflat}-{Opx.forceflatx}" + f"the two operators have contrasting forceflat {Op.forceflat}-{Opx.forceflat}" ) else: Op.forceflat = self.forceflat From 3e9ead498f212eaaf7d4959282877374bdfc51f6 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Mon, 27 Mar 2023 15:22:35 +0300 Subject: [PATCH 05/69] feature: added forceflat on HStack and VStack --- pylops/basicoperators/hstack.py | 13 ++++++++++++- pylops/basicoperators/vstack.py | 13 ++++++++++++- pylops/linearoperator.py | 5 ++++- 3 files changed, 28 insertions(+), 3 deletions(-) diff --git a/pylops/basicoperators/hstack.py b/pylops/basicoperators/hstack.py index c97fac1c..775ecace 100644 --- a/pylops/basicoperators/hstack.py +++ b/pylops/basicoperators/hstack.py @@ -44,6 +44,10 @@ class HStack(LinearOperator): nproc : :obj:`int`, optional Number of processes used to evaluate the N operators in parallel using ``multiprocessing``. If ``nproc=1``, work in serial mode. + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after matvec. dtype : :obj:`str`, optional Type of elements in input array. @@ -107,6 +111,7 @@ def __init__( self, ops: Sequence[LinearOperator], nproc: int = 1, + forceflat: bool = None, dtype: Optional[str] = None, ) -> None: self.ops = ops @@ -129,7 +134,13 @@ def __init__( dtype = _get_dtype(self.ops) if dtype is None else np.dtype(dtype) clinear = all([getattr(oper, "clinear", True) for oper in self.ops]) - super().__init__(dtype=dtype, shape=(self.nops, self.mops), clinear=clinear) + super().__init__( + dtype=dtype, + shape=(self.nops, self.mops), + dimsd=ops[0].dimsd, + clinear=clinear, + forceflat=forceflat, + ) @property def nproc(self) -> int: diff --git a/pylops/basicoperators/vstack.py b/pylops/basicoperators/vstack.py index dc6734f8..d897e0f5 100644 --- a/pylops/basicoperators/vstack.py +++ b/pylops/basicoperators/vstack.py @@ -44,6 +44,10 @@ class VStack(LinearOperator): nproc : :obj:`int`, optional Number of processes used to evaluate the N operators in parallel using ``multiprocessing``. If ``nproc=1``, work in serial mode. + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after rmatvec. dtype : :obj:`str`, optional Type of elements in input array. @@ -107,6 +111,7 @@ def __init__( self, ops: Sequence[LinearOperator], nproc: int = 1, + forceflat: bool = None, dtype: Optional[DTypeLike] = None, ) -> None: self.ops = ops @@ -129,7 +134,13 @@ def __init__( dtype = _get_dtype(self.ops) if dtype is None else np.dtype(dtype) clinear = all([getattr(oper, "clinear", True) for oper in self.ops]) - super().__init__(dtype=dtype, shape=(self.nops, self.mops), clinear=clinear) + super().__init__( + dtype=dtype, + shape=(self.nops, self.mops), + dims=ops[0].dims, + clinear=clinear, + forceflat=forceflat, + ) @property def nproc(self) -> int: diff --git a/pylops/linearoperator.py b/pylops/linearoperator.py index c3d793a1..e397dfc3 100644 --- a/pylops/linearoperator.py +++ b/pylops/linearoperator.py @@ -281,7 +281,10 @@ def forceflat(self): @forceflat.setter def forceflat(self, new_forceflat: bool) -> None: - self._forceflat = bool(new_forceflat) + # note that this can also be None so we check before forcing bool + self._forceflat = ( + new_forceflat if new_forceflat is None else bool(new_forceflat) + ) @forceflat.deleter def forceflat(self): From 6626213af80aec9bd6d6114ec9a55495ae9700a3 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Mon, 27 Mar 2023 16:07:34 +0300 Subject: [PATCH 06/69] feature: added forceflat to HStack and VStack --- pylops/basicoperators/hstack.py | 13 ++++++++++++- pylops/basicoperators/vstack.py | 14 ++++++++++++-- pylops/linearoperator.py | 5 ++++- 3 files changed, 28 insertions(+), 4 deletions(-) diff --git a/pylops/basicoperators/hstack.py b/pylops/basicoperators/hstack.py index c97fac1c..775ecace 100644 --- a/pylops/basicoperators/hstack.py +++ b/pylops/basicoperators/hstack.py @@ -44,6 +44,10 @@ class HStack(LinearOperator): nproc : :obj:`int`, optional Number of processes used to evaluate the N operators in parallel using ``multiprocessing``. If ``nproc=1``, work in serial mode. + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after matvec. dtype : :obj:`str`, optional Type of elements in input array. @@ -107,6 +111,7 @@ def __init__( self, ops: Sequence[LinearOperator], nproc: int = 1, + forceflat: bool = None, dtype: Optional[str] = None, ) -> None: self.ops = ops @@ -129,7 +134,13 @@ def __init__( dtype = _get_dtype(self.ops) if dtype is None else np.dtype(dtype) clinear = all([getattr(oper, "clinear", True) for oper in self.ops]) - super().__init__(dtype=dtype, shape=(self.nops, self.mops), clinear=clinear) + super().__init__( + dtype=dtype, + shape=(self.nops, self.mops), + dimsd=ops[0].dimsd, + clinear=clinear, + forceflat=forceflat, + ) @property def nproc(self) -> int: diff --git a/pylops/basicoperators/vstack.py b/pylops/basicoperators/vstack.py index dc6734f8..d9d5d207 100644 --- a/pylops/basicoperators/vstack.py +++ b/pylops/basicoperators/vstack.py @@ -44,6 +44,10 @@ class VStack(LinearOperator): nproc : :obj:`int`, optional Number of processes used to evaluate the N operators in parallel using ``multiprocessing``. If ``nproc=1``, work in serial mode. + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after rmatvec. dtype : :obj:`str`, optional Type of elements in input array. @@ -107,6 +111,7 @@ def __init__( self, ops: Sequence[LinearOperator], nproc: int = 1, + forceflat: bool = None, dtype: Optional[DTypeLike] = None, ) -> None: self.ops = ops @@ -126,10 +131,15 @@ def __init__( self.pool = None if self.nproc > 1: self.pool = mp.Pool(processes=nproc) - dtype = _get_dtype(self.ops) if dtype is None else np.dtype(dtype) clinear = all([getattr(oper, "clinear", True) for oper in self.ops]) - super().__init__(dtype=dtype, shape=(self.nops, self.mops), clinear=clinear) + super().__init__( + dtype=dtype, + shape=(self.nops, self.mops), + dims=ops[0].dims, + clinear=clinear, + forceflat=forceflat, + ) @property def nproc(self) -> int: diff --git a/pylops/linearoperator.py b/pylops/linearoperator.py index c3d793a1..e397dfc3 100644 --- a/pylops/linearoperator.py +++ b/pylops/linearoperator.py @@ -281,7 +281,10 @@ def forceflat(self): @forceflat.setter def forceflat(self, new_forceflat: bool) -> None: - self._forceflat = bool(new_forceflat) + # note that this can also be None so we check before forcing bool + self._forceflat = ( + new_forceflat if new_forceflat is None else bool(new_forceflat) + ) @forceflat.deleter def forceflat(self): From f156ac433de44b893dcbe0a340888d084be969cc Mon Sep 17 00:00:00 2001 From: mrava87 Date: Tue, 28 Mar 2023 08:48:47 +0300 Subject: [PATCH 07/69] feature: set forceflat=True for HStack/VStack if Ops have different dimsd/dims --- pylops/basicoperators/hstack.py | 11 +++++++++-- pylops/basicoperators/vstack.py | 10 +++++++++- 2 files changed, 18 insertions(+), 3 deletions(-) diff --git a/pylops/basicoperators/hstack.py b/pylops/basicoperators/hstack.py index 775ecace..5cfbbec0 100644 --- a/pylops/basicoperators/hstack.py +++ b/pylops/basicoperators/hstack.py @@ -126,18 +126,25 @@ def __init__( raise ValueError("operators have different number of rows") self.nops = int(nops[0]) self.mmops = np.insert(np.cumsum(mops), 0, 0) + # define dimsd (check if all operators have the same, + # otherwise make same as self.nops and forceflat=True) + dimsd = [op.dimsd for op in self.ops] + if len(set(dimsd)) == 1: + dimsd = dimsd[0] + else: + dimsd = (self.nops,) + forceflat = True # create pool for multiprocessing self._nproc = nproc self.pool = None if self.nproc > 1: self.pool = mp.Pool(processes=nproc) - dtype = _get_dtype(self.ops) if dtype is None else np.dtype(dtype) clinear = all([getattr(oper, "clinear", True) for oper in self.ops]) super().__init__( dtype=dtype, shape=(self.nops, self.mops), - dimsd=ops[0].dimsd, + dimsd=dimsd, clinear=clinear, forceflat=forceflat, ) diff --git a/pylops/basicoperators/vstack.py b/pylops/basicoperators/vstack.py index d9d5d207..812b1a7e 100644 --- a/pylops/basicoperators/vstack.py +++ b/pylops/basicoperators/vstack.py @@ -126,6 +126,14 @@ def __init__( raise ValueError("operators have different number of columns") self.mops = int(mops[0]) self.nnops = np.insert(np.cumsum(nops), 0, 0) + # define dims (check if all operators have the same, + # otherwise make same as self.mops and forceflat=True) + dims = [op.dims for op in self.ops] + if len(set(dims)) == 1: + dims = dims[0] + else: + dims = (self.mops,) + forceflat = True # create pool for multiprocessing self._nproc = nproc self.pool = None @@ -136,7 +144,7 @@ def __init__( super().__init__( dtype=dtype, shape=(self.nops, self.mops), - dims=ops[0].dims, + dims=dims, clinear=clinear, forceflat=forceflat, ) From b5233773fb1967a2de89205c04c3ad46363cef22 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Tue, 28 Mar 2023 09:11:09 +0300 Subject: [PATCH 08/69] feature: added forceflat to BlockDiag --- pylops/basicoperators/blockdiag.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/pylops/basicoperators/blockdiag.py b/pylops/basicoperators/blockdiag.py index 2447ffcb..2b279ee1 100644 --- a/pylops/basicoperators/blockdiag.py +++ b/pylops/basicoperators/blockdiag.py @@ -44,6 +44,10 @@ class BlockDiag(LinearOperator): nproc : :obj:`int`, optional Number of processes used to evaluate the N operators in parallel using ``multiprocessing``. If ``nproc=1``, work in serial mode. + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after matvec and rmatvec. dtype : :obj:`str`, optional Type of elements in input array. @@ -108,6 +112,7 @@ def __init__( self, ops: Sequence[LinearOperator], nproc: int = 1, + forceflat: bool = None, dtype: Optional[DTypeLike] = None, ) -> None: self.ops = ops @@ -130,7 +135,12 @@ def __init__( dtype = _get_dtype(ops) if dtype is None else np.dtype(dtype) clinear = all([getattr(oper, "clinear", True) for oper in self.ops]) - super().__init__(dtype=dtype, shape=(self.nops, self.mops), clinear=clinear) + super().__init__( + dtype=dtype, + shape=(self.nops, self.mops), + clinear=clinear, + forceflat=forceflat, + ) @property def nproc(self) -> int: From bd59748c756451cf461bba118ebbb80fe1f0733f Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sat, 1 Apr 2023 14:02:39 +0300 Subject: [PATCH 09/69] doc: fixes https://github.com/PyLops/pylops/issues/503 --- tutorials/solvers.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tutorials/solvers.py b/tutorials/solvers.py index 018b5f85..c0bfe72c 100755 --- a/tutorials/solvers.py +++ b/tutorials/solvers.py @@ -178,8 +178,8 @@ # derivative operator # # .. math:: -# \mathbf{x} = (\mathbf{R^TR}+\epsilon_\nabla^2\nabla^T\nabla)^{-1} -# \mathbf{R^Ty} +# \mathbf{x} = (\mathbf{R}^T\mathbf{R}+\epsilon_\nabla^2\nabla^T\nabla)^{-1} +# \mathbf{R}^T\mathbf{y} # Create regularization operator D2op = pylops.SecondDerivative(N, dtype="float64") From 169bd8def0d0adb372381e6d5149dc5ce11298aa Mon Sep 17 00:00:00 2001 From: mrava87 Date: Fri, 7 Apr 2023 23:20:25 +0300 Subject: [PATCH 10/69] feature: Added forceflat to BlockDiag and MatrixMult --- pylops/basicoperators/blockdiag.py | 20 +++++++++++++++++++- pylops/basicoperators/matrixmult.py | 23 ++++++++++++++++++++++- pylops/basicoperators/sum.py | 17 ++++++++--------- 3 files changed, 49 insertions(+), 11 deletions(-) diff --git a/pylops/basicoperators/blockdiag.py b/pylops/basicoperators/blockdiag.py index 2b279ee1..16fd476e 100644 --- a/pylops/basicoperators/blockdiag.py +++ b/pylops/basicoperators/blockdiag.py @@ -127,6 +127,22 @@ def __init__( self.mops = int(mops.sum()) self.nnops = np.insert(np.cumsum(nops), 0, 0) self.mmops = np.insert(np.cumsum(mops), 0, 0) + # define dims (check if all operators have the same, + # otherwise make same as self.mops and forceflat=True) + dims = [op.dims for op in self.ops] + if len(set(dims)) == 1: + dims = (len(ops), *dims[0]) + else: + dims = (self.mops,) + forceflat = True + # define dimsd (check if all operators have the same, + # otherwise make same as self.nops and forceflat=True) + dimsd = [op.dimsd for op in self.ops] + if len(set(dimsd)) == 1: + dimsd = (len(ops), *dimsd[0]) + else: + dimsd = (self.nops,) + forceflat = True # create pool for multiprocessing self._nproc = nproc self.pool: Optional[mp.pool.Pool] = None @@ -137,10 +153,12 @@ def __init__( clinear = all([getattr(oper, "clinear", True) for oper in self.ops]) super().__init__( dtype=dtype, - shape=(self.nops, self.mops), + dims=dims, + dimsd=dimsd, clinear=clinear, forceflat=forceflat, ) + print(dims, dimsd) @property def nproc(self) -> int: diff --git a/pylops/basicoperators/matrixmult.py b/pylops/basicoperators/matrixmult.py index 65a62d4f..9aa44e73 100644 --- a/pylops/basicoperators/matrixmult.py +++ b/pylops/basicoperators/matrixmult.py @@ -29,6 +29,12 @@ class MatrixMult(LinearOperator): Number of samples for each other dimension of model (model/data will be reshaped and ``A`` applied multiple times to each column of the model/data). + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after matvec and rmatvec. Note that this is only + required when `otherdims=None`, otherwise pylops will detect whether to + return a 1d or nd array. dtype : :obj:`str`, optional Type of elements in input array. name : :obj:`str`, optional @@ -56,6 +62,7 @@ def __init__( self, A: NDArray, otherdims: Optional[Union[int, InputDimsLike]] = None, + forceflat: bool = None, dtype: DTypeLike = "float64", name: str = "M", ) -> None: @@ -81,12 +88,26 @@ def __init__( self.reshape = True explicit = False + # Check if forceflat is needed and set it back to None otherwise + if otherdims is not None and forceflat is not None: + logging.warning( + f"setting forceflat=None since otherdims!=None. " + f"PyLops will automatically detect whether to return " + f"a 1d or nd array based on the shape of the input" + f"array." + ) + forceflat = None # Check dtype for correctness (upcast to complex when A is complex) if np.iscomplexobj(A) and not np.iscomplexobj(np.ones(1, dtype=dtype)): dtype = A.dtype logging.warning("Matrix A is a complex object, dtype cast to %s" % dtype) super().__init__( - dtype=np.dtype(dtype), dims=dims, dimsd=dimsd, explicit=explicit, name=name + dtype=np.dtype(dtype), + dims=dims, + dimsd=dimsd, + explicit=explicit, + forceflat=forceflat, + name=name, ) def _matvec(self, x: NDArray) -> NDArray: diff --git a/pylops/basicoperators/sum.py b/pylops/basicoperators/sum.py index c3066514..855cf521 100644 --- a/pylops/basicoperators/sum.py +++ b/pylops/basicoperators/sum.py @@ -83,15 +83,14 @@ def __init__( dimsd = list(dims).copy() dimsd.pop(self.axis) # check if forceflat is needed and set it back to None otherwise - if len(dims) > 2: - if forceflat is not None: - logging.warning( - f"setting forceflat=None since len(dims)={len(dims)}>2. " - f"PyLops will automatically detect whether to return " - f"a 1d or nd array based on the shape of the input" - f"array." - ) - forceflat = None + if len(dims) > 2 and forceflat is not None: + logging.warning( + f"setting forceflat=None since len(dims)={len(dims)}>2. " + f"PyLops will automatically detect whether to return " + f"a 1d or nd array based on the shape of the input" + f"array." + ) + forceflat = None super().__init__( dtype=np.dtype(dtype), dims=dims, From 0ed9bff42688e8d8a28cf6f6fdf754fc5a78edca Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sat, 8 Apr 2023 21:39:39 +0300 Subject: [PATCH 11/69] feature: changed ISTA and FISTA to use matvec/rmatvec (or *mat) --- pylops/optimization/cls_sparsity.py | 46 +++++++++++++++++++---------- 1 file changed, 30 insertions(+), 16 deletions(-) diff --git a/pylops/optimization/cls_sparsity.py b/pylops/optimization/cls_sparsity.py index 80aae105..9c15f40f 100644 --- a/pylops/optimization/cls_sparsity.py +++ b/pylops/optimization/cls_sparsity.py @@ -1137,8 +1137,8 @@ def setup( Parameters ---------- y : :obj:`np.ndarray` - Data of size :math:`[N \times 1]` or :math:`[N \times R]` where for a solution for multiple right-hand-side - is found when ``R>1``. + Data of size :math:`[N \times 1]` or :math:`[N \times R]` where + a solution for multiple right-hand-side is found when ``R>1``. x0: :obj:`numpy.ndarray`, optional Initial guess niter : :obj:`int` @@ -1192,6 +1192,21 @@ def setup( self.ncp = get_array_module(y) + # choose matvec/rmatvec or matmat/rmatmat based on R + if y.ndim > 1 and y.shape[1] > 1: + self.Opmatvec = self.Op.matmat + self.Oprmatvec = self.Op.rmatmat + if self.SOp is not None: + self.SOpmatvec = self.SOp.matmat + self.SOprmatvec = self.SOp.rmatmat + else: + self.Opmatvec = self.Op.matvec + self.Oprmatvec = self.Op.rmatvec + if self.SOp is not None: + self.SOpmatvec = self.SOp.matvec + self.SOprmatvec = self.SOp.rmatvec + + # choose thresholding function if threshkind not in [ "hard", "soft", @@ -1214,7 +1229,6 @@ def setup( "soft-percentile, or half-percentile thresholding" ) - # choose thresholding function self.threshf: Callable[[NDArray, float], NDArray] if threshkind == "soft": self.threshf = _softthreshold @@ -1238,7 +1252,7 @@ def setup( self.alpha = alpha elif not hasattr(self, "alpha"): # compute largest eigenvalues of Op^H * Op - Op1 = self.Op.H * self.Op + Op1 = self.Op.H @ self.Op if get_module_name(self.ncp) == "numpy": maxeig: float = np.abs( Op1.eigs( @@ -1316,7 +1330,7 @@ def step(self, x: NDArray, show: bool = False) -> Tuple[NDArray, float]: # store old vector xold = x.copy() # compute residual - res: NDArray = self.y - self.Op @ x + res: NDArray = self.y - self.Opmatvec(x) if self.monitorres: self.normres = np.linalg.norm(res) if self.normres > self.normresold: @@ -1328,19 +1342,19 @@ def step(self, x: NDArray, show: bool = False) -> Tuple[NDArray, float]: else: self.normresold = self.normres - # compute gradient (must use .H @ for the solver to work with either single or multiple rhs) - grad: NDArray = self.alpha * (self.Op.H @ res) + # compute gradient + grad: NDArray = self.alpha * (self.Oprmatvec(res)) - # update inverted model (must use @ / .H @ for the solver to work with either single or multiple rhs) + # update inverted model x_unthesh: NDArray = x + grad if self.SOp is not None: - x_unthesh = self.SOp.H @ x_unthesh + x_unthesh = self.SOprmatvec(x_unthesh) if self.perc is None and self.decay is not None: x = self.threshf(x_unthesh, self.decay[self.iiter] * self.thresh) elif self.perc is not None: x = self.threshf(x_unthesh, 100 - self.perc) if self.SOp is not None: - x = self.SOp @ x + x = self.SOpmatvec(x) # model update xupdate = np.linalg.norm(x - xold) @@ -1594,7 +1608,7 @@ def step(self, x: NDArray, z: NDArray, show: bool = False) -> NDArray: # store old vector xold = x.copy() # compute residual - resz: NDArray = self.y - self.Op @ z + resz: NDArray = self.y - self.Opmatvec(z) if self.monitorres: self.normres = np.linalg.norm(resz) if self.normres > self.normresold: @@ -1606,19 +1620,19 @@ def step(self, x: NDArray, z: NDArray, show: bool = False) -> NDArray: else: self.normresold = self.normres - # compute gradient (must use .H @ for the solver to work with either single or multiple rhs - grad: NDArray = self.alpha * (self.Op.H @ resz) + # compute gradient + grad: NDArray = self.alpha * (self.Oprmatvec(resz)) - # update inverted model (must use @ / .H @ for the solver to work with either single or multiple rhs) + # update inverted model x_unthesh: NDArray = z + grad if self.SOp is not None: - x_unthesh = self.SOp.H @ x_unthesh + x_unthesh = self.SOprmatvec(x_unthesh) if self.perc is None and self.decay is not None: x = self.threshf(x_unthesh, self.decay[self.iiter] * self.thresh) elif self.perc is not None: x = self.threshf(x_unthesh, 100 - self.perc) if self.SOp is not None: - x = self.SOp @ x + x = self.SOpmatvec(x) # update auxiliary coefficients told = self.t From b3be22c618c302fc1705fdd67573f73c29f8ef11 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sat, 8 Apr 2023 21:41:03 +0300 Subject: [PATCH 12/69] minor: fixed flake8 --- pylops/basicoperators/matrixmult.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pylops/basicoperators/matrixmult.py b/pylops/basicoperators/matrixmult.py index 9aa44e73..a5f713ab 100644 --- a/pylops/basicoperators/matrixmult.py +++ b/pylops/basicoperators/matrixmult.py @@ -91,10 +91,10 @@ def __init__( # Check if forceflat is needed and set it back to None otherwise if otherdims is not None and forceflat is not None: logging.warning( - f"setting forceflat=None since otherdims!=None. " - f"PyLops will automatically detect whether to return " - f"a 1d or nd array based on the shape of the input" - f"array." + "setting forceflat=None since otherdims!=None. " + "PyLops will automatically detect whether to return " + "a 1d or nd array based on the shape of the input " + "array." ) forceflat = None # Check dtype for correctness (upcast to complex when A is complex) From 6180b5afdcefb4ca34cb1a8c0036b0f2265fea56 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sat, 8 Apr 2023 22:32:39 +0300 Subject: [PATCH 13/69] feature: added ndarray capabilities to Identity --- pylops/basicoperators/identity.py | 111 +++++++++++++++++++++++++----- 1 file changed, 93 insertions(+), 18 deletions(-) diff --git a/pylops/basicoperators/identity.py b/pylops/basicoperators/identity.py index a41f5ce9..2c6965d2 100644 --- a/pylops/basicoperators/identity.py +++ b/pylops/basicoperators/identity.py @@ -1,13 +1,14 @@ __all__ = ["Identity"] -from typing import Optional +from typing import Optional, Union import numpy as np from pylops import LinearOperator from pylops.utils.backend import get_array_module -from pylops.utils.typing import DTypeLike, NDArray +from pylops.utils.decorators import reshaped +from pylops.utils.typing import DTypeLike, InputDimsLike, NDArray class Identity(LinearOperator): @@ -19,16 +20,29 @@ class Identity(LinearOperator): removes last :math:`N - M` elements from data in adjoint and pads with :math:`0` in forward. + Note that the identity operator can handle both 1d and nd arrays; in the + case of nd arrays, all elements of N must be larger or equal than those of M + (or all elements of M must be larger or equal than those of N). + Parameters ---------- - N : :obj:`int` + N : :obj:`int` or :obj:`tuple` Number of samples in data (and model, if ``M`` is not provided). - M : :obj:`int`, optional - Number of samples in model. + If a tuple is provided, this is interpreted as the data (and model) + are nd-arrays. + M : :obj:`int` or :obj:`tuple`, optional + Number of samples in model. If a tuple is provided, this is interpreted + as the model is an nd-array. Note that when `M` is a tuple, `N` must be + also a tuple with the same number of elements. inplace : :obj:`bool`, optional Work inplace (``True``) or make a new copy (``False``). By default, data is a reference to the model (in forward) and model is a reference to the data (in adjoint). + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after matvec and rmatvec. Note that this is only + required when `N` and `M` are tuples (input and output arrays are nd-arrays. dtype : :obj:`str`, optional Type of elements in input array. name : :obj:`str`, optional @@ -44,6 +58,15 @@ class Identity(LinearOperator): Operator contains a matrix that can be solved explicitly (``True``) or not (``False``) + Raises + ------ + ValueError + - If ``M`` is a tuple with different number of elements of ``N`` + - If ``N`` ``M`` are non-identical tuples and some values are largers + in ``N`` and some in ``M`` + NotImplemented + If ``N`` or ``M`` have type different from int or tuple/list + Notes ----- For :math:`M = N`, an *Identity* operator simply moves the model @@ -87,38 +110,90 @@ class Identity(LinearOperator): def __init__( self, - N: int, - M: Optional[int] = None, + N: Union[int, InputDimsLike], + M: Optional[Union[int, InputDimsLike]] = None, inplace: bool = True, + forceflat: bool = None, dtype: DTypeLike = "float64", name: str = "I", ) -> None: M = N if M is None else M - super().__init__(dtype=np.dtype(dtype), shape=(N, M), name=name) + if isinstance(N, int) and isinstance(M, int): + # N and M are scalars (1d-arrays) + super().__init__( + dtype=np.dtype(dtype), + dims=(M,), + dimsd=(N,), + forceflat=forceflat, + name=name, + ) + # identify behaviour for matvec/rmatvec: 'same' for N=M, + # 'data' for N>M, and 'model' for M>N + if N == M: + self.mode = "same" + elif N < M: + self.mode = "model" + self.sliceN = slice(0, N) + self.sliceM = slice(0, M) + else: + self.mode = "data" + self.sliceN = slice(0, N) + self.sliceM = slice(0, M) + elif isinstance(N, (tuple, list)) and isinstance(M, (tuple, list)): + # N and M are tuples (nd-arrays) + # First check that all elements in N and M are the same or that + # all elements of either N or M are bigger than the other one and + # raise error is not the case + if np.array_equal(N, M): + self.mode = "same" + elif np.array_equal(M, np.maximum(N, M)): + self.mode = "model" + self.sliceN = tuple([slice(0, n) for n in N]) + self.sliceM = tuple([slice(0, m) for m in M]) + elif np.array_equal(N, np.maximum(N, M)): + self.mode = "data" + self.sliceN = tuple([slice(0, n) for n in N]) + self.sliceM = tuple([slice(0, m) for m in M]) + else: + raise ValueError( + "N and M are not identical, " + "and some values are larger in N and some in M" + ) + super().__init__( + dtype=np.dtype(dtype), dims=M, dimsd=N, forceflat=forceflat, name=name + ) + else: + raise NotImplemented( + f"N and M must have same type and equal to " + f"int, tuple, or list, instead their types" + f" are type(N)={type(N)} and type(M)={type(M)}" + ) self.inplace = inplace + @reshaped def _matvec(self, x: NDArray) -> NDArray: ncp = get_array_module(x) if not self.inplace: x = x.copy() - if self.shape[0] == self.shape[1]: + if self.mode == "same": y = x - elif self.shape[0] < self.shape[1]: - y = x[: self.shape[0]] + elif self.mode == "model": + y = x[self.sliceN] else: - y = ncp.zeros(self.shape[0], dtype=self.dtype) - y[: self.shape[1]] = x + y = ncp.zeros(self.dimsd, dtype=self.dtype) + y[self.sliceM] = x return y + @reshaped def _rmatvec(self, x: NDArray) -> NDArray: ncp = get_array_module(x) if not self.inplace: x = x.copy() - if self.shape[0] == self.shape[1]: + if self.mode == "same": y = x - elif self.shape[0] < self.shape[1]: - y = ncp.zeros(self.shape[1], dtype=self.dtype) - y[: self.shape[0]] = x + elif self.mode == "model": + y = ncp.zeros(self.dims, dtype=self.dtype) + y[self.sliceN] = x else: - y = x[: self.shape[1]] + y = x[self.sliceM] return y From d37d3a4f51bfe06920cbe994b5b07c87c777edc7 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sat, 8 Apr 2023 22:35:45 +0300 Subject: [PATCH 14/69] minor: fix flake8 in Identity --- pylops/basicoperators/identity.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pylops/basicoperators/identity.py b/pylops/basicoperators/identity.py index 2c6965d2..64e17eb4 100644 --- a/pylops/basicoperators/identity.py +++ b/pylops/basicoperators/identity.py @@ -64,7 +64,7 @@ class Identity(LinearOperator): - If ``M`` is a tuple with different number of elements of ``N`` - If ``N`` ``M`` are non-identical tuples and some values are largers in ``N`` and some in ``M`` - NotImplemented + NotImplementedError If ``N`` or ``M`` have type different from int or tuple/list Notes @@ -163,7 +163,7 @@ def __init__( dtype=np.dtype(dtype), dims=M, dimsd=N, forceflat=forceflat, name=name ) else: - raise NotImplemented( + raise NotImplementedError( f"N and M must have same type and equal to " f"int, tuple, or list, instead their types" f" are type(N)={type(N)} and type(M)={type(M)}" From a9afeb35b7ab23be223b0b4c9f10db75025056e1 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sat, 15 Apr 2023 22:10:24 +0300 Subject: [PATCH 15/69] feature: added forceflat to Block --- pylops/basicoperators/block.py | 45 +++++++++++++++++++++++++--------- 1 file changed, 33 insertions(+), 12 deletions(-) diff --git a/pylops/basicoperators/block.py b/pylops/basicoperators/block.py index 88e989bd..15eea4a8 100644 --- a/pylops/basicoperators/block.py +++ b/pylops/basicoperators/block.py @@ -1,7 +1,6 @@ __all__ = ["Block"] import multiprocessing as mp - from typing import Iterable, Optional from pylops import LinearOperator @@ -15,12 +14,18 @@ class _Block(LinearOperator): Used to be able to provide operators from different libraries to Block. """ - def __init__(self, ops: Iterable[Iterable[LinearOperator]], - dtype: Optional[DTypeLike] = None, - _HStack=HStack, - _VStack=VStack, - args_HStack: Optional[dict] = None, - args_VStack: Optional[dict] = None, name: str = 'B'): + + def __init__( + self, + ops: Iterable[Iterable[LinearOperator]], + forceflat: bool = None, + dtype: Optional[DTypeLike] = None, + _HStack=HStack, + _VStack=VStack, + args_HStack: Optional[dict] = None, + args_VStack: Optional[dict] = None, + name: str = "B", + ): if args_HStack is None: self.args_HStack = {} else: @@ -30,7 +35,12 @@ def __init__(self, ops: Iterable[Iterable[LinearOperator]], else: self.args_VStack = args_VStack hblocks = [_HStack(hblock, dtype=dtype, **self.args_HStack) for hblock in ops] - super().__init__(Op=_VStack(ops=hblocks, dtype=dtype, **self.args_VStack), name=name) + super().__init__( + Op=_VStack( + ops=hblocks, forceflat=forceflat, dtype=dtype, **self.args_VStack + ), + name=name, + ) def _matvec(self, x: NDArray) -> NDArray: return super()._matvec(x) @@ -53,6 +63,10 @@ class Block(_Block): nproc : :obj:`int`, optional Number of processes used to evaluate the N operators in parallel using ``multiprocessing``. If ``nproc=1``, work in serial mode. + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after rmatvec. dtype : :obj:`str`, optional Type of elements in input array. @@ -123,9 +137,16 @@ class Block(_Block): \end{bmatrix} """ - def __init__(self, ops: Iterable[Iterable[LinearOperator]], - nproc: int = 1, - dtype: Optional[DTypeLike] = None): + + def __init__( + self, + ops: Iterable[Iterable[LinearOperator]], + nproc: int = 1, + forceflat: bool = None, + dtype: Optional[DTypeLike] = None, + ): if nproc > 1: self.pool = mp.Pool(processes=nproc) - super().__init__(ops=ops, dtype=dtype, args_VStack={"nproc": nproc}) + super().__init__( + ops=ops, forceflat=forceflat, dtype=dtype, args_VStack={"nproc": nproc} + ) From 200973f9cc5dcc8b329033a4b57694c6647a5e97 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sat, 15 Apr 2023 23:36:40 +0300 Subject: [PATCH 16/69] feature: allow NDarray for Zero --- pylops/basicoperators/zero.py | 42 +++++++++++++++++++++++++++-------- 1 file changed, 33 insertions(+), 9 deletions(-) diff --git a/pylops/basicoperators/zero.py b/pylops/basicoperators/zero.py index 04e5f12c..68100c1f 100644 --- a/pylops/basicoperators/zero.py +++ b/pylops/basicoperators/zero.py @@ -1,12 +1,12 @@ __all__ = ["Zero"] -from typing import Optional +from typing import Optional, Union import numpy as np from pylops import LinearOperator from pylops.utils.backend import get_array_module -from pylops.utils.typing import DTypeLike, NDArray +from pylops.utils.typing import DTypeLike, InputDimsLike, NDArray class Zero(LinearOperator): @@ -17,10 +17,14 @@ class Zero(LinearOperator): Parameters ---------- - N : :obj:`int` - Number of samples in data (and model in M is not provided). - M : :obj:`int`, optional - Number of samples in model. + N : :obj:`int` or :obj:`tuple` + Number of samples in data (and model, if ``M`` is not provided). + If a tuple is provided, this is interpreted as the data (and model) + are nd-arrays. + M : :obj:`int` or :obj:`tuple`, optional + Number of samples in model. If a tuple is provided, this is interpreted + as the model is an nd-array. Note that when `M` is a tuple, `N` must be + also a tuple with the same number of elements. dtype : :obj:`str`, optional Type of elements in input array. name : :obj:`str`, optional @@ -55,13 +59,33 @@ class Zero(LinearOperator): def __init__( self, - N: int, - M: Optional[int] = None, + N: Union[int, InputDimsLike], + M: Optional[Union[int, InputDimsLike]] = None, + forceflat: bool = None, dtype: DTypeLike = "float64", name: str = "Z", ) -> None: M = N if M is None else M - super().__init__(dtype=np.dtype(dtype), shape=(N, M), name=name) + if isinstance(N, int) and isinstance(M, int): + # N and M are scalars (1d-arrays) + super().__init__( + dtype=np.dtype(dtype), + dims=(M,), + dimsd=(N,), + forceflat=forceflat, + name=name, + ) + elif isinstance(N, (tuple, list)) and isinstance(M, (tuple, list)): + # N and M are tuples (nd-arrays) + super().__init__( + dtype=np.dtype(dtype), dims=M, dimsd=N, forceflat=forceflat, name=name + ) + else: + raise NotImplementedError( + f"N and M must have same type and equal to " + f"int, tuple, or list, instead their types" + f" are type(N)={type(N)} and type(M)={type(M)}" + ) def _matvec(self, x: NDArray) -> NDArray: ncp = get_array_module(x) From 114be117f3f42c191fc3631748f265174001577b Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sun, 16 Apr 2023 09:41:04 +0300 Subject: [PATCH 17/69] minor: reduced redundant code in Zero init --- pylops/basicoperators/zero.py | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/pylops/basicoperators/zero.py b/pylops/basicoperators/zero.py index 68100c1f..23839a9d 100644 --- a/pylops/basicoperators/zero.py +++ b/pylops/basicoperators/zero.py @@ -68,24 +68,23 @@ def __init__( M = N if M is None else M if isinstance(N, int) and isinstance(M, int): # N and M are scalars (1d-arrays) - super().__init__( - dtype=np.dtype(dtype), - dims=(M,), - dimsd=(N,), - forceflat=forceflat, - name=name, - ) + dims, dimsd = (M,), (N,) elif isinstance(N, (tuple, list)) and isinstance(M, (tuple, list)): # N and M are tuples (nd-arrays) - super().__init__( - dtype=np.dtype(dtype), dims=M, dimsd=N, forceflat=forceflat, name=name - ) + dims, dimsd = M, N else: raise NotImplementedError( f"N and M must have same type and equal to " f"int, tuple, or list, instead their types" f" are type(N)={type(N)} and type(M)={type(M)}" ) + super().__init__( + dtype=np.dtype(dtype), + dims=dims, + dimsd=dimsd, + forceflat=forceflat, + name=name, + ) def _matvec(self, x: NDArray) -> NDArray: ncp = get_array_module(x) From cfcdb55b2aa7faf127da29c99ab2d695346f3cc7 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sun, 16 Apr 2023 09:50:15 +0300 Subject: [PATCH 18/69] minor: added forceflat to docstring of Zero --- pylops/basicoperators/identity.py | 2 +- pylops/basicoperators/zero.py | 5 +++++ 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/pylops/basicoperators/identity.py b/pylops/basicoperators/identity.py index 64e17eb4..c2d05a30 100644 --- a/pylops/basicoperators/identity.py +++ b/pylops/basicoperators/identity.py @@ -42,7 +42,7 @@ class Identity(LinearOperator): .. versionadded:: 2.2.0 Force an array to be flattened after matvec and rmatvec. Note that this is only - required when `N` and `M` are tuples (input and output arrays are nd-arrays. + required when `N` and `M` are tuples (input and output arrays are nd-arrays). dtype : :obj:`str`, optional Type of elements in input array. name : :obj:`str`, optional diff --git a/pylops/basicoperators/zero.py b/pylops/basicoperators/zero.py index 23839a9d..eb375d9a 100644 --- a/pylops/basicoperators/zero.py +++ b/pylops/basicoperators/zero.py @@ -25,6 +25,11 @@ class Zero(LinearOperator): Number of samples in model. If a tuple is provided, this is interpreted as the model is an nd-array. Note that when `M` is a tuple, `N` must be also a tuple with the same number of elements. + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after matvec and rmatvec. Note that this is only + required when `N` and `M` are tuples (input and output arrays are nd-arrays). dtype : :obj:`str`, optional Type of elements in input array. name : :obj:`str`, optional From 89e55d3ee72892537982bb9371c723003edcdb35 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sun, 16 Apr 2023 11:09:24 +0300 Subject: [PATCH 19/69] minor: remove unwanted print from BlockDiag --- pylops/basicoperators/blockdiag.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pylops/basicoperators/blockdiag.py b/pylops/basicoperators/blockdiag.py index 16fd476e..e13ed026 100644 --- a/pylops/basicoperators/blockdiag.py +++ b/pylops/basicoperators/blockdiag.py @@ -158,7 +158,6 @@ def __init__( clinear=clinear, forceflat=forceflat, ) - print(dims, dimsd) @property def nproc(self) -> int: From c18002463e2c68e604070340588988a5b3a7ffe2 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Tue, 2 May 2023 20:13:47 +0200 Subject: [PATCH 20/69] feature: added conv with long filters --- pylops/signalprocessing/convolve1d.py | 248 +++++++++++++++++++------- 1 file changed, 184 insertions(+), 64 deletions(-) diff --git a/pylops/signalprocessing/convolve1d.py b/pylops/signalprocessing/convolve1d.py index 9921ccf9..ffac372b 100644 --- a/pylops/signalprocessing/convolve1d.py +++ b/pylops/signalprocessing/convolve1d.py @@ -1,5 +1,6 @@ __all__ = ["Convolve1D"] +from functools import partial from typing import Callable, Tuple, Union import numpy as np @@ -8,6 +9,7 @@ from pylops import LinearOperator from pylops.utils._internal import _value_or_sized_to_tuple from pylops.utils.backend import ( + get_array_module, get_convolve, get_fftconvolve, get_oaconvolve, @@ -18,7 +20,10 @@ def _choose_convfunc( - x: npt.ArrayLike, method: Union[None, str], dims + x: npt.ArrayLike, + method: Union[None, str], + dims: Union[int, InputDimsLike], + axis: int = -1, ) -> Tuple[Callable, str]: """Choose convolution function @@ -34,15 +39,181 @@ def _choose_convfunc( if method is None: method = "fft" if method == "fft": - convfunc = get_fftconvolve(x) + convfunc = partial(get_fftconvolve(x), axes=axis) elif method == "overlapadd": - convfunc = get_oaconvolve(x) + convfunc = partial(get_oaconvolve(x), axes=axis)(x) else: raise NotImplementedError("method must be fft or overlapadd") return convfunc, method -class Convolve1D(LinearOperator): +def _pad_along_axis(array: np.ndarray, pad_size: tuple, axis: int = 0) -> np.ndarray: + + npad = [(0, 0)] * array.ndim + npad[axis] = pad_size + return np.pad(array, pad_width=npad) + + +class _Convolve1Dshort(LinearOperator): + r"""1D convolution operator with compact filter (shorter than input signal)""" + + def __init__( + self, + dims: Union[int, InputDimsLike], + h: NDArray, + offset: int = 0, + axis: int = -1, + method: str = None, + dtype: DTypeLike = "float64", + name: str = "C", + ) -> None: + dims = _value_or_sized_to_tuple(dims) + super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dims, name=name) + self.axis = axis + self.nh = h.size if h.ndim == 1 else h.shape[axis] + if offset > self.nh - 1: + raise ValueError("offset must be smaller than h.shape[axis] - 1") + self.h = h + self.offset = 2 * (self.nh // 2 - int(offset)) + if self.nh % 2 == 0: + self.offset -= 1 + if self.offset != 0: + self.h = _pad_along_axis( + self.h, + ( + self.offset if self.offset > 0 else 0, + -self.offset if self.offset < 0 else 0, + ), + axis=-1 if h.ndim == 1 else axis, + ) + self.hstar = np.flip(self.h, axis=-1) + + # add dimensions to filter to match dimensions of model and data + if self.h.ndim == 1: + hdims = np.ones(len(self.dims), dtype=int) + hdims[self.axis] = len(self.h) + self.h = self.h.reshape(hdims) + self.hstar = self.hstar.reshape(hdims) + + # choose method and function handle + self.convfunc, self.method = _choose_convfunc(h, method, self.dims, self.axis) + + @reshaped + def _matvec(self, x: NDArray) -> NDArray: + if type(self.h) != type(x): + self.h = to_cupy_conditional(x, self.h) + self.convfunc, self.method = _choose_convfunc( + self.h, self.method, self.dims + ) + return self.convfunc(x, self.h, mode="same") + + @reshaped + def _rmatvec(self, x: NDArray) -> NDArray: + if type(self.hstar) != type(x): + self.hstar = to_cupy_conditional(x, self.hstar) + self.convfunc, self.method = _choose_convfunc( + self.hstar, self.method, self.dims + ) + return self.convfunc(x, self.hstar, mode="same") + + +class _Convolve1Dlong(LinearOperator): + """1D convolution operator with extended filter (larger than input signal)""" + + def __init__( + self, + dims: Union[int, InputDimsLike], + h: NDArray, + offset: int = 0, + axis: int = -1, + method: str = None, + dtype: DTypeLike = "float64", + name: str = "C", + ) -> None: + dims = _value_or_sized_to_tuple(dims) + dimsd = h.shape + super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dimsd, name=name) + + # create filter + self.axis = axis + if offset > self.dims[self.axis] - 1: + raise ValueError("offset must be smaller than self.dims[self.axis] - 1") + self.nh = h.size if h.ndim == 1 else h.shape[axis] + self.h = h + self.offset = 2 * (self.dims[self.axis] // 2 - int(offset)) + if self.dims[self.axis] % 2 == 0: + self.offset -= 1 + self.hstar = np.flip(self.h, axis=-1) + + self.pad = np.zeros((len(dims), 2), dtype=int) + self.pad[self.axis, 0] = self.offset if self.offset > 0 else 0 + self.pad[self.axis, 1] = -self.offset if self.offset < 0 else 0 + + self.padd = np.zeros((len(dims), 2), dtype=int) + self.padd[self.axis, 1] = self.offset if self.offset > 0 else 0 + self.padd[self.axis, 0] = -self.offset if self.offset < 0 else 0 + + # add dimensions to filter to match dimensions of model and data + if self.h.ndim == 1: + hdims = np.ones(len(self.dims), dtype=int) + hdims[self.axis] = len(self.h) + self.h = self.h.reshape(hdims) + self.hstar = self.hstar.reshape(hdims) + + # choose method and function handle + self.convfunc, self.method = _choose_convfunc(h, method, self.dims, self.axis) + + @reshaped + def _matvec(self, x: NDArray) -> NDArray: + if type(self.h) != type(x): + self.h = to_cupy_conditional(x, self.h) + self.convfunc, self.method = _choose_convfunc( + self.h, self.method, self.dims + ) + x = np.pad(x, self.pad) + y = self.convfunc(self.h, x, mode="same") + return y + + @reshaped + def _rmatvec(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + if type(self.h) != type(x): + self.hstar = to_cupy_conditional(x, self.hstar) + self.convfunc, self.method = _choose_convfunc( + self.hstar, self.method, self.dims + ) + x = np.pad(x, self.padd) + y = self.convfunc(self.hstar, x) + if self.dims[self.axis] % 2 == 0: + y = ncp.take( + y, + range( + len(y) // 2 - self.dims[self.axis] // 2, + len(y) // 2 + self.dims[self.axis] // 2, + ), + axis=self.axis, + ) + else: + y = ncp.take( + y, + range( + len(y) // 2 - self.dims[self.axis] // 2, + len(y) // 2 + self.dims[self.axis] // 2 + 1, + ), + axis=self.axis, + ) + return y + + +def Convolve1D( + dims: Union[int, InputDimsLike], + h: NDArray, + offset: int = 0, + axis: int = -1, + method: str = None, + dtype: DTypeLike = "float64", + name: str = "C", +) -> LinearOperator: r"""1D convolution operator. Apply one-dimensional convolution with a compact filter to model (and data) @@ -126,63 +297,12 @@ class Convolve1D(LinearOperator): y(t) = \mathscr{F}^{-1} (H(f)^* * X(f)) """ - - def __init__( - self, - dims: Union[int, InputDimsLike], - h: NDArray, - offset: int = 0, - axis: int = -1, - method: str = None, - dtype: DTypeLike = "float64", - name: str = "C", - ) -> None: - dims = _value_or_sized_to_tuple(dims) - super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dims, name=name) - - self.axis = axis - if offset > len(h) - 1: - raise ValueError("offset must be smaller than len(h) - 1") - self.h = h - self.hstar = np.flip(self.h, axis=-1) - self.nh = len(h) - self.offset = 2 * (self.nh // 2 - int(offset)) - if self.nh % 2 == 0: - self.offset -= 1 - if self.offset != 0: - self.h = np.pad( - self.h, - ( - self.offset if self.offset > 0 else 0, - -self.offset if self.offset < 0 else 0, - ), - mode="constant", - ) - self.hstar = np.flip(self.h, axis=-1) - - # add dimensions to filter to match dimensions of model and data - hdims = np.ones(len(self.dims), dtype=int) - hdims[self.axis] = len(self.h) - self.h = self.h.reshape(hdims) - self.hstar = self.hstar.reshape(hdims) - - # choose method and function handle - self.convfunc, self.method = _choose_convfunc(h, method, self.dims) - - @reshaped - def _matvec(self, x: NDArray) -> NDArray: - if type(self.h) is not type(x): - self.h = to_cupy_conditional(x, self.h) - self.convfunc, self.method = _choose_convfunc( - self.h, self.method, self.dims - ) - return self.convfunc(x, self.h, mode="same") - - @reshaped - def _rmatvec(self, x: NDArray) -> NDArray: - if type(self.hstar) is not type(x): - self.hstar = to_cupy_conditional(x, self.hstar) - self.convfunc, self.method = _choose_convfunc( - self.hstar, self.method, self.dims - ) - return self.convfunc(x, self.hstar, mode="same") + nh = h.size if h.ndim == 1 else h.shape[axis] + if nh <= _value_or_sized_to_tuple(dims)[axis]: + convop = _Convolve1Dshort + else: + convop = _Convolve1Dlong + c = convop( + dims=dims, h=h, offset=offset, axis=axis, method=method, dtype=dtype, name=name + ) + return c From b410a9fcda3ede0544d580a95fa7cbca22328759 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Tue, 2 May 2023 20:27:33 +0200 Subject: [PATCH 21/69] minor: converted Convolve1d into a class --- pylops/signalprocessing/convolve1d.py | 50 +++++++++++++++++---------- 1 file changed, 32 insertions(+), 18 deletions(-) diff --git a/pylops/signalprocessing/convolve1d.py b/pylops/signalprocessing/convolve1d.py index ffac372b..18cac5a9 100644 --- a/pylops/signalprocessing/convolve1d.py +++ b/pylops/signalprocessing/convolve1d.py @@ -205,15 +205,7 @@ def _rmatvec(self, x: NDArray) -> NDArray: return y -def Convolve1D( - dims: Union[int, InputDimsLike], - h: NDArray, - offset: int = 0, - axis: int = -1, - method: str = None, - dtype: DTypeLike = "float64", - name: str = "C", -) -> LinearOperator: +class Convolve1D(LinearOperator): r"""1D convolution operator. Apply one-dimensional convolution with a compact filter to model (and data) @@ -297,12 +289,34 @@ def Convolve1D( y(t) = \mathscr{F}^{-1} (H(f)^* * X(f)) """ - nh = h.size if h.ndim == 1 else h.shape[axis] - if nh <= _value_or_sized_to_tuple(dims)[axis]: - convop = _Convolve1Dshort - else: - convop = _Convolve1Dlong - c = convop( - dims=dims, h=h, offset=offset, axis=axis, method=method, dtype=dtype, name=name - ) - return c + + def __init__( + self, + dims: Union[int, InputDimsLike], + h: NDArray, + offset: int = 0, + axis: int = -1, + method: str = None, + dtype: DTypeLike = "float64", + name: str = "C", + ) -> None: + nh = h.size if h.ndim == 1 else h.shape[axis] + if nh <= _value_or_sized_to_tuple(dims)[axis]: + convop = _Convolve1Dshort + else: + convop = _Convolve1Dlong + Op = convop( + dims=dims, + h=h, + offset=offset, + axis=axis, + method=method, + dtype=dtype, + ) + super().__init__(Op=Op, name=name) + + def _matvec(self, x: NDArray) -> NDArray: + return super()._matvec(x) + + def _rmatvec(self, x: NDArray) -> NDArray: + return super()._rmatvec(x) From a0fdd0633f97a88b1b601f05400b9fe8c07f7531 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Tue, 2 May 2023 20:38:14 +0200 Subject: [PATCH 22/69] test: shorten wavelet to have nt>ntwav --- pytests/test_kirchhoff.py | 2 +- pytests/test_lsm.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/pytests/test_kirchhoff.py b/pytests/test_kirchhoff.py index 02deeafd..317821c3 100755 --- a/pytests/test_kirchhoff.py +++ b/pytests/test_kirchhoff.py @@ -49,7 +49,7 @@ r2d = np.vstack((rx, 2 * np.ones(PAR["nrx"]))) r3d = np.vstack((ryy.ravel(), rxx.ravel(), 2 * np.ones(PAR["nrx"] * PAR["nry"]))) -wav, _, wavc = ricker(t[:41], f0=40) +wav, _, wavc = ricker(t[:21], f0=40) par1 = {"mode": "analytic", "dynamic": False} par2 = {"mode": "eikonal", "dynamic": False} diff --git a/pytests/test_lsm.py b/pytests/test_lsm.py index 33024710..7f39cefd 100755 --- a/pytests/test_lsm.py +++ b/pytests/test_lsm.py @@ -48,7 +48,7 @@ r2d = np.vstack((rx, 2 * np.ones(PAR["nrx"]))) r3d = np.vstack((ryy.ravel(), rxx.ravel(), 2 * np.ones(PAR["nrx"] * PAR["nry"]))) -wav, _, wavc = ricker(t[:41], f0=40) +wav, _, wavc = ricker(t[:21], f0=40) par1 = {"mode": "analytic", "dynamic": False} par2 = {"mode": "eikonal", "dynamic": False} From 6478d56783cedf4e8d5cb1b84071ff5b30f47a50 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sun, 7 May 2023 22:33:35 +0300 Subject: [PATCH 23/69] doc: added example with convolution with large filter --- examples/plot_convolve.py | 117 +++++++++++++++++++++++++- pylops/signalprocessing/convolve1d.py | 13 +-- pytests/test_convolve.py | 39 +++++++-- 3 files changed, 156 insertions(+), 13 deletions(-) diff --git a/examples/plot_convolve.py b/examples/plot_convolve.py index 19431ff3..4afa9604 100644 --- a/examples/plot_convolve.py +++ b/examples/plot_convolve.py @@ -22,6 +22,7 @@ """ import matplotlib.pyplot as plt import numpy as np +import scipy as sp from scipy.sparse.linalg import lsqr import pylops @@ -132,7 +133,7 @@ plt.subplots_adjust(top=0.8) ############################################################################### -# Finally we do the same using three dimensional signals and +# We now do the same using three dimensional signals and # filters taking advantage of the # :py:class:`pylops.signalprocessing.ConvolveND` operator. ny, nx, nz = 13, 10, 7 @@ -149,7 +150,7 @@ xlsqr = xlsqr.reshape(Cop.dims) fig, axs = plt.subplots(3, 3, figsize=(10, 12)) -fig.suptitle("Convolve 3d data", y=0.95, fontsize=14, fontweight="bold") +fig.suptitle("Convolve 3d data", y=0.98, fontsize=14, fontweight="bold") axs[0][0].imshow(x[ny // 3], cmap="gray", vmin=-1, vmax=1) axs[0][1].imshow(y[ny // 3], cmap="gray", vmin=-1, vmax=1) axs[0][2].imshow(xlsqr[ny // 3], cmap="gray", vmin=-1, vmax=1) @@ -172,3 +173,115 @@ axs[2][1].axis("tight") axs[2][2].axis("tight") plt.tight_layout() + +############################################################################### +# Up until now, we have only considered the case where the filter is compact +# and much shorter of the input data. There are however scenarios where the +# opposite (i.e., filter is longer than the signal) is desired. This is for +# example the case when one wants to estimate a filter (:math:`\mathbf{h}`) +# to match two signals (:math:`\mathbf{x}` and :math:`\mathbf{y}`): +# +# .. math:: +# J = || \mathbf{y} - \mathbf{X} \mathbf{h} ||_2^2 +# +# Such a scenario is very commonly used in so-called adaptive substraction +# techniques. We will try now to use :py:class:`pylops.signalprocessing.Convolve1D` +# to match two signals that have both a phase and amplitude mismatch. + +# Define input signal +nt = 101 +dt = 0.004 +t = np.arange(nt) * dt + +x = np.zeros(nt) +x[[int(nt / 4), int(nt / 2), int(2 * nt / 3)]] = [3, -2, 1] +h, th, hcenter = ricker(t[:41], f0=20) + +Cop = pylops.signalprocessing.Convolve1D(nt, h=h, offset=hcenter, dtype="float32") +x = Cop * x + +# Phase and amplitude corrupt the input +amp = 0.9 +phase = 40 + +y = amp * ( + x * np.cos(np.deg2rad(phase)) + + np.imag(sp.signal.hilbert(x)) * np.sin(np.deg2rad(phase)) +) + +# Define convolution operator +nh = 21 +th = np.arange(nh) * dt - dt * nh // 2 +Yop = pylops.signalprocessing.Convolve1D(nh, h=y, offset=nh // 2) + +# Find filter to match x to y +h = Yop / x +ymatched = Yop @ h + +# Find sparse filter to match x to y +hsparse = pylops.optimization.sparsity.fista(Yop, x, niter=100, eps=1e-1)[0] +ymatchedsparse = Yop @ hsparse + +fig, ax = plt.subplots(1, 1, figsize=(10, 3)) +ax.plot(t, x, "k", lw=2, label=r"$x$") +ax.plot(t, y, "r", lw=2, label=r"$y$") +ax.plot(t, ymatched, "--g", lw=2, label=r"$y_{matched}$") +ax.plot(t, x - ymatched, "--k", lw=2, label=r"$x-y_{matched,sparse}$") +ax.plot(t, ymatchedsparse, "--m", lw=2, label=r"$y_{matched}$") +ax.plot(t, x - ymatchedsparse, "--k", lw=2, label=r"$x-y_{matched,sparse}$") +ax.set_title("Signals to match", fontsize=14, fontweight="bold") +ax.legend(loc="upper right") +plt.tight_layout() + +fig, axs = plt.subplots(1, 2, figsize=(10, 3)) +axs[0].plot(th, h, "k", lw=2) +axs[0].set_title("Matching filter", fontsize=14, fontweight="bold") +axs[1].plot(th, hsparse, "k", lw=2) +axs[1].set_title("Sparse Matching filter", fontsize=14, fontweight="bold") +plt.tight_layout() + + +############################################################################### +# Finally, in some cases one wants to convolve (or correlate) two signals of +# the same size. This can also be obtained using +# :py:class:`pylops.signalprocessing.Convolve1D`. We will see here a case +# where the operator is used to trace-wise auto-correlate signals from +# a 2-dimensional array representing a seismic dataset. + +# Create data +par = {"ox": -140, "dx": 2, "nx": 140, "ot": 0, "dt": 0.004, "nt": 200, "f0": 20} + +v = 1500 +t0 = [0.2, 0.4, 0.5] +px = [0, 0, 0] +pxx = [1e-5, 5e-6, 1e-20] +amp = [1.0, -2, 0.5] + +t, t2, x, y = pylops.utils.seismicevents.makeaxis(par) +wav = pylops.utils.wavelets.ricker(t[:41], f0=par["f0"])[0] +_, data = pylops.utils.seismicevents.parabolic2d(x, t, t0, px, pxx, amp, wav) + +# Convolution operator +Xop = pylops.signalprocessing.Convolve1D( + (par["nx"], par["nt"]), h=data, offset=par["nt"] // 2, axis=-1, method="fft" +) + +corr = Xop.H @ data + +fig, axs = plt.subplots(1, 2, figsize=(10, 6)) +axs[0].imshow(data.T, cmap="gray", vmin=-1, vmax=1, extent=(x[0], x[-1], t[-1], t[0])) +axs[0].set_xlabel("Rec (m)", fontsize=14, fontweight="bold") +axs[0].set_ylabel("T (s)", fontsize=14, fontweight="bold") +axs[0].set_title("Data", fontsize=14, fontweight="bold") +axs[0].axis("tight") +axs[1].imshow( + corr.T, + cmap="gray", + vmin=-10, + vmax=10, + extent=(x[0], x[-1], t[par["nt"] // 2], -t[par["nt"] // 2]), +) +axs[1].set_xlabel("Rec (m)", fontsize=14, fontweight="bold") +axs[1].set_title("Auto-correlation", fontsize=14, fontweight="bold") +axs[1].axis("tight") +plt.tight_layout() diff --git a/pylops/signalprocessing/convolve1d.py b/pylops/signalprocessing/convolve1d.py index 18cac5a9..8208d0ce 100644 --- a/pylops/signalprocessing/convolve1d.py +++ b/pylops/signalprocessing/convolve1d.py @@ -208,17 +208,18 @@ def _rmatvec(self, x: NDArray) -> NDArray: class Convolve1D(LinearOperator): r"""1D convolution operator. - Apply one-dimensional convolution with a compact filter to model (and data) - along an ``axis`` of a multi-dimensional array. + Apply one-dimensional convolution with i) a compact filter (shorter than input signal) or + ii) an extended filter (larger than input signal) to model (and data) along an ``axis`` + of a multi-dimensional array. Parameters ---------- dims : :obj:`list` or :obj:`int` - Number of samples for each dimension + Number of samples for each dimension of the model h : :obj:`numpy.ndarray` - 1d compact filter to be convolved to input signal + 1d filter to be convolved to input signal offset : :obj:`int` - Index of the center of the compact filter + Index of the center of the filter axis : :obj:`int`, optional .. versionadded:: 2.0.0 @@ -253,7 +254,7 @@ class Convolve1D(LinearOperator): Notes ----- The Convolve1D operator applies convolution between the input signal - :math:`x(t)` and a compact filter kernel :math:`h(t)` in forward model: + :math:`x(t)` and a filter kernel :math:`h(t)` in forward model: .. math:: y(t) = \int\limits_{-\infty}^{\infty} h(t-\tau) x(\tau) \,\mathrm{d}\tau diff --git a/pytests/test_convolve.py b/pytests/test_convolve.py index 4cddec63..d26938ba 100755 --- a/pytests/test_convolve.py +++ b/pytests/test_convolve.py @@ -133,7 +133,9 @@ def test_Convolve1D(par): x = np.zeros((par["nx"])) x[par["nx"] // 2] = 1.0 - xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr( + Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0 + )[0] assert_array_almost_equal(x, xlsqr, decimal=1) # 1D on 2D @@ -153,7 +155,9 @@ def test_Convolve1D(par): int(par["nx"] / 2 - 3) : int(par["nx"] / 2 + 3), ] = 1.0 x = x.ravel() - xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr( + Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0 + )[0] assert_array_almost_equal(x, xlsqr, decimal=1) # 1D on 3D @@ -175,10 +179,31 @@ def test_Convolve1D(par): int(par["nx"] / 2 - 3) : int(par["nx"] / 2 + 3), ] = 1.0 x = x.ravel() - xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0)[ + 0 + ] assert_array_almost_equal(x, xlsqr, decimal=1) +@pytest.mark.parametrize( + "par", [(par1_1d), (par2_1d), (par3_1d), (par4_1d), (par5_1d), (par6_1d)] +) +def test_Convolve1D_long(par): + """Dot-test and inversion for Convolve1D operator with long filter""" + np.random.seed(10) + # 1D + if par["axis"] == 0: + x = np.zeros((par["nx"])) + x[par["nx"] // 2] = 1.0 + Xop = Convolve1D(nfilt[0], h=x, offset=nfilt[0] // 2, dtype="float64") + assert dottest(Xop, par["nx"], nfilt[0]) + + h1lsqr = lsqr( + Xop, Xop * h1, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0 + )[0] + assert_array_almost_equal(h1, h1lsqr, decimal=1) + + @pytest.mark.parametrize( "par", [(par1_2d), (par2_2d), (par3_2d), (par4_2d), (par5_2d), (par6_2d)] ) @@ -200,7 +225,9 @@ def test_Convolve2D(par): int(par["nx"] / 2 - 3) : int(par["nx"] / 2 + 3), ] = 1.0 x = x.ravel() - xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr( + Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0 + )[0] assert_array_almost_equal(x, xlsqr, decimal=1) # 2D on 3D @@ -224,7 +251,9 @@ def test_Convolve2D(par): int(par["nx"] / 2 - 3) : int(par["nx"] / 2 + 3), ] = 1.0 x = x.ravel() - xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0)[ + 0 + ] # due to ringing in solution we cannot use assert_array_almost_equal assert np.linalg.norm(xlsqr - x) / np.linalg.norm(xlsqr) < 2e-1 From 33e72ec7f3b7f5f94e5842ee9cc9e67cc146a7b0 Mon Sep 17 00:00:00 2001 From: cako Date: Sun, 7 May 2023 19:31:34 -0700 Subject: [PATCH 24/69] feat: forceflat test for Sum --- pytests/test_basicoperators.py | 54 +++++++++++++++++++++++++++++++--- 1 file changed, 50 insertions(+), 4 deletions(-) diff --git a/pytests/test_basicoperators.py b/pytests/test_basicoperators.py index 02ae6946..43cb15f2 100755 --- a/pytests/test_basicoperators.py +++ b/pytests/test_basicoperators.py @@ -63,7 +63,9 @@ def test_LinearRegression(par): assert dottest(LRop, par["ny"], 2) x = np.array([1.0, 2.0], dtype=np.float32) - xlsqr = lsqr(LRop, LRop * x, damp=1e-10, iter_lim=300, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr( + LRop, LRop * x, damp=1e-10, iter_lim=300, atol=1e-8, btol=1e-8, show=0 + )[0] assert_array_almost_equal(x, xlsqr, decimal=3) y = LRop * x @@ -82,7 +84,9 @@ def test_MatrixMult(par): assert dottest(Gop, par["ny"], par["nx"], complexflag=0 if par["imag"] == 0 else 3) x = np.ones(par["nx"]) + par["imag"] * np.ones(par["nx"]) - xlsqr = lsqr(Gop, Gop * x, damp=1e-20, iter_lim=300, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr(Gop, Gop * x, damp=1e-20, iter_lim=300, atol=1e-8, btol=1e-8, show=0)[ + 0 + ] assert_array_almost_equal(x, xlsqr, decimal=4) @@ -100,7 +104,9 @@ def test_MatrixMult_sparse(par): assert dottest(Gop, par["ny"], par["nx"], complexflag=0 if par["imag"] == 1 else 3) x = np.ones(par["nx"]) + par["imag"] * np.ones(par["nx"]) - xlsqr = lsqr(Gop, Gop * x, damp=1e-20, iter_lim=300, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr(Gop, Gop * x, damp=1e-20, iter_lim=300, atol=1e-8, btol=1e-8, show=0)[ + 0 + ] assert_array_almost_equal(x, xlsqr, decimal=4) @@ -133,7 +139,9 @@ def test_MatrixMult_repeated(par): ) x = (np.ones((par["nx"], 5)) + par["imag"] * np.ones((par["nx"], 5))).ravel() - xlsqr = lsqr(Gop, Gop * x, damp=1e-20, iter_lim=300, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr(Gop, Gop * x, damp=1e-20, iter_lim=300, atol=1e-8, btol=1e-8, show=0)[ + 0 + ] assert_array_almost_equal(x, xlsqr, decimal=4) @@ -460,6 +468,44 @@ def test_Roll3D(par): assert_array_almost_equal(x[str(axis)].ravel(), xinv, decimal=3) +@pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j), (par3)]) +def test_Sum2D_forceflat(par): + np.random.seed(10) + flat_dimsd = par["ny"] + flat_dims = par["ny"] * par["nx"] + x = np.random.randn(flat_dims) + par["imag"] * np.random.randn(flat_dims) + + Sop_True = Sum((par["ny"], par["nx"]), axis=-1, forceflat=True) + y = Sop_True @ x + xadj = Sop_True.H @ y + assert y.shape == (flat_dimsd,) + assert xadj.shape == (flat_dims,) + + Sop_None = Sum((par["ny"], par["nx"]), axis=-1) + y = Sop_None @ x + xadj = Sop_None.H @ y + assert y.shape == (par["ny"],) + assert xadj.shape == (par["ny"], par["nx"]) + + Sop_False = Sum((par["ny"], par["nx"]), axis=-1, forceflat=False) + y = Sop_False @ x + xadj = Sop_False.H @ y + assert y.shape == (par["ny"],) + assert xadj.shape == (par["ny"], par["nx"]) + + with pytest.raises(ValueError): + Sop_True * Sop_False.H + + Sop = Sop_True * Sop_None.H + assert Sop.forceflat is True + + Sop = Sop_False * Sop_None.H + assert Sop.forceflat is False + + Sop = Sop_None * Sop_None.H + assert Sop.forceflat is None + + @pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j), (par3)]) def test_Sum2D(par): """Dot-test for Sum operator on 2d signal""" From 5ee83c4e78da8589e15ed24b2208d66e9e318dca Mon Sep 17 00:00:00 2001 From: cako Date: Sun, 7 May 2023 19:32:05 -0700 Subject: [PATCH 25/69] fix: use matvec everywhere --- pylops/optimization/cls_sparsity.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/pylops/optimization/cls_sparsity.py b/pylops/optimization/cls_sparsity.py index 9c15f40f..6ab4a30d 100644 --- a/pylops/optimization/cls_sparsity.py +++ b/pylops/optimization/cls_sparsity.py @@ -491,7 +491,7 @@ def step(self, x: NDArray, show: bool = False, **kwargs_solver) -> NDArray: x = self._step(x, **kwargs_solver) # compute residual - self.r: NDArray = self.y - self.Op * x + self.r: NDArray = self.y - self.Op.matvec(x) self.rnorm = self.ncp.linalg.norm(self.r) self.iiter += 1 @@ -537,7 +537,7 @@ def run( nouter = nouter if self.nouter is None else self.nouter if x is not None: self.x0 = x.copy() - self.y = self.y - self.Op * x + self.y = self.y - self.Op.matvec(x) # choose xold to ensure tolerance test is passed initially xold = x.copy() + np.inf while self.iiter < nouter and self.ncp.linalg.norm(x - xold) >= self.tolIRLS: @@ -1881,7 +1881,7 @@ def run( """ pinv, _, _, info = ext_spgl1( - self.Op if self.SOp is None else self.Op * self.SOp.H, + self.Op if self.SOp is None else self.Op @ self.SOp.H, self.y, tau=self.tau, sigma=self.sigma, @@ -2248,13 +2248,15 @@ def step( )[0] # Shrinkage self.d = [ - _softthreshold(self.RegsL1[ireg] * x + self.b[ireg], self.epsRL1s[ireg]) + _softthreshold( + self.RegsL1[ireg].matvec(x) + self.b[ireg], self.epsRL1s[ireg] + ) for ireg in range(self.nregsL1) ] # Bregman update self.b = [ - self.b[ireg] + self.tau * (self.RegsL1[ireg] * x - self.d[ireg]) + self.b[ireg] + self.tau * (self.RegsL1[ireg].matvec(x) - self.d[ireg]) for ireg in range(self.nregsL1) ] From 7b974510af8a26952a79615692a10ef84b22d891 Mon Sep 17 00:00:00 2001 From: cako Date: Sun, 7 May 2023 19:32:41 -0700 Subject: [PATCH 26/69] refactor: simplify forceflat logic --- pylops/linearoperator.py | 40 ++++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/pylops/linearoperator.py b/pylops/linearoperator.py index e397dfc3..67ffa59f 100644 --- a/pylops/linearoperator.py +++ b/pylops/linearoperator.py @@ -95,9 +95,10 @@ class LinearOperator(_LinearOperator): forceflat : :obj:`bool`, optional .. versionadded:: 2.2.0 - Force an array to be flattened after matvec/rmatvec if the input is ambiguous (i.e., is a 1D array both when - operating with ND arrays and with 1D arrays. Defaults to ``None`` for operators that have no ambiguity (and - to ``False`` for those with ambiguity) + Force an array to be flattened after matvec/rmatvec if the input is ambiguous + (i.e., is a 1D array both when operating with ND arrays and with 1D arrays). + Defaults to ``None`` for operators that have no ambiguity or to ``False`` + for those with ambiguity. name : :obj:`str`, optional .. versionadded:: 2.0.0 @@ -361,25 +362,26 @@ def __add__(self, x: LinearOperator) -> LinearOperator: ) Op.clinear = Op.clinear and Opx.clinear Op.explicit = False - # Define forceflat (if differing, raise error) - if isinstance(self.forceflat, bool) and isinstance(Opx.forceflat, bool): + if self.forceflat is None and Opx.forceflat is None: + Op.forceflat = None + elif self.forceflat is not None and Opx.forceflat is not None: + # Define forceflat only if differing, otherwise raise error if self.forceflat != Opx.forceflat: raise ValueError( - f"the two operators have contrasting forceflat {Op.forceflat}-{Opx.forceflat}" + f"Operators have conflicting forceflat {Op.forceflat} != {Opx.forceflat}" ) - else: - Op.forceflat = self.forceflat - if isinstance(self.forceflat, bool) or isinstance(Opx.forceflat, bool): + Op.forceflat = self.forceflat + else: # Only one of them is None Op.forceflat = ( self.forceflat if self.forceflat is not None else Opx.forceflat ) - else: - Op.forceflat = None + # Replace if shape-like if len(self.dims) == 1: Op.dims = Opx.dims if len(self.dimsd) == 1: Op.dimsd = Opx.dimsd + return Op else: return NotImplemented @@ -627,20 +629,19 @@ def dot(self, x: NDArray) -> NDArray: self._copy_attributes(Op, exclude=["dims", "explicit", "forceflat", "name"]) Op.clinear = Op.clinear and Opx.clinear Op.explicit = False - # Define forceflat (if differing, raise error) - if isinstance(self.forceflat, bool) and isinstance(Opx.forceflat, bool): + if self.forceflat is None and Opx.forceflat is None: + Op.forceflat = None + elif self.forceflat is not None and Opx.forceflat is not None: + # Define forceflat only if differing, otherwise raise error if self.forceflat != Opx.forceflat: raise ValueError( - f"the two operators have contrasting forceflat {Op.forceflat}-{Opx.forceflat}" + f"Operators have conflicting forceflat {Op.forceflat} != {Opx.forceflat}" ) - else: - Op.forceflat = self.forceflat - if isinstance(self.forceflat, bool) or isinstance(Opx.forceflat, bool): + Op.forceflat = self.forceflat + else: # Only one of them is None Op.forceflat = ( self.forceflat if self.forceflat is not None else Opx.forceflat ) - else: - Op.forceflat = None Op.dims = Opx.dims return Op elif np.isscalar(x): @@ -840,7 +841,6 @@ def tosparse(self) -> NDArray: # loop through columns of self for i in range(n): - # make a unit vector for the ith column unit_i = np.zeros(n) unit_i[i] = 1 From cbd42602b3d0188590abc6f37e0e5b98ba7e297c Mon Sep 17 00:00:00 2001 From: mrava87 Date: Mon, 8 May 2023 21:15:45 +0300 Subject: [PATCH 27/69] feat: forceflat test for Bilinear --- pytests/test_basicoperators.py | 21 +++++++++++---------- pytests/test_interpolation.py | 24 ++++++++++++++++++++++++ 2 files changed, 35 insertions(+), 10 deletions(-) diff --git a/pytests/test_basicoperators.py b/pytests/test_basicoperators.py index 43cb15f2..8f3f0528 100755 --- a/pytests/test_basicoperators.py +++ b/pytests/test_basicoperators.py @@ -468,8 +468,19 @@ def test_Roll3D(par): assert_array_almost_equal(x[str(axis)].ravel(), xinv, decimal=3) +@pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j), (par3)]) +def test_Sum2D(par): + """Dot-test for Sum operator on 2d signal""" + for axis in [0, 1]: + dim_d = [par["ny"], par["nx"]] + dim_d.pop(axis) + Sop = Sum(dims=(par["ny"], par["nx"]), axis=axis, dtype=par["dtype"]) + assert dottest(Sop, np.prod(dim_d), par["ny"] * par["nx"]) + + @pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j), (par3)]) def test_Sum2D_forceflat(par): + """Dot-test for Sum operator on 2d signal with forceflat""" np.random.seed(10) flat_dimsd = par["ny"] flat_dims = par["ny"] * par["nx"] @@ -506,16 +517,6 @@ def test_Sum2D_forceflat(par): assert Sop.forceflat is None -@pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j), (par3)]) -def test_Sum2D(par): - """Dot-test for Sum operator on 2d signal""" - for axis in [0, 1]: - dim_d = [par["ny"], par["nx"]] - dim_d.pop(axis) - Sop = Sum(dims=(par["ny"], par["nx"]), axis=axis, dtype=par["dtype"]) - assert dottest(Sop, np.prod(dim_d), par["ny"] * par["nx"]) - - @pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j), (par3)]) def test_Sum3D(par): """Dot-test, forward and adjoint for Sum operator on 3d signal""" diff --git a/pytests/test_interpolation.py b/pytests/test_interpolation.py index a691ba16..14473ba0 100755 --- a/pytests/test_interpolation.py +++ b/pytests/test_interpolation.py @@ -405,6 +405,30 @@ def test_Bilinear_2dsignal(par): assert_array_almost_equal(y, x[iava[0], iava[1]]) +@pytest.mark.parametrize("par", [(par1), (par2)]) +def test_Bilinear_2dsignal_flatten(par): + """Dot-test and forward for Interp operator for 2d signal with forceflat""" + np.random.seed(1) + flat_dimsd = par["ny"] + flat_dims = par["nx"] * par["nt"] + dims = (par["nx"], par["nt"]) + + x = np.random.normal(0, 1, dims) + par["imag"] * np.random.normal(0, 1, dims) + + iava = np.vstack((np.arange(0, 10), np.arange(0, 10))) + Iop_True = Bilinear(iava, dims=dims, dtype=par["dtype"], forceflat=True) + y = Iop_True @ x + xadj = Iop_True.H @ y + assert y.shape == (flat_dimsd,) + assert xadj.shape == (flat_dims,) + + Iop_None = Bilinear(iava, dims=dims, dtype=par["dtype"]) + y = Iop_None @ x + xadj = Iop_None.H @ y + assert y.shape == (par["ny"],) + assert xadj.shape == dims + + @pytest.mark.parametrize("par", [(par1), (par2)]) def test_Bilinear_3dsignal(par): """Dot-test and forward for Interp operator for 3d signal""" From 01bfce768a08f935e555482531acc01e28ed6444 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Mon, 8 May 2023 21:33:02 +0300 Subject: [PATCH 28/69] fix: fix forceflat test for Bilinear --- pytests/test_interpolation.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/pytests/test_interpolation.py b/pytests/test_interpolation.py index 14473ba0..a22c04ef 100755 --- a/pytests/test_interpolation.py +++ b/pytests/test_interpolation.py @@ -409,23 +409,23 @@ def test_Bilinear_2dsignal(par): def test_Bilinear_2dsignal_flatten(par): """Dot-test and forward for Interp operator for 2d signal with forceflat""" np.random.seed(1) - flat_dimsd = par["ny"] - flat_dims = par["nx"] * par["nt"] dims = (par["nx"], par["nt"]) + flat_dims = par["nx"] * par["nt"] + dimsd = 10 x = np.random.normal(0, 1, dims) + par["imag"] * np.random.normal(0, 1, dims) - iava = np.vstack((np.arange(0, 10), np.arange(0, 10))) + iava = np.vstack((np.arange(0, dimsd), np.arange(0, dimsd))) Iop_True = Bilinear(iava, dims=dims, dtype=par["dtype"], forceflat=True) y = Iop_True @ x xadj = Iop_True.H @ y - assert y.shape == (flat_dimsd,) + assert y.shape == (dimsd,) assert xadj.shape == (flat_dims,) Iop_None = Bilinear(iava, dims=dims, dtype=par["dtype"]) y = Iop_None @ x xadj = Iop_None.H @ y - assert y.shape == (par["ny"],) + assert y.shape == (dimsd,) assert xadj.shape == dims From f2c919437d714c0e7eecb95b2b085c8eea613afd Mon Sep 17 00:00:00 2001 From: Yuxi Hong Date: Tue, 16 May 2023 16:00:58 +0300 Subject: [PATCH 29/69] minor: fix a typo in Torch description in docs --- docs/source/installation.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/installation.rst b/docs/source/installation.rst index 2e22eb2d..3b237f23 100755 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -487,7 +487,7 @@ or via ``pip`` with Torch ----- -`Torch `_ used to allow seamless integration between PyLops and PyTorch operators. +`Torch `_ is used to allow seamless integration between PyLops and PyTorch operators. Install it via ``conda`` with: From 30b57644f18c058724be10c9f38ce7086901daf7 Mon Sep 17 00:00:00 2001 From: Wei Zhang Date: Mon, 22 May 2023 15:25:02 +0300 Subject: [PATCH 30/69] feature: added nonstatconvolve3d operator --- .../_nonstatconvolve3d_cuda.py | 153 +++++++++ pylops/signalprocessing/nonstatconvolve3d.py | 310 ++++++++++++++++++ 2 files changed, 463 insertions(+) create mode 100644 pylops/signalprocessing/_nonstatconvolve3d_cuda.py create mode 100644 pylops/signalprocessing/nonstatconvolve3d.py diff --git a/pylops/signalprocessing/_nonstatconvolve3d_cuda.py b/pylops/signalprocessing/_nonstatconvolve3d_cuda.py new file mode 100644 index 00000000..8602d272 --- /dev/null +++ b/pylops/signalprocessing/_nonstatconvolve3d_cuda.py @@ -0,0 +1,153 @@ +from math import floor + +from numba import cuda + + +@cuda.jit(max_registers=40) +def _matvec_rmatvec(x, y, hs, hshape, xdims, ohx, ohy, ohz, dhx, dhy, dhz, nhx, nhy, nhz, rmatvec): + """Cuda kernels for NonStationaryConvolve3D operator + + Cuda implementation of matvec and rmatvec for NonStationaryConvolve3D operator. See + :class:`pylops.signalprocessing.NonStationaryConvolve3D` for details about input parameters. + + """ + ix, iy, iz = cuda.grid(3) + + if ix < xdims[0] and iy < xdims[1] and iz < xdims[2]: + + ## find closest filters and interpolate h + ihx_l = int(floor((ix - ohx) / dhx)) #id number of left for hs_arr + ihy_b = int(floor((iy - ohy) / dhy)) #id number of back for hs_arr + ihz_t = int(floor((iz - ohz) / dhz)) #id number of top for hs_arr + + dhx_r = (ix - ohx) / dhx - ihx_l #weight for right psfs, left 1-dhx_r + dhy_f = (iy - ohy) / dhy - ihy_b #weight for front psfs, back 1-dhy_f + dhz_d = (iz - ohz) / dhz - ihz_t #weight for down psfs, top 1-dhz_d + + if ihx_l < 0: + ihx_l = ihx_r = 0 + dhx_l = dhx_r = 0.5 + elif ihx_l >= nhx - 1: + ihx_l = ihx_r = nhx - 1 + dhx_l = dhx_r = 0.5 + else: + ihx_r = ihx_l + 1 + dhx_l = 1.0 - dhx_r + + if ihy_b < 0: + ihy_b = ihy_f = 0 + dhy_b = dhy_f = 0.5 + elif ihy_b >= nhy - 1: + ihy_b = ihy_f = nhy - 1 + dhy_b = dhy_f = 0.5 + else: + ihy_f = ihy_b + 1 + dhy_b = 1.0 - dhy_f + + if ihz_t < 0: + ihz_t = ihz_d = 0 + dhz_t = dhz_d = 0.5 + elif ihz_t >= nhz - 1: + ihz_t = ihz_d = nhz - 1 + dhz_t = dhz_d = 0.5 + else: + ihz_d = ihz_t + 1 + dhz_t = 1.0 - dhz_d + + h_lbt = hs[ihx_l, ihy_b, ihz_t] + h_lbd = hs[ihx_l, ihy_b, ihz_d] + h_lft = hs[ihx_l, ihy_f, ihz_t] + h_lfd = hs[ihx_l, ihy_f, ihz_d] + + h_rbt = hs[ihx_r, ihy_b, ihz_t] + h_rbd = hs[ihx_r, ihy_b, ihz_d] + h_rft = hs[ihx_r, ihy_f, ihz_t] + h_rfd = hs[ihx_r, ihy_f, ihz_d] + + ## find extremes of model where to apply h (in case h is going out of model) + xextremes = ( + max(0, ix - hshape[0] // 2), + min(ix + hshape[0] // 2 + 1, xdims[0]), + ) + yextremes = ( + max(0, iy - hshape[1] // 2), + min(iy + hshape[1] // 2 + 1, xdims[1]), + ) + zextremes = ( + max(0, iz - hshape[2] // 2), + min(iz + hshape[2] // 2 + 1, xdims[2]), + ) + + + ## find extremes of h (in case h is going out of model) + hxextremes = ( + max(0, -ix + hshape[0] // 2), + min(hshape[0], hshape[0] // 2 + (xdims[0] - ix)), + ) + hyextremes = ( + max(0, -iy + hshape[1] // 2), + min(hshape[1], hshape[1] // 2 + (xdims[1] - iy)), + ) + hzextremes = ( + max(0, -iz + hshape[2] // 2), + min(hshape[2], hshape[2] // 2 + (xdims[2] - iz)), + ) + + # place filter in output + for ixx, hxx in zip( + range(xextremes[0], xextremes[1]), range(hxextremes[0], hxextremes[1]) + ): + for iyy, hyy in zip( + range(yextremes[0], yextremes[1]), range(hyextremes[0], hyextremes[1]) + ): + for izz, hzz in zip( + range(zextremes[0], zextremes[1]), + range(hzextremes[0], hzextremes[1]), + ): + h = ( + dhx_l * dhy_b * dhz_t * h_lbt[hxx, hyy, hzz] + + dhx_l * dhy_b * dhz_d * h_lbd[hxx, hyy, hzz] + + dhx_l * dhy_f * dhz_t * h_lft[hxx, hyy, hzz] + + dhx_l * dhy_f * dhz_d * h_lfd[hxx, hyy, hzz] + + dhx_r * dhy_b * dhz_t * h_rbt[hxx, hyy, hzz] + + dhx_r * dhy_b * dhz_d * h_rbd[hxx, hyy, hzz] + + dhx_r * dhy_f * dhz_t * h_rft[hxx, hyy, hzz] + + dhx_r * dhy_f * dhz_d * h_rfd[hxx, hyy, hzz] + ) + + if rmatvec: + y[ix, iy, iz] += h * x[ixx, iyy, izz] + else: + cuda.atomic.add(y, (ixx, iyy, izz), x[ix, iy, iz] * h) + + +def _matvec_rmatvec_call( + x, + y, + hs, + hshape, + xdims, + ohx, + ohy, + ohz, + dhx, + dhy, + dhz, + nhx, + nhy, + nhz, + rmatvec= False, + dim_block=(1,32,16), + dim_grid=(24,24,24), +): + """Caller for NonStationaryConvolve3D operator + + Caller for cuda implementation of matvec and rmatvec for NonStationaryConvolve3D operator, with same signature + as numpy/numba counterparts. See :class:`pylops.signalprocessing.NonStationaryConvolve3D` for details about + input parameters. + + """ + _matvec_rmatvec[dim_grid, dim_block]( + x, y, hs, hshape, xdims, ohx, ohy, ohz, dhx, dhy, dhz, nhx, nhy, nhz, rmatvec + ) + return y diff --git a/pylops/signalprocessing/nonstatconvolve3d.py b/pylops/signalprocessing/nonstatconvolve3d.py new file mode 100644 index 00000000..61628d9d --- /dev/null +++ b/pylops/signalprocessing/nonstatconvolve3d.py @@ -0,0 +1,310 @@ +__all__ = ["NonStationaryConvolve3D"] + +import os +from typing import Tuple, Union + +import numpy as np + +from pylops import LinearOperator +from pylops.utils import deps +from pylops.utils._internal import _value_or_sized_to_tuple +from pylops.utils.backend import get_array_module +from pylops.utils.decorators import reshaped +from pylops.utils.typing import DTypeLike, InputDimsLike, NDArray + +jit_message = deps.numba_import("the nonstatconvolve3d module") +# jit_message = None +if jit_message is None: + from numba import jit, prange + + from _nonstatconvolve3d_cuda import _matvec_rmatvec_call as _Convolve3D_kernel + + # detect whether to use parallel or not + numba_threads = int(os.getenv("NUMBA_NUM_THREADS", "1")) + parallel = True if numba_threads != 1 else False +else: + prange = range + + + + +class NonStationaryConvolve3D(LinearOperator): + r"""3D non-stationary convolution operator. + + Apply non-stationary three-dimensional convolution. A varying compact filter + is provided on a coarser grid and on-the-fly interpolation is applied + in forward and adjoint modes. + + Parameters + ---------- + dims : :obj:`list` or :obj:`int` + Number of samples for each dimension + hs : :obj:`numpy.ndarray` + Bank of 3d compact filters of size + :math:`n_{\text{filts},x} \times n_{\text{filts},y} \times n_{\text{filts},z} \times n_h \times n_{h,x} \times n_{h,y} \times n_{h,z}`. + Filters must have odd number of samples and are assumed to be + centered in the middle of the filter support. + ihx : :obj:`tuple` + Indices of the x locations of the filters ``hs`` in the model (and data). Note + that the filters must be regularly sampled, i.e. :math:`dh_x=\text{diff}(ihx)=\text{const.}` + ihy : :obj:`tuple` + Indices of the y locations of the filters ``hs`` in the model (and data). Note + that the filters must be regularly sampled, i.e. :math:`dh_y=\text{diff}(ihy)=\text{const.}` + ihz : :obj:`tuple` + Indices of the z locations of the filters ``hs`` in the model (and data). Note + that the filters must be regularly sampled, i.e. :math:`dh_z=\text{diff}(ihz)=\text{const.}` + engine : :obj:`str`, optional + Engine used for spread computation (``numpy``, ``numba``, or ``cuda``) + dim_block : :obj:`tuple`, optional + Number of threads in each block (only when ``engine=cuda``) + dtype : :obj:`str`, optional + Type of elements in input array. + name : :obj:`str`, optional + Name of operator (to be used by :func:`pylops.utils.describe.describe`) + + Attributes + ---------- + shape : :obj:`tuple` + Operator shape + explicit : :obj:`bool` + Operator contains a matrix that can be solved explicitly (``True``) or + not (``False``) + + Raises + ------ + ValueError + If filters ``hs`` have even size + ValueError + If ``ihx``, ``ihy`` or ``ihz`` is not regularly sampled + NotImplementedError + If ``engine`` is neither ``numpy``, ``fftw``, nor ``scipy``. + """ + def __init__( + self, + dims: Union[int, InputDimsLike], + hs: NDArray, + ihx: InputDimsLike, + ihy: InputDimsLike, + ihz: InputDimsLike, + engine: str = "numpy", + dim_block: Tuple[int, int] = (2, 16, 16), + dtype: DTypeLike = "float64", + name: str = "C") -> None: + if engine not in ["numpy", "numba", "cuda"]: + raise NotImplementedError("engine must be numpy or numba or cuda") + if hs.shape[3] % 2 == 0 or hs.shape[4] % 2 == 0 or hs.shape[5] % 2 == 0: + raise ValueError("filters hs must have odd length") + #init + self.hs = hs + self.hshape = hs.shape[3:] + self.ohx, self.dhx, self.nhx = ihx[0], ihx[1] - ihx[0], len(ihx) + self.ohy, self.dhy, self.nhy = ihy[0], ihy[1] - ihy[0], len(ihy) + self.ohz, self.dhz, self.nhz = ihz[0], ihz[1] - ihz[0], len(ihz) + self.ehx, self.ehx, self.ehz = ihx[-1], ihy[-1], ihz[-1] + # self.dims = tuple(dims) + self.dims = _value_or_sized_to_tuple(dims) + self.engine = engine + # self.dtype = dtype + # self.shape = dims + + super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dims, name=name) + + + # create additional input parameters for engine=cuda + self.kwargs_cuda = {} + if engine == "cuda": + self.kwargs_cuda["dim_block"] = dim_block + + gridx = (self.dims[0] + dim_block[0] - 1)//dim_block[0] + gridy = (self.dims[1] + dim_block[1] - 1)//dim_block[1] + gridz = (self.dims[2] + dim_block[2] - 1)//dim_block[2] + + self.kwargs_cuda["dim_grid"] = (gridx, gridy, gridz) + + self._register_multiplications(engine) + + def _register_multiplications(self, engine: str) -> None: + if engine == "numba": + numba_opts = dict(nopython=True, fastmath=True, nogil=True, parallel=True) + self._mvrmv = jit(**numba_opts)(self._matvec_rmatvec) + elif engine == "cuda": + self._mvrmv = _Convolve3D_kernel + else: + self._mvrmv = self._matvec_rmatvec + + @staticmethod + def _matvec_rmatvec( + x: NDArray, + y: NDArray, + hs: NDArray, + hshape: Tuple[int, int, int], + xdims: Tuple[int, int, int], + ohx: float, + ohy: float, + ohz: float, + dhx: float, + dhy: float, + dhz: float, + nhx: int, + nhy: int, + nhz: int, + rmatvec: bool = False, + ) -> NDArray: + for ix in prange(xdims[0]): + for iy in range(xdims[1]): + for iz in range(xdims[2]): + # find closest filters and interpolate h + ihx_l = int(np.floor((ix - ohx) / dhx)) #id number of left for hs_arr + ihy_b = int(np.floor((iy - ohy) / dhy)) #id number of back for hs_arr + ihz_t = int(np.floor((iz - ohz) / dhz)) #id number of top for hs_arr + + dhx_r = (ix - ohx) / dhx - ihx_l #weight for right psfs, left 1-ihz_t + dhy_f = (iy - ohy) / dhy - ihy_b #weight for front psfs, left 1-ihz_t + dhz_d = (iz - ohz) / dhz - ihz_t #weight for down psfs, top 1-dhz_d + + if ihx_l < 0: + ihx_l = ihx_r = 0 + dhx_l = dhx_r = 0.5 + elif ihx_l >= nhx - 1: + ihx_l = ihx_r = nhx - 1 + dhx_l = dhx_r = 0.5 + else: + ihx_r = ihx_l + 1 + dhx_l = 1.0 - dhx_r + + if ihy_b < 0: + ihy_b = ihy_f = 0 + dhy_b = dhy_f = 0.5 + elif ihy_b >= nhy - 1: + ihy_b = ihy_f = nhy - 1 + dhy_b = dhy_f = 0.5 + else: + ihy_f = ihy_b + 1 + dhy_b = 1.0 - dhy_f + + if ihz_t < 0: + ihz_t = ihz_d = 0 + dhz_t = dhz_d = 0.5 + elif ihz_t >= nhz - 1: + ihz_t = ihz_d = nhz - 1 + dhz_t = dhz_d = 0.5 + else: + ihz_d = ihz_t + 1 + dhz_t = 1.0 - dhz_d + + h_lbt = hs[ihx_l, ihy_b, ihz_t] + h_lbd = hs[ihx_l, ihy_b, ihz_d] + h_lft = hs[ihx_l, ihy_f, ihz_t] + h_lfd = hs[ihx_l, ihy_f, ihz_d] + + h_rbt = hs[ihx_r, ihy_b, ihz_t] + h_rbd = hs[ihx_r, ihy_b, ihz_d] + h_rft = hs[ihx_r, ihy_f, ihz_t] + h_rfd = hs[ihx_r, ihy_f, ihz_d] + + h = ( + dhx_l * dhy_b * dhz_t * h_lbt + + dhx_l * dhy_b * dhz_d * h_lbd + + dhx_l * dhy_f * dhz_t * h_lft + + dhx_l * dhy_f * dhz_d * h_lfd + + dhx_r * dhy_b * dhz_t * h_rbt + + dhx_r * dhy_b * dhz_d * h_rbd + + dhx_r * dhy_f * dhz_t * h_rft + + dhx_r * dhy_f * dhz_d * h_rfd + ) + + # find extremes of model where to apply h (in case h is going out of model) + xextremes = ( + max(0, ix - hshape[0] // 2), + min(ix + hshape[0] // 2 + 1, xdims[0]), + ) + yextremes = ( + max(0, iy - hshape[1] // 2), + min(iy + hshape[1] // 2 + 1, xdims[1]), + ) + zextremes = ( + max(0, iz - hshape[2] // 2), + min(iz + hshape[2] // 2 + 1, xdims[2]), + ) + + + # find extremes of h (in case h is going out of model) + hxextremes = ( + max(0, -ix + hshape[0] // 2), + min(hshape[0], hshape[0] // 2 + (xdims[0] - ix)), + ) + hyextremes = ( + max(0, -iy + hshape[1] // 2), + min(hshape[1], hshape[1] // 2 + (xdims[1] - iy)), + ) + hzextremes = ( + max(0, -iz + hshape[2] // 2), + min(hshape[2], hshape[2] // 2 + (xdims[2] - iz)), + ) + + if not rmatvec: + y[ xextremes[0] : xextremes[1], yextremes[0] : yextremes[1], zextremes[0] : zextremes[1] ] += ( + x[ix, iy, iz] + * h[ hxextremes[0] : hxextremes[1], hyextremes[0] : hyextremes[1], hzextremes[0] : hzextremes[1] ] + ) + else: + y[ix, iy, iz] = np.sum( + h[ hxextremes[0] : hxextremes[1], hyextremes[0] : hyextremes[1], hzextremes[0] : hzextremes[1] ] + * x[ xextremes[0] : xextremes[1], yextremes[0] : yextremes[1], zextremes[0] : zextremes[1] ] + ) + return y + + # def _forward(self, x): + @reshaped #reshape to 1D array + def _matvec(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + # ncp = cp.get_array_module(x) + # y = ncp.zeros(self.dims, dtype=x.dtype) + y = ncp.zeros(self.dims, dtype=self.dtype) + y = self._mvrmv( + x, + y, + self.hs, + self.hshape, + self.dims, + self.ohx, + self.ohy, + self.ohz, + self.dhx, + self.dhy, + self.dhz, + self.nhx, + self.nhy, + self.nhz, + rmatvec=False, + **self.kwargs_cuda + ) + return y + + + # def _adjoint(self, x): + @reshaped #reshape to 1D array + def _rmatvec(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + # ncp = cp.get_array_module(x) + # y = ncp.zeros(self.dims, dtype=x.dtype) + y = ncp.zeros(self.dims, dtype=self.dtype) + y = self._mvrmv( + x, + y, + self.hs, + self.hshape, + self.dims, + self.ohx, + self.ohy, + self.ohz, + self.dhx, + self.dhy, + self.dhz, + self.nhx, + self.nhy, + self.nhz, + rmatvec=True, + **self.kwargs_cuda + ) + return y From 3f65b0247de0a1028aed41aa867e178a207fc034 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Fri, 26 May 2023 21:04:48 +0300 Subject: [PATCH 31/69] bug: fix size of wavelet in plot_ista to be smaller than that of input --- examples/plot_ista.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/plot_ista.py b/examples/plot_ista.py index e37f2247..d6359905 100755 --- a/examples/plot_ista.py +++ b/examples/plot_ista.py @@ -106,7 +106,7 @@ x[int(nt / 2)] = 1 x[nt - 20] = 0.5 -h, th, hcenter = pylops.utils.wavelets.ricker(t[:101], f0=20) +h, th, hcenter = pylops.utils.wavelets.ricker(t[:21], f0=20) Cop = pylops.signalprocessing.Convolve1D(nt, h=h, offset=hcenter, dtype="float32") y = Cop * x From 7c34a7b7a2dd22ea3a65797f338d1765cf67e394 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Thu, 1 Jun 2023 14:17:45 +0300 Subject: [PATCH 32/69] patch: allow psnr to use range instead of max --- pylops/utils/metrics.py | 15 +++++++++++---- pytests/test_metrics.py | 4 ++-- 2 files changed, 13 insertions(+), 6 deletions(-) diff --git a/pylops/utils/metrics.py b/pylops/utils/metrics.py index fa71967b..e1e7a55c 100644 --- a/pylops/utils/metrics.py +++ b/pylops/utils/metrics.py @@ -75,10 +75,10 @@ def snr(xref, xcmp): return snr -def psnr(xref, xcmp, xmax=None): +def psnr(xref, xcmp, xmax=None, xmin=0.0): """Peak Signal to Noise Ratio (PSNR) - Compute Peak Signal to Noise Ratio between two vectors. + Compute Peak Signal to Noise Ratio between two vectors Parameters ---------- @@ -89,6 +89,10 @@ def psnr(xref, xcmp, xmax=None): xmax : :obj:`float`, optional Maximum value to use. If ``None``, the actual maximum of the reference vector is used + xmin : :obj:`float`, optional + Minimum value to use. If ``None``, the actual minimum of + the reference vector is used (``0`` is default for + backward compatibility) Returns ------- @@ -98,5 +102,8 @@ def psnr(xref, xcmp, xmax=None): """ if xmax is None: xmax = xref.max() - psrn = 10.0 * np.log10(xmax**2 / mse(xref, xcmp)) - return psrn + if xmin is None: + xmin = xref.min() + xrange = xmax - xmin + psnr = 10.0 * np.log10(xrange**2 / mse(xref, xcmp)) + return psnr diff --git a/pytests/test_metrics.py b/pytests/test_metrics.py index 7ca20f67..420434af 100755 --- a/pytests/test_metrics.py +++ b/pytests/test_metrics.py @@ -49,7 +49,7 @@ def test_psnr(par): xref = np.ones(par["nx"]) xcmp = np.zeros(par["nx"]) - psnrsame = psnr(xref, xref, xmax=1.0) - psnrcmp = psnr(xref, xcmp, xmax=1.0) + psnrsame = psnr(xref, xref, xmax=1.0, xmin=0.0) + psnrcmp = psnr(xref, xcmp, xmax=1.0, xmin=0.0) assert psnrsame == np.inf assert psnrcmp == 0.0 From 4d0399e1e9ecb9177da1e7671510d15e3efca8bb Mon Sep 17 00:00:00 2001 From: mrava87 Date: Fri, 2 Jun 2023 21:29:28 +0300 Subject: [PATCH 33/69] minor: simplify if else statements for offset/pad --- pylops/signalprocessing/convolve1d.py | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/pylops/signalprocessing/convolve1d.py b/pylops/signalprocessing/convolve1d.py index 8208d0ce..9093f526 100644 --- a/pylops/signalprocessing/convolve1d.py +++ b/pylops/signalprocessing/convolve1d.py @@ -80,10 +80,7 @@ def __init__( if self.offset != 0: self.h = _pad_along_axis( self.h, - ( - self.offset if self.offset > 0 else 0, - -self.offset if self.offset < 0 else 0, - ), + (max(self.offset, 0), -min(self.offset, 0)), axis=-1 if h.ndim == 1 else axis, ) self.hstar = np.flip(self.h, axis=-1) @@ -146,12 +143,12 @@ def __init__( self.hstar = np.flip(self.h, axis=-1) self.pad = np.zeros((len(dims), 2), dtype=int) - self.pad[self.axis, 0] = self.offset if self.offset > 0 else 0 - self.pad[self.axis, 1] = -self.offset if self.offset < 0 else 0 + self.pad[self.axis, 0] = max(self.offset, 0) + self.pad[self.axis, 1] = -min(self.offset, 0) self.padd = np.zeros((len(dims), 2), dtype=int) - self.padd[self.axis, 1] = self.offset if self.offset > 0 else 0 - self.padd[self.axis, 0] = -self.offset if self.offset < 0 else 0 + self.padd[self.axis, 1] = max(self.offset, 0) + self.padd[self.axis, 0] = -min(self.offset, 0) # add dimensions to filter to match dimensions of model and data if self.h.ndim == 1: From b5eefc707ec7083b970ac90cd23abbef68c2fd4f Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sun, 2 Jul 2023 13:21:25 +0300 Subject: [PATCH 34/69] bug: change jit into jit_message in checks --- pylops/waveeqprocessing/kirchhoff.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pylops/waveeqprocessing/kirchhoff.py b/pylops/waveeqprocessing/kirchhoff.py index 595664df..c34db994 100644 --- a/pylops/waveeqprocessing/kirchhoff.py +++ b/pylops/waveeqprocessing/kirchhoff.py @@ -1221,7 +1221,7 @@ def _ampsrcrec_kirch_rmatvec( def _register_multiplications(self, engine: str) -> None: if engine not in ["numpy", "numba"]: raise KeyError("engine must be numpy or numba") - if engine == "numba" and jit is not None: + if engine == "numba" and jit_message is None: # numba numba_opts = dict( nopython=True, nogil=True, parallel=parallel @@ -1240,7 +1240,7 @@ def _register_multiplications(self, engine: str) -> None: self._kirch_rmatvec = jit(**numba_opts)(self._trav_kirch_rmatvec) else: - if engine == "numba" and jit is None: + if engine == "numba" and jit_message is not None: logging.warning(jit_message) if self.dynamic and self.travsrcrec: self._kirch_matvec = self._ampsrcrec_kirch_matvec From 8aba0a60f9db464c7e31c17cc7900e083ac1e846 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Thu, 6 Jul 2023 20:34:55 +0300 Subject: [PATCH 35/69] bug: fixed skfmm check in kirchhoff --- pylops/waveeqprocessing/kirchhoff.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pylops/waveeqprocessing/kirchhoff.py b/pylops/waveeqprocessing/kirchhoff.py index c34db994..a743ec6d 100644 --- a/pylops/waveeqprocessing/kirchhoff.py +++ b/pylops/waveeqprocessing/kirchhoff.py @@ -578,7 +578,7 @@ def _traveltime_table( dist_recs = trav_recs * vel elif mode == "eikonal": - if skfmm is not None: + if skfmm_message is None: dist_srcs = np.zeros((ny * nx * nz, ns)) dist_recs = np.zeros((ny * nx * nz, nr)) trav_srcs = np.zeros((ny * nx * nz, ns)) From 5d2f600cd355fa2bb13db4562fd2991a73d3d705 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Wed, 19 Jul 2023 16:56:58 +0300 Subject: [PATCH 36/69] minor: added tests for NonStationaryConvolve3D and finalized cleaning of files --- README.md | 1 + docs/source/api/index.rst | 1 + docs/source/credits.rst | 3 +- pylops/signalprocessing/__init__.py | 3 + .../_nonstatconvolve3d_cuda.py | 6 +- pylops/signalprocessing/nonstatconvolve2d.py | 6 +- pylops/signalprocessing/nonstatconvolve3d.py | 97 +++++++++---------- pytests/test_nonstatconvolve.py | 74 ++++++++++++++ 8 files changed, 134 insertions(+), 57 deletions(-) diff --git a/README.md b/README.md index a9258eef..a39436be 100644 --- a/README.md +++ b/README.md @@ -147,3 +147,4 @@ A list of video tutorials to learn more about PyLops: * Juan Daniel Romero, jdromerom * Aniket Singh Rawat, dikwickley * Rohan Babbar, rohanbabbar04 +* Wei Zhang, ZhangWeiGeo diff --git a/docs/source/api/index.rst b/docs/source/api/index.rst index ea25d5bf..9d7c3f66 100755 --- a/docs/source/api/index.rst +++ b/docs/source/api/index.rst @@ -91,6 +91,7 @@ Signal processing ConvolveND NonStationaryConvolve1D NonStationaryConvolve2D + NonStationaryConvolve3D NonStationaryFilters1D NonStationaryFilters2D Interp diff --git a/docs/source/credits.rst b/docs/source/credits.rst index 2a8558d1..036554b8 100755 --- a/docs/source/credits.rst +++ b/docs/source/credits.rst @@ -18,4 +18,5 @@ Contributors * `Muhammad Izzatullah `_, izzatum * `Juan Daniel Romero `_, jdromerom * `Aniket Singh Rawat `_, dikwickley -* `Rohan Babbar `_, rohanbabbar04 \ No newline at end of file +* `Rohan Babbar `_, rohanbabbar04 +* `Wei Zhang `_, ZhangWeiGeo \ No newline at end of file diff --git a/pylops/signalprocessing/__init__.py b/pylops/signalprocessing/__init__.py index 2ce1fed1..8efa532e 100755 --- a/pylops/signalprocessing/__init__.py +++ b/pylops/signalprocessing/__init__.py @@ -12,6 +12,7 @@ ConvolveND ND convolution operator. NonStationaryConvolve1D 1D nonstationary convolution operator. NonStationaryConvolve2D 2D nonstationary convolution operator. + NonStationaryConvolve3D 3D nonstationary convolution operator. NonStationaryFilters1D 1D nonstationary filter estimation operator. NonStationaryFilters2D 2D nonstationary filter estimation operator. Interp Interpolation operator. @@ -43,6 +44,7 @@ from .convolve2d import * from .nonstatconvolve1d import * from .nonstatconvolve2d import * +from .nonstatconvolve3d import * from .shift import * from .interp import * from .bilinear import * @@ -71,6 +73,7 @@ "Convolve2D", "NonStationaryConvolve1D", "NonStationaryConvolve2D", + "NonStationaryConvolve3D", "NonStationaryFilters1D", "NonStationaryFilters2D", "Interp", diff --git a/pylops/signalprocessing/_nonstatconvolve3d_cuda.py b/pylops/signalprocessing/_nonstatconvolve3d_cuda.py index 8602d272..a5810f03 100644 --- a/pylops/signalprocessing/_nonstatconvolve3d_cuda.py +++ b/pylops/signalprocessing/_nonstatconvolve3d_cuda.py @@ -137,8 +137,8 @@ def _matvec_rmatvec_call( nhy, nhz, rmatvec= False, - dim_block=(1,32,16), - dim_grid=(24,24,24), + num_blocks=(1,32,16), + num_threads_per_blocks=(24,24,24), ): """Caller for NonStationaryConvolve3D operator @@ -147,7 +147,7 @@ def _matvec_rmatvec_call( input parameters. """ - _matvec_rmatvec[dim_grid, dim_block]( + _matvec_rmatvec[num_blocks, num_threads_per_blocks]( x, y, hs, hshape, xdims, ohx, ohy, ohz, dhx, dhy, dhz, nhx, nhy, nhz, rmatvec ) return y diff --git a/pylops/signalprocessing/nonstatconvolve2d.py b/pylops/signalprocessing/nonstatconvolve2d.py index 569736d6..8cd025c1 100644 --- a/pylops/signalprocessing/nonstatconvolve2d.py +++ b/pylops/signalprocessing/nonstatconvolve2d.py @@ -36,15 +36,15 @@ class NonStationaryConvolve2D(LinearOperator): Apply non-stationary two-dimensional convolution. A varying compact filter is provided on a coarser grid and on-the-fly interpolation is applied in forward and adjoint modes. - Both input and output have size :math`n_x \time n_z`. + Both input and output have size :math`n_x \times n_z`. Parameters ---------- dims : :obj:`list` or :obj:`int` - Number of samples for each dimension (which we refer to as :math`n_x \time n_z`). + Number of samples for each dimension (which we refer to as :math`n_x \times n_z`). hs : :obj:`numpy.ndarray` Bank of 2d compact filters of size - :math:`n_{\text{filts},x} \times n_{\text{filts},z} \times n_h \times n_{h,x} \times n_{h,z}`. + :math:`n_{\text{filts},x} \times n_{\text{filts},z} \times n_{h,x} \times n_{h,z}`. Filters must have odd number of samples and are assumed to be centered in the middle of the filter support. ihx : :obj:`tuple` diff --git a/pylops/signalprocessing/nonstatconvolve3d.py b/pylops/signalprocessing/nonstatconvolve3d.py index 61628d9d..1b7eaf3b 100644 --- a/pylops/signalprocessing/nonstatconvolve3d.py +++ b/pylops/signalprocessing/nonstatconvolve3d.py @@ -13,11 +13,13 @@ from pylops.utils.typing import DTypeLike, InputDimsLike, NDArray jit_message = deps.numba_import("the nonstatconvolve3d module") -# jit_message = None + if jit_message is None: from numba import jit, prange - from _nonstatconvolve3d_cuda import _matvec_rmatvec_call as _Convolve3D_kernel + from ._nonstatconvolve2d_cuda import ( + _matvec_rmatvec_call as _matvec_rmatvec_cuda_call, + ) # detect whether to use parallel or not numba_threads = int(os.getenv("NUMBA_NUM_THREADS", "1")) @@ -26,22 +28,20 @@ prange = range - - class NonStationaryConvolve3D(LinearOperator): r"""3D non-stationary convolution operator. Apply non-stationary three-dimensional convolution. A varying compact filter is provided on a coarser grid and on-the-fly interpolation is applied - in forward and adjoint modes. + in forward and adjoint modes. Both input and output have size :math`n_x \times n_y \times n_z`. Parameters ---------- dims : :obj:`list` or :obj:`int` - Number of samples for each dimension + Number of samples for each dimension (which we refer to as :math`n_x \times n_y \times n_z`). hs : :obj:`numpy.ndarray` Bank of 3d compact filters of size - :math:`n_{\text{filts},x} \times n_{\text{filts},y} \times n_{\text{filts},z} \times n_h \times n_{h,x} \times n_{h,y} \times n_{h,z}`. + :math:`n_{\text{filts},x} \times n_{\text{filts},y} \times n_{\text{filts},z} \times n_{h,x} \times n_{h,y} \times n_{h,z}`. Filters must have odd number of samples and are assumed to be centered in the middle of the filter support. ihx : :obj:`tuple` @@ -55,7 +55,7 @@ class NonStationaryConvolve3D(LinearOperator): that the filters must be regularly sampled, i.e. :math:`dh_z=\text{diff}(ihz)=\text{const.}` engine : :obj:`str`, optional Engine used for spread computation (``numpy``, ``numba``, or ``cuda``) - dim_block : :obj:`tuple`, optional + num_threads_per_blocks : :obj:`tuple`, optional Number of threads in each block (only when ``engine=cuda``) dtype : :obj:`str`, optional Type of elements in input array. @@ -78,49 +78,54 @@ class NonStationaryConvolve3D(LinearOperator): If ``ihx``, ``ihy`` or ``ihz`` is not regularly sampled NotImplementedError If ``engine`` is neither ``numpy``, ``fftw``, nor ``scipy``. + + Notes + ----- + See :class:`pylops.signalprocessing.NonStationaryConvolve2D`. + """ def __init__( - self, - dims: Union[int, InputDimsLike], - hs: NDArray, - ihx: InputDimsLike, - ihy: InputDimsLike, - ihz: InputDimsLike, - engine: str = "numpy", - dim_block: Tuple[int, int] = (2, 16, 16), - dtype: DTypeLike = "float64", - name: str = "C") -> None: + self, + dims: Union[int, InputDimsLike], + hs: NDArray, + ihx: InputDimsLike, + ihy: InputDimsLike, + ihz: InputDimsLike, + engine: str = "numpy", + num_threads_per_blocks: Tuple[int, int, int] = (2, 16, 16), + dtype: DTypeLike = "float64", + name: str = "C", + ) -> None: if engine not in ["numpy", "numba", "cuda"]: raise NotImplementedError("engine must be numpy or numba or cuda") if hs.shape[3] % 2 == 0 or hs.shape[4] % 2 == 0 or hs.shape[5] % 2 == 0: raise ValueError("filters hs must have odd length") - #init + if len(np.unique(np.diff(ihx))) > 1 or len(np.unique(np.diff(ihy))) > 1 or len(np.unique(np.diff(ihz))) > 1: + raise ValueError( + "the indices of filters 'ih' are must be regularly sampled" + ) + if min(ihx) < 0 or min(ihy) < 0 or min(ihz) < 0 or max(ihx) >= dims[0] or max(ihy) >= dims[1] or max(ihz) >= dims[2]: + raise ValueError( + "the indices of filters 'ih' must be larger than 0 and smaller than `dims`" + ) self.hs = hs self.hshape = hs.shape[3:] self.ohx, self.dhx, self.nhx = ihx[0], ihx[1] - ihx[0], len(ihx) self.ohy, self.dhy, self.nhy = ihy[0], ihy[1] - ihy[0], len(ihy) self.ohz, self.dhz, self.nhz = ihz[0], ihz[1] - ihz[0], len(ihz) self.ehx, self.ehx, self.ehz = ihx[-1], ihy[-1], ihz[-1] - # self.dims = tuple(dims) self.dims = _value_or_sized_to_tuple(dims) self.engine = engine - # self.dtype = dtype - # self.shape = dims - super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dims, name=name) - # create additional input parameters for engine=cuda self.kwargs_cuda = {} if engine == "cuda": - self.kwargs_cuda["dim_block"] = dim_block - - gridx = (self.dims[0] + dim_block[0] - 1)//dim_block[0] - gridy = (self.dims[1] + dim_block[1] - 1)//dim_block[1] - gridz = (self.dims[2] + dim_block[2] - 1)//dim_block[2] - - self.kwargs_cuda["dim_grid"] = (gridx, gridy, gridz) - + self.kwargs_cuda["num_threads_per_blocks"] = num_threads_per_blocks + num_blocks_x = (self.dims[0] + num_threads_per_blocks[0] - 1) // num_threads_per_blocks[0] + num_blocks_y = (self.dims[1] + num_threads_per_blocks[1] - 1) // num_threads_per_blocks[1] + num_blocks_z = (self.dims[2] + num_threads_per_blocks[2] - 1) // num_threads_per_blocks[2] + self.kwargs_cuda["num_blocks"] = (num_blocks_x, num_blocks_y, num_blocks_z) self._register_multiplications(engine) def _register_multiplications(self, engine: str) -> None: @@ -128,7 +133,7 @@ def _register_multiplications(self, engine: str) -> None: numba_opts = dict(nopython=True, fastmath=True, nogil=True, parallel=True) self._mvrmv = jit(**numba_opts)(self._matvec_rmatvec) elif engine == "cuda": - self._mvrmv = _Convolve3D_kernel + self._mvrmv = _matvec_rmatvec_cuda_call else: self._mvrmv = self._matvec_rmatvec @@ -154,13 +159,13 @@ def _matvec_rmatvec( for iy in range(xdims[1]): for iz in range(xdims[2]): # find closest filters and interpolate h - ihx_l = int(np.floor((ix - ohx) / dhx)) #id number of left for hs_arr - ihy_b = int(np.floor((iy - ohy) / dhy)) #id number of back for hs_arr - ihz_t = int(np.floor((iz - ohz) / dhz)) #id number of top for hs_arr + ihx_l = int(np.floor((ix - ohx) / dhx)) # id number of left for hs_arr + ihy_b = int(np.floor((iy - ohy) / dhy)) # id number of back for hs_arr + ihz_t = int(np.floor((iz - ohz) / dhz)) # id number of top for hs_arr - dhx_r = (ix - ohx) / dhx - ihx_l #weight for right psfs, left 1-ihz_t - dhy_f = (iy - ohy) / dhy - ihy_b #weight for front psfs, left 1-ihz_t - dhz_d = (iz - ohz) / dhz - ihz_t #weight for down psfs, top 1-dhz_d + dhx_r = (ix - ohx) / dhx - ihx_l # weight for right psfs, left 1-ihz_t + dhy_f = (iy - ohy) / dhy - ihy_b # weight for front psfs, left 1-ihz_t + dhz_d = (iz - ohz) / dhz - ihz_t # weight for down psfs, top 1-dhz_d if ihx_l < 0: ihx_l = ihx_r = 0 @@ -226,8 +231,7 @@ def _matvec_rmatvec( max(0, iz - hshape[2] // 2), min(iz + hshape[2] // 2 + 1, xdims[2]), ) - - + # find extremes of h (in case h is going out of model) hxextremes = ( max(0, -ix + hshape[0] // 2), @@ -254,12 +258,9 @@ def _matvec_rmatvec( ) return y - # def _forward(self, x): - @reshaped #reshape to 1D array + @reshaped def _matvec(self, x: NDArray) -> NDArray: ncp = get_array_module(x) - # ncp = cp.get_array_module(x) - # y = ncp.zeros(self.dims, dtype=x.dtype) y = ncp.zeros(self.dims, dtype=self.dtype) y = self._mvrmv( x, @@ -281,13 +282,9 @@ def _matvec(self, x: NDArray) -> NDArray: ) return y - - # def _adjoint(self, x): - @reshaped #reshape to 1D array + @reshaped def _rmatvec(self, x: NDArray) -> NDArray: ncp = get_array_module(x) - # ncp = cp.get_array_module(x) - # y = ncp.zeros(self.dims, dtype=x.dtype) y = ncp.zeros(self.dims, dtype=self.dtype) y = self._mvrmv( x, diff --git a/pytests/test_nonstatconvolve.py b/pytests/test_nonstatconvolve.py index 0a90c489..201567fd 100755 --- a/pytests/test_nonstatconvolve.py +++ b/pytests/test_nonstatconvolve.py @@ -7,8 +7,10 @@ from pylops.signalprocessing import ( Convolve1D, Convolve2D, + ConvolveND, NonStationaryConvolve1D, NonStationaryConvolve2D, + NonStationaryConvolve3D, NonStationaryFilters1D, NonStationaryFilters2D, ) @@ -16,10 +18,15 @@ # filters nfilts = (5, 7) +nfilts3 = (5, 5, 7) + h1 = triang(nfilts[0], sym=True) h2 = np.outer(triang(nfilts[0], sym=True), triang(nfilts[1], sym=True)) +h3 = np.outer(triang(nfilts[0], sym=True), np.outer(triang(nfilts[0], sym=True), triang(nfilts[1], sym=True)[np.newaxis]).T).reshape(nfilts3) + h1stat = np.vstack([h1, h1, h1]) h1ns = np.vstack([h1, -h1, 2 * h1]) + h2stat = np.vstack( [ h2.ravel(), @@ -41,6 +48,25 @@ ] ).reshape(3, 2, nfilts[0], nfilts[1]) +h3stat = np.vstack( + [ + h3.ravel(), + ] * 8 +).reshape(2, 2, 2, nfilts[0], nfilts[0], nfilts[1]) +h3ns = np.vstack( + [ + 2 * h3.ravel(), + h3.ravel(), + h3.ravel(), + h3.ravel(), + h3.ravel(), + h3.ravel(), + -h3.ravel(), + 2 * h3.ravel(), + ] +).reshape(2, 2, 2, nfilts[0], nfilts[0], nfilts[1]) + + par1_1d = { "nz": 21, "nx": 31, @@ -275,3 +301,51 @@ def test_NonStationaryFilters2D(par): h2lsqr = lsqr(Cop, Cop * h2ns.ravel(), damp=1e-20, iter_lim=400, show=0)[0] assert_array_almost_equal(h2ns.ravel(), h2lsqr, decimal=1) + + +@pytest.mark.parametrize("par", [(par_2d)]) +def test_NonStationaryConvolve3D(par): + """Dot-test and inversion for NonStationaryConvolve3D operator""" + Cop = NonStationaryConvolve3D( + dims=(par["nx"], par["nx"], par["nz"]), + hs=h3ns, + ihx=(int(par["nx"] // 4), int(3 * par["nx"] // 4)), + ihy=(int(par["nx"] // 4), int(3 * par["nx"] // 4)), + ihz=(int(par["nz"] // 4), int(3 * par["nz"] // 4)), + dtype="float64", + ) + assert dottest(Cop, par["nx"] * par["nx"] * par["nz"], par["nx"] * par["nx"] * par["nz"]) + + x = np.zeros((par["nx"], par["nx"], par["nz"])) + x[ + int(par["nx"] / 2 - 3) : int(par["nx"] / 2 + 3), + int(par["nx"] / 2 - 3) : int(par["nx"] / 2 + 3), + int(par["nz"] / 2 - 3) : int(par["nz"] / 2 + 3), + ] = 1.0 + x = x.ravel() + xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=40, atol=1e-8, btol=1e-8, show=0)[0] + # given the size of the problem, we can only run few iterations and test accuracy up to 30% + assert np.linalg.norm(x - xlsqr) / np.linalg.norm(x) < 0.3 + + +@pytest.mark.parametrize("par", [(par_2d)]) +def test_StationaryConvolve3D(par): + """Check that Convolve3D and NonStationaryConvolve3D return same result for + stationary filter""" + Cop = NonStationaryConvolve3D( + dims=(par["nx"], par["nx"], par["nz"]), + hs=h3stat, + ihx=(int(par["nx"] // 4), int(3 * par["nx"] // 4)), + ihy=(int(par["nx"] // 4), int(3 * par["nx"] // 4)), + ihz=(int(par["nz"] // 4), int(3 * par["nz"] // 4)), + dtype="float64", + ) + Cop_stat = ConvolveND( + dims=(par["nx"], par["nx"], par["nz"]), + h=h3, + offset=(nfilts[0] // 2, nfilts[0] // 2, nfilts[1] // 2), + dtype="float64", + ) + x = np.random.normal(0, 1, (par["nx"], par["nx"], par["nz"])) + + assert_array_almost_equal(Cop_stat * x, Cop * x, decimal=10) From 65dcd9c656a97963f9e174890d8ca929883ea5c7 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Thu, 20 Jul 2023 09:13:17 +0300 Subject: [PATCH 37/69] minor: fixed flake8 issues --- .../_nonstatconvolve3d_cuda.py | 43 +++---- pylops/signalprocessing/nonstatconvolve3d.py | 117 ++++++++++++------ 2 files changed, 104 insertions(+), 56 deletions(-) diff --git a/pylops/signalprocessing/_nonstatconvolve3d_cuda.py b/pylops/signalprocessing/_nonstatconvolve3d_cuda.py index a5810f03..7d84a988 100644 --- a/pylops/signalprocessing/_nonstatconvolve3d_cuda.py +++ b/pylops/signalprocessing/_nonstatconvolve3d_cuda.py @@ -4,7 +4,9 @@ @cuda.jit(max_registers=40) -def _matvec_rmatvec(x, y, hs, hshape, xdims, ohx, ohy, ohz, dhx, dhy, dhz, nhx, nhy, nhz, rmatvec): +def _matvec_rmatvec( + x, y, hs, hshape, xdims, ohx, ohy, ohz, dhx, dhy, dhz, nhx, nhy, nhz, rmatvec +): """Cuda kernels for NonStationaryConvolve3D operator Cuda implementation of matvec and rmatvec for NonStationaryConvolve3D operator. See @@ -14,16 +16,16 @@ def _matvec_rmatvec(x, y, hs, hshape, xdims, ohx, ohy, ohz, dhx, dhy, dhz, nhx, ix, iy, iz = cuda.grid(3) if ix < xdims[0] and iy < xdims[1] and iz < xdims[2]: - + ## find closest filters and interpolate h - ihx_l = int(floor((ix - ohx) / dhx)) #id number of left for hs_arr - ihy_b = int(floor((iy - ohy) / dhy)) #id number of back for hs_arr - ihz_t = int(floor((iz - ohz) / dhz)) #id number of top for hs_arr - - dhx_r = (ix - ohx) / dhx - ihx_l #weight for right psfs, left 1-dhx_r - dhy_f = (iy - ohy) / dhy - ihy_b #weight for front psfs, back 1-dhy_f - dhz_d = (iz - ohz) / dhz - ihz_t #weight for down psfs, top 1-dhz_d - + ihx_l = int(floor((ix - ohx) / dhx)) # id number of left for hs_arr + ihy_b = int(floor((iy - ohy) / dhy)) # id number of back for hs_arr + ihz_t = int(floor((iz - ohz) / dhz)) # id number of top for hs_arr + + dhx_r = (ix - ohx) / dhx - ihx_l # weight for right psfs, left 1-dhx_r + dhy_f = (iy - ohy) / dhy - ihy_b # weight for front psfs, back 1-dhy_f + dhz_d = (iz - ohz) / dhz - ihz_t # weight for down psfs, top 1-dhz_d + if ihx_l < 0: ihx_l = ihx_r = 0 dhx_l = dhx_r = 0.5 @@ -33,7 +35,7 @@ def _matvec_rmatvec(x, y, hs, hshape, xdims, ohx, ohy, ohz, dhx, dhy, dhz, nhx, else: ihx_r = ihx_l + 1 dhx_l = 1.0 - dhx_r - + if ihy_b < 0: ihy_b = ihy_f = 0 dhy_b = dhy_f = 0.5 @@ -58,13 +60,13 @@ def _matvec_rmatvec(x, y, hs, hshape, xdims, ohx, ohy, ohz, dhx, dhy, dhz, nhx, h_lbd = hs[ihx_l, ihy_b, ihz_d] h_lft = hs[ihx_l, ihy_f, ihz_t] h_lfd = hs[ihx_l, ihy_f, ihz_d] - + h_rbt = hs[ihx_r, ihy_b, ihz_t] h_rbd = hs[ihx_r, ihy_b, ihz_d] h_rft = hs[ihx_r, ihy_f, ihz_t] h_rfd = hs[ihx_r, ihy_f, ihz_d] - ## find extremes of model where to apply h (in case h is going out of model) + ## find extremes of model where to apply h (in case h is going out of model) xextremes = ( max(0, ix - hshape[0] // 2), min(ix + hshape[0] // 2 + 1, xdims[0]), @@ -77,8 +79,7 @@ def _matvec_rmatvec(x, y, hs, hshape, xdims, ohx, ohy, ohz, dhx, dhy, dhz, nhx, max(0, iz - hshape[2] // 2), min(iz + hshape[2] // 2 + 1, xdims[2]), ) - - + ## find extremes of h (in case h is going out of model) hxextremes = ( max(0, -ix + hshape[0] // 2), @@ -92,7 +93,7 @@ def _matvec_rmatvec(x, y, hs, hshape, xdims, ohx, ohy, ohz, dhx, dhy, dhz, nhx, max(0, -iz + hshape[2] // 2), min(hshape[2], hshape[2] // 2 + (xdims[2] - iz)), ) - + # place filter in output for ixx, hxx in zip( range(xextremes[0], xextremes[1]), range(hxextremes[0], hxextremes[1]) @@ -105,7 +106,7 @@ def _matvec_rmatvec(x, y, hs, hshape, xdims, ohx, ohy, ohz, dhx, dhy, dhz, nhx, range(hzextremes[0], hzextremes[1]), ): h = ( - dhx_l * dhy_b * dhz_t * h_lbt[hxx, hyy, hzz] + dhx_l * dhy_b * dhz_t * h_lbt[hxx, hyy, hzz] + dhx_l * dhy_b * dhz_d * h_lbd[hxx, hyy, hzz] + dhx_l * dhy_f * dhz_t * h_lft[hxx, hyy, hzz] + dhx_l * dhy_f * dhz_d * h_lfd[hxx, hyy, hzz] @@ -114,7 +115,7 @@ def _matvec_rmatvec(x, y, hs, hshape, xdims, ohx, ohy, ohz, dhx, dhy, dhz, nhx, + dhx_r * dhy_f * dhz_t * h_rft[hxx, hyy, hzz] + dhx_r * dhy_f * dhz_d * h_rfd[hxx, hyy, hzz] ) - + if rmatvec: y[ix, iy, iz] += h * x[ixx, iyy, izz] else: @@ -136,9 +137,9 @@ def _matvec_rmatvec_call( nhx, nhy, nhz, - rmatvec= False, - num_blocks=(1,32,16), - num_threads_per_blocks=(24,24,24), + rmatvec=False, + num_blocks=(1, 32, 16), + num_threads_per_blocks=(24, 24, 24), ): """Caller for NonStationaryConvolve3D operator diff --git a/pylops/signalprocessing/nonstatconvolve3d.py b/pylops/signalprocessing/nonstatconvolve3d.py index 1b7eaf3b..853bb144 100644 --- a/pylops/signalprocessing/nonstatconvolve3d.py +++ b/pylops/signalprocessing/nonstatconvolve3d.py @@ -41,7 +41,8 @@ class NonStationaryConvolve3D(LinearOperator): Number of samples for each dimension (which we refer to as :math`n_x \times n_y \times n_z`). hs : :obj:`numpy.ndarray` Bank of 3d compact filters of size - :math:`n_{\text{filts},x} \times n_{\text{filts},y} \times n_{\text{filts},z} \times n_{h,x} \times n_{h,y} \times n_{h,z}`. + :math:`n_{\text{filts},x} \times n_{\text{filts},y} \times + n_{\text{filts},z} \times n_{h,x} \times n_{h,y} \times n_{h,z}`. Filters must have odd number of samples and are assumed to be centered in the middle of the filter support. ihx : :obj:`tuple` @@ -83,7 +84,8 @@ class NonStationaryConvolve3D(LinearOperator): ----- See :class:`pylops.signalprocessing.NonStationaryConvolve2D`. - """ + """ + def __init__( self, dims: Union[int, InputDimsLike], @@ -100,34 +102,51 @@ def __init__( raise NotImplementedError("engine must be numpy or numba or cuda") if hs.shape[3] % 2 == 0 or hs.shape[4] % 2 == 0 or hs.shape[5] % 2 == 0: raise ValueError("filters hs must have odd length") - if len(np.unique(np.diff(ihx))) > 1 or len(np.unique(np.diff(ihy))) > 1 or len(np.unique(np.diff(ihz))) > 1: + if ( + len(np.unique(np.diff(ihx))) > 1 + or len(np.unique(np.diff(ihy))) > 1 + or len(np.unique(np.diff(ihz))) > 1 + ): raise ValueError( "the indices of filters 'ih' are must be regularly sampled" ) - if min(ihx) < 0 or min(ihy) < 0 or min(ihz) < 0 or max(ihx) >= dims[0] or max(ihy) >= dims[1] or max(ihz) >= dims[2]: + if ( + min(ihx) < 0 + or min(ihy) < 0 + or min(ihz) < 0 + or max(ihx) >= dims[0] + or max(ihy) >= dims[1] + or max(ihz) >= dims[2] + ): raise ValueError( "the indices of filters 'ih' must be larger than 0 and smaller than `dims`" ) self.hs = hs self.hshape = hs.shape[3:] - self.ohx, self.dhx, self.nhx = ihx[0], ihx[1] - ihx[0], len(ihx) - self.ohy, self.dhy, self.nhy = ihy[0], ihy[1] - ihy[0], len(ihy) - self.ohz, self.dhz, self.nhz = ihz[0], ihz[1] - ihz[0], len(ihz) + self.ohx, self.dhx, self.nhx = ihx[0], ihx[1] - ihx[0], len(ihx) + self.ohy, self.dhy, self.nhy = ihy[0], ihy[1] - ihy[0], len(ihy) + self.ohz, self.dhz, self.nhz = ihz[0], ihz[1] - ihz[0], len(ihz) self.ehx, self.ehx, self.ehz = ihx[-1], ihy[-1], ihz[-1] self.dims = _value_or_sized_to_tuple(dims) self.engine = engine super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dims, name=name) - + # create additional input parameters for engine=cuda self.kwargs_cuda = {} if engine == "cuda": self.kwargs_cuda["num_threads_per_blocks"] = num_threads_per_blocks - num_blocks_x = (self.dims[0] + num_threads_per_blocks[0] - 1) // num_threads_per_blocks[0] - num_blocks_y = (self.dims[1] + num_threads_per_blocks[1] - 1) // num_threads_per_blocks[1] - num_blocks_z = (self.dims[2] + num_threads_per_blocks[2] - 1) // num_threads_per_blocks[2] + num_blocks_x = ( + self.dims[0] + num_threads_per_blocks[0] - 1 + ) // num_threads_per_blocks[0] + num_blocks_y = ( + self.dims[1] + num_threads_per_blocks[1] - 1 + ) // num_threads_per_blocks[1] + num_blocks_z = ( + self.dims[2] + num_threads_per_blocks[2] - 1 + ) // num_threads_per_blocks[2] self.kwargs_cuda["num_blocks"] = (num_blocks_x, num_blocks_y, num_blocks_z) self._register_multiplications(engine) - + def _register_multiplications(self, engine: str) -> None: if engine == "numba": numba_opts = dict(nopython=True, fastmath=True, nogil=True, parallel=True) @@ -136,7 +155,7 @@ def _register_multiplications(self, engine: str) -> None: self._mvrmv = _matvec_rmatvec_cuda_call else: self._mvrmv = self._matvec_rmatvec - + @staticmethod def _matvec_rmatvec( x: NDArray, @@ -159,14 +178,26 @@ def _matvec_rmatvec( for iy in range(xdims[1]): for iz in range(xdims[2]): # find closest filters and interpolate h - ihx_l = int(np.floor((ix - ohx) / dhx)) # id number of left for hs_arr - ihy_b = int(np.floor((iy - ohy) / dhy)) # id number of back for hs_arr - ihz_t = int(np.floor((iz - ohz) / dhz)) # id number of top for hs_arr - - dhx_r = (ix - ohx) / dhx - ihx_l # weight for right psfs, left 1-ihz_t - dhy_f = (iy - ohy) / dhy - ihy_b # weight for front psfs, left 1-ihz_t - dhz_d = (iz - ohz) / dhz - ihz_t # weight for down psfs, top 1-dhz_d - + ihx_l = int( + np.floor((ix - ohx) / dhx) + ) # id number of left for hs_arr + ihy_b = int( + np.floor((iy - ohy) / dhy) + ) # id number of back for hs_arr + ihz_t = int( + np.floor((iz - ohz) / dhz) + ) # id number of top for hs_arr + + dhx_r = ( + ix - ohx + ) / dhx - ihx_l # weight for right psfs, left 1-ihz_t + dhy_f = ( + iy - ohy + ) / dhy - ihy_b # weight for front psfs, left 1-ihz_t + dhz_d = ( + iz - ohz + ) / dhz - ihz_t # weight for down psfs, top 1-dhz_d + if ihx_l < 0: ihx_l = ihx_r = 0 dhx_l = dhx_r = 0.5 @@ -176,7 +207,7 @@ def _matvec_rmatvec( else: ihx_r = ihx_l + 1 dhx_l = 1.0 - dhx_r - + if ihy_b < 0: ihy_b = ihy_f = 0 dhy_b = dhy_f = 0.5 @@ -186,7 +217,7 @@ def _matvec_rmatvec( else: ihy_f = ihy_b + 1 dhy_b = 1.0 - dhy_f - + if ihz_t < 0: ihz_t = ihz_d = 0 dhz_t = dhz_d = 0.5 @@ -196,19 +227,19 @@ def _matvec_rmatvec( else: ihz_d = ihz_t + 1 dhz_t = 1.0 - dhz_d - + h_lbt = hs[ihx_l, ihy_b, ihz_t] h_lbd = hs[ihx_l, ihy_b, ihz_d] h_lft = hs[ihx_l, ihy_f, ihz_t] h_lfd = hs[ihx_l, ihy_f, ihz_d] - + h_rbt = hs[ihx_r, ihy_b, ihz_t] h_rbd = hs[ihx_r, ihy_b, ihz_d] h_rft = hs[ihx_r, ihy_f, ihz_t] h_rfd = hs[ihx_r, ihy_f, ihz_d] - + h = ( - dhx_l * dhy_b * dhz_t * h_lbt + dhx_l * dhy_b * dhz_t * h_lbt + dhx_l * dhy_b * dhz_d * h_lbd + dhx_l * dhy_f * dhz_t * h_lft + dhx_l * dhy_f * dhz_d * h_lfd @@ -217,7 +248,7 @@ def _matvec_rmatvec( + dhx_r * dhy_f * dhz_t * h_rft + dhx_r * dhy_f * dhz_d * h_rfd ) - + # find extremes of model where to apply h (in case h is going out of model) xextremes = ( max(0, ix - hshape[0] // 2), @@ -245,19 +276,35 @@ def _matvec_rmatvec( max(0, -iz + hshape[2] // 2), min(hshape[2], hshape[2] // 2 + (xdims[2] - iz)), ) - + if not rmatvec: - y[ xextremes[0] : xextremes[1], yextremes[0] : yextremes[1], zextremes[0] : zextremes[1] ] += ( + y[ + xextremes[0] : xextremes[1], + yextremes[0] : yextremes[1], + zextremes[0] : zextremes[1], + ] += ( x[ix, iy, iz] - * h[ hxextremes[0] : hxextremes[1], hyextremes[0] : hyextremes[1], hzextremes[0] : hzextremes[1] ] + * h[ + hxextremes[0] : hxextremes[1], + hyextremes[0] : hyextremes[1], + hzextremes[0] : hzextremes[1], + ] ) else: y[ix, iy, iz] = np.sum( - h[ hxextremes[0] : hxextremes[1], hyextremes[0] : hyextremes[1], hzextremes[0] : hzextremes[1] ] - * x[ xextremes[0] : xextremes[1], yextremes[0] : yextremes[1], zextremes[0] : zextremes[1] ] + h[ + hxextremes[0] : hxextremes[1], + hyextremes[0] : hyextremes[1], + hzextremes[0] : hzextremes[1], + ] + * x[ + xextremes[0] : xextremes[1], + yextremes[0] : yextremes[1], + zextremes[0] : zextremes[1], + ] ) - return y - + return y + @reshaped def _matvec(self, x: NDArray) -> NDArray: ncp = get_array_module(x) From 53e727cfad71b570ca874eb19aa06d48b5eca7c8 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Thu, 20 Jul 2023 09:16:26 +0300 Subject: [PATCH 38/69] minor: fixed flake8 issues in test --- .../_nonstatconvolve3d_cuda.py | 6 ++--- pytests/test_nonstatconvolve.py | 26 ++++++++++++++----- 2 files changed, 22 insertions(+), 10 deletions(-) diff --git a/pylops/signalprocessing/_nonstatconvolve3d_cuda.py b/pylops/signalprocessing/_nonstatconvolve3d_cuda.py index 7d84a988..5b309050 100644 --- a/pylops/signalprocessing/_nonstatconvolve3d_cuda.py +++ b/pylops/signalprocessing/_nonstatconvolve3d_cuda.py @@ -17,7 +17,7 @@ def _matvec_rmatvec( if ix < xdims[0] and iy < xdims[1] and iz < xdims[2]: - ## find closest filters and interpolate h + # find closest filters and interpolate h ihx_l = int(floor((ix - ohx) / dhx)) # id number of left for hs_arr ihy_b = int(floor((iy - ohy) / dhy)) # id number of back for hs_arr ihz_t = int(floor((iz - ohz) / dhz)) # id number of top for hs_arr @@ -66,7 +66,7 @@ def _matvec_rmatvec( h_rft = hs[ihx_r, ihy_f, ihz_t] h_rfd = hs[ihx_r, ihy_f, ihz_d] - ## find extremes of model where to apply h (in case h is going out of model) + # find extremes of model where to apply h (in case h is going out of model) xextremes = ( max(0, ix - hshape[0] // 2), min(ix + hshape[0] // 2 + 1, xdims[0]), @@ -80,7 +80,7 @@ def _matvec_rmatvec( min(iz + hshape[2] // 2 + 1, xdims[2]), ) - ## find extremes of h (in case h is going out of model) + # find extremes of h (in case h is going out of model) hxextremes = ( max(0, -ix + hshape[0] // 2), min(hshape[0], hshape[0] // 2 + (xdims[0] - ix)), diff --git a/pytests/test_nonstatconvolve.py b/pytests/test_nonstatconvolve.py index 201567fd..688fe39d 100755 --- a/pytests/test_nonstatconvolve.py +++ b/pytests/test_nonstatconvolve.py @@ -22,7 +22,10 @@ h1 = triang(nfilts[0], sym=True) h2 = np.outer(triang(nfilts[0], sym=True), triang(nfilts[1], sym=True)) -h3 = np.outer(triang(nfilts[0], sym=True), np.outer(triang(nfilts[0], sym=True), triang(nfilts[1], sym=True)[np.newaxis]).T).reshape(nfilts3) +h3 = np.outer( + triang(nfilts[0], sym=True), + np.outer(triang(nfilts[0], sym=True), triang(nfilts[1], sym=True)[np.newaxis]).T, +).reshape(nfilts3) h1stat = np.vstack([h1, h1, h1]) h1ns = np.vstack([h1, -h1, 2 * h1]) @@ -51,7 +54,8 @@ h3stat = np.vstack( [ h3.ravel(), - ] * 8 + ] + * 8 ).reshape(2, 2, 2, nfilts[0], nfilts[0], nfilts[1]) h3ns = np.vstack( [ @@ -170,7 +174,9 @@ def test_NonStationaryConvolve1D(par): x = np.zeros((par["nx"])) x[par["nx"] // 2] = 1.0 - xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr( + Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0 + )[0] assert_array_almost_equal(x, xlsqr, decimal=1) # 1D on 2D @@ -190,7 +196,9 @@ def test_NonStationaryConvolve1D(par): int(par["nz"] / 2 - 3) : int(par["nz"] / 2 + 3), ] = 1.0 x = x.ravel() - xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=400, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=400, atol=1e-8, btol=1e-8, show=0)[ + 0 + ] assert_array_almost_equal(x, xlsqr, decimal=1) @@ -234,7 +242,9 @@ def test_NonStationaryConvolve2D(par): int(par["nz"] / 2 - 3) : int(par["nz"] / 2 + 3), ] = 1.0 x = x.ravel() - xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=400, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=400, atol=1e-8, btol=1e-8, show=0)[ + 0 + ] assert_array_almost_equal(x, xlsqr, decimal=1) @@ -266,7 +276,7 @@ def test_StationaryConvolve2D(par): (par1_1d), ], ) -def test_NonStationaryFilters2D(par): +def test_NonStationaryFilters1D(par): """Dot-test and inversion for NonStationaryFilters2D operator""" x = np.zeros((par["nx"])) x[par["nx"] // 4], x[par["nx"] // 2], x[3 * par["nx"] // 4] = 1.0, 1.0, 1.0 @@ -314,7 +324,9 @@ def test_NonStationaryConvolve3D(par): ihz=(int(par["nz"] // 4), int(3 * par["nz"] // 4)), dtype="float64", ) - assert dottest(Cop, par["nx"] * par["nx"] * par["nz"], par["nx"] * par["nx"] * par["nz"]) + assert dottest( + Cop, par["nx"] * par["nx"] * par["nz"], par["nx"] * par["nx"] * par["nz"] + ) x = np.zeros((par["nx"], par["nx"], par["nz"])) x[ From adfd2083d37fdb992970b0382cc859d894185b5b Mon Sep 17 00:00:00 2001 From: mrava87 Date: Thu, 20 Jul 2023 10:19:35 +0300 Subject: [PATCH 39/69] bug: fix wrong import in nonstatconvolve3d.py --- pylops/signalprocessing/nonstatconvolve3d.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pylops/signalprocessing/nonstatconvolve3d.py b/pylops/signalprocessing/nonstatconvolve3d.py index 853bb144..026423a1 100644 --- a/pylops/signalprocessing/nonstatconvolve3d.py +++ b/pylops/signalprocessing/nonstatconvolve3d.py @@ -17,7 +17,7 @@ if jit_message is None: from numba import jit, prange - from ._nonstatconvolve2d_cuda import ( + from ._nonstatconvolve3d_cuda import ( _matvec_rmatvec_call as _matvec_rmatvec_cuda_call, ) @@ -110,6 +110,7 @@ def __init__( raise ValueError( "the indices of filters 'ih' are must be regularly sampled" ) + if ( min(ihx) < 0 or min(ihy) < 0 From d88aaf18bef0a79fce85eb76482d59e5752bc045 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Thu, 20 Jul 2023 11:18:46 +0300 Subject: [PATCH 40/69] doc: fixed some typos in nonstatconvolve files --- pylops/signalprocessing/nonstatconvolve2d.py | 6 +++--- pylops/signalprocessing/nonstatconvolve3d.py | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/pylops/signalprocessing/nonstatconvolve2d.py b/pylops/signalprocessing/nonstatconvolve2d.py index 8cd025c1..a7056986 100644 --- a/pylops/signalprocessing/nonstatconvolve2d.py +++ b/pylops/signalprocessing/nonstatconvolve2d.py @@ -36,12 +36,12 @@ class NonStationaryConvolve2D(LinearOperator): Apply non-stationary two-dimensional convolution. A varying compact filter is provided on a coarser grid and on-the-fly interpolation is applied in forward and adjoint modes. - Both input and output have size :math`n_x \times n_z`. + Both input and output have size :math:`n_x \times n_z`. Parameters ---------- dims : :obj:`list` or :obj:`int` - Number of samples for each dimension (which we refer to as :math`n_x \times n_z`). + Number of samples for each dimension (which we refer to as :math:`n_x \times n_z`). hs : :obj:`numpy.ndarray` Bank of 2d compact filters of size :math:`n_{\text{filts},x} \times n_{\text{filts},z} \times n_{h,x} \times n_{h,z}`. @@ -424,7 +424,7 @@ def _register_multiplications(self, engine: str) -> None: self._mv = self.__matvec self._rmv = self.__rmatvec - # use same matvec method as inNonStationaryConvolve2D + # use same matvec method as in NonStationaryConvolve2D __matvec = staticmethod(NonStationaryConvolve2D._matvec_rmatvec) @staticmethod diff --git a/pylops/signalprocessing/nonstatconvolve3d.py b/pylops/signalprocessing/nonstatconvolve3d.py index 026423a1..236851f8 100644 --- a/pylops/signalprocessing/nonstatconvolve3d.py +++ b/pylops/signalprocessing/nonstatconvolve3d.py @@ -33,12 +33,12 @@ class NonStationaryConvolve3D(LinearOperator): Apply non-stationary three-dimensional convolution. A varying compact filter is provided on a coarser grid and on-the-fly interpolation is applied - in forward and adjoint modes. Both input and output have size :math`n_x \times n_y \times n_z`. + in forward and adjoint modes. Both input and output have size :math:`n_x \times n_y \times n_z`. Parameters ---------- dims : :obj:`list` or :obj:`int` - Number of samples for each dimension (which we refer to as :math`n_x \times n_y \times n_z`). + Number of samples for each dimension (which we refer to as :math:`n_x \times n_y \times n_z`). hs : :obj:`numpy.ndarray` Bank of 3d compact filters of size :math:`n_{\text{filts},x} \times n_{\text{filts},y} \times From 8ae98bb09c3dd17339f578d0aa2bd38e7b5db400 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Thu, 20 Jul 2023 14:17:36 +0300 Subject: [PATCH 41/69] bug: add forceflat to restriction to avoid unexpected unflattening of adjoint --- pylops/basicoperators/restriction.py | 26 +++++++++++++++++++++++++- 1 file changed, 25 insertions(+), 1 deletion(-) diff --git a/pylops/basicoperators/restriction.py b/pylops/basicoperators/restriction.py index 4aaf10ea..c2e51a31 100644 --- a/pylops/basicoperators/restriction.py +++ b/pylops/basicoperators/restriction.py @@ -1,5 +1,7 @@ __all__ = ["Restriction"] +import logging + from typing import Sequence, Union import numpy as np @@ -11,6 +13,8 @@ from pylops.utils.backend import get_array_module, to_cupy_conditional from pylops.utils.typing import DTypeLike, InputDimsLike, IntNDArray, NDArray +logging.basicConfig(format="%(levelname)s: %(message)s", level=logging.WARNING) + def _compute_iavamask(dims, axis, iava, ncp): """Compute restriction mask when using cupy arrays""" @@ -46,6 +50,12 @@ class Restriction(LinearOperator): Work inplace (``True``) or make a new copy (``False``). By default, data is a reference to the model (in forward) and model is a reference to the data (in adjoint). + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after rmatvec. Note that this is only + required when `len(dims)=2`, otherwise pylops will detect whether to + return a 1d or nd array. dtype : :obj:`str`, optional Type of elements in input array. name : :obj:`str`, optional @@ -97,6 +107,7 @@ def __init__( iava: Union[IntNDArray, Sequence[int]], axis: int = -1, inplace: bool = True, + forceflat: bool = None, dtype: DTypeLike = "float64", name: str = "R", ) -> None: @@ -105,7 +116,20 @@ def __init__( axis = normalize_axis_index(axis, len(dims)) dimsd = list(dims) # data dimensions dimsd[axis] = len(iava) - super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dimsd, name=name) + + # check if forceflat is needed and set it back to None otherwise + if len(dims) > 2: + if forceflat is not None: + logging.warning( + f"setting forceflat=None since len(dims)={len(dims)}>2. " + f"PyLops will automatically detect whether to return " + f"a 1d or nd array based on the shape of the input" + f"array." + ) + forceflat = None + + super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dimsd, + forceflat=forceflat, name=name) iavareshape = np.ones(len(self.dims), dtype=int) iavareshape[axis] = len(iava) From 5e263f7f4c87705c3248ae84d3e4f27198df5c48 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Mon, 7 Aug 2023 15:15:25 +0200 Subject: [PATCH 42/69] doc: updated readthedocs yaml file --- .readthedocs.yaml | 23 +++++++++++++++++++++++ readthedocs.yml | 14 -------------- 2 files changed, 23 insertions(+), 14 deletions(-) create mode 100644 .readthedocs.yaml delete mode 100644 readthedocs.yml diff --git a/.readthedocs.yaml b/.readthedocs.yaml new file mode 100644 index 00000000..146206bb --- /dev/null +++ b/.readthedocs.yaml @@ -0,0 +1,23 @@ +# .readthedocs.yaml +# Read the Docs configuration file +# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details + +# Required +version: 2 + +# Set the version of Python and other tools you might need +build: + os: ubuntu-20.04 + tools: + python: "3.9" + +# Build documentation in the docs/ directory with Sphinx +sphinx: + configuration: docs/source/conf.py + +# Declare the Python requirements required to build your docs +python: + install: + - requirements: requirements-dev.txt + - method: setuptools + path: . diff --git a/readthedocs.yml b/readthedocs.yml deleted file mode 100644 index a4cf3385..00000000 --- a/readthedocs.yml +++ /dev/null @@ -1,14 +0,0 @@ -# .readthedocs.yml - -version: 2 - -build: - os: ubuntu-20.04 - tools: - python: "3.9" - -python: - install: - - requirements: requirements-dev.txt - - method: setuptools - path: . From a64fba1281089e2cc1d5aec6e5939d75b6fafae6 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Mon, 7 Aug 2023 15:36:25 +0200 Subject: [PATCH 43/69] minor: fix flake8 issue in convolve1d --- pylops/signalprocessing/convolve1d.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pylops/signalprocessing/convolve1d.py b/pylops/signalprocessing/convolve1d.py index 9093f526..74057aa8 100644 --- a/pylops/signalprocessing/convolve1d.py +++ b/pylops/signalprocessing/convolve1d.py @@ -97,7 +97,7 @@ def __init__( @reshaped def _matvec(self, x: NDArray) -> NDArray: - if type(self.h) != type(x): + if type(self.h) is not type(x): self.h = to_cupy_conditional(x, self.h) self.convfunc, self.method = _choose_convfunc( self.h, self.method, self.dims @@ -106,7 +106,7 @@ def _matvec(self, x: NDArray) -> NDArray: @reshaped def _rmatvec(self, x: NDArray) -> NDArray: - if type(self.hstar) != type(x): + if type(self.hstar) is not type(x): self.hstar = to_cupy_conditional(x, self.hstar) self.convfunc, self.method = _choose_convfunc( self.hstar, self.method, self.dims @@ -162,7 +162,7 @@ def __init__( @reshaped def _matvec(self, x: NDArray) -> NDArray: - if type(self.h) != type(x): + if type(self.h) is not type(x): self.h = to_cupy_conditional(x, self.h) self.convfunc, self.method = _choose_convfunc( self.h, self.method, self.dims @@ -174,7 +174,7 @@ def _matvec(self, x: NDArray) -> NDArray: @reshaped def _rmatvec(self, x: NDArray) -> NDArray: ncp = get_array_module(x) - if type(self.h) != type(x): + if type(self.h) is not type(x): self.hstar = to_cupy_conditional(x, self.hstar) self.convfunc, self.method = _choose_convfunc( self.hstar, self.method, self.dims From a892de6460a6d0700422558440a5c322cb5c08c7 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Mon, 7 Aug 2023 22:13:55 +0200 Subject: [PATCH 44/69] bug: fixed convmtx as discussed in https://github.com/PyLops/pylops/issues/517 --- pylops/avo/poststack.py | 2 +- pylops/avo/prestack.py | 4 +- pylops/utils/signalprocessing.py | 42 ++++++++++++------ pytests/test_signalutils.py | 74 +++++++++++++++++++++++++------- 4 files changed, 91 insertions(+), 31 deletions(-) diff --git a/pylops/avo/poststack.py b/pylops/avo/poststack.py index 11c49682..7e514707 100644 --- a/pylops/avo/poststack.py +++ b/pylops/avo/poststack.py @@ -102,7 +102,7 @@ def _PoststackLinearModelling( # Create wavelet operator if len(wav.shape) == 1: - C = convmtx(wav, nt0)[:, len(wav) // 2 : -len(wav) // 2 + 1] + C = convmtx(wav, nt0, len(wav) // 2)[:nt0] else: C = nonstationary_convmtx(wav, nt0, hc=wav.shape[1] // 2, pad=(nt0, nt0)) # Combine operators diff --git a/pylops/avo/prestack.py b/pylops/avo/prestack.py index 491dfb9c..8630bc9a 100644 --- a/pylops/avo/prestack.py +++ b/pylops/avo/prestack.py @@ -191,7 +191,7 @@ def PrestackLinearModelling( D = get_block_diag(theta)(*([D] * nG)) # Create wavelet operator - C = ncp.asarray(convmtx(wav, nt0))[:, len(wav) // 2 : -len(wav) // 2 + 1] + C = ncp.asarray(convmtx(wav, nt0, len(wav) // 2)[:nt0]) C = [C] * ntheta C = get_block_diag(theta)(*C) @@ -346,7 +346,7 @@ def PrestackWaveletModelling( M = ncp.dot(G, ncp.dot(D, m.T.ravel())).reshape(ntheta, nt0) Mconv = VStack( [ - MatrixMult(convmtx(M[itheta], nwav)[wavc : -nwav + wavc + 1], dtype=dtype) + MatrixMult(convmtx(M[itheta], nwav, wavc)[:nt0], dtype=dtype) for itheta in range(ntheta) ] ) diff --git a/pylops/utils/signalprocessing.py b/pylops/utils/signalprocessing.py index 17917021..0aa82e3a 100644 --- a/pylops/utils/signalprocessing.py +++ b/pylops/utils/signalprocessing.py @@ -5,6 +5,7 @@ "dip_estimate", ] +import warnings from typing import Tuple import numpy as np @@ -15,39 +16,56 @@ from pylops.utils.typing import NDArray -def convmtx(h: npt.ArrayLike, n: int) -> NDArray: +def convmtx(h: npt.ArrayLike, n: int, offset: int = 0) -> NDArray: r"""Convolution matrix - Equivalent of `MATLAB's convmtx function - `_ . Makes a dense convolution matrix :math:`\mathbf{C}` such that the dot product ``np.dot(C, x)`` is the convolution of - the filter :math:`h` and the input signal :math:`x`. + the filter :math:`h` centered on `offset` and the input signal :math:`x`. + + Equivalent of `MATLAB's convmtx function + `_ for + ``offset=0``. Parameters ---------- h : :obj:`np.ndarray` Convolution filter (1D array) n : :obj:`int` - Number of columns (if :math:`\text{len}(h) < n`) or rows - (if :math:`\text{len}(h) \geq n`) of convolution matrix + Number of columns of convolution matrix + offset : :obj:`int` + Index of the center of the filter Returns ------- C : :obj:`np.ndarray` Convolution matrix of size :math:`\text{len}(h)+n-1 \times n` - (if :math:`\text{len}(h) < n`) or :math:`n \times \text{len}(h)+n-1` - (if :math:`\text{len}(h) \geq n`) """ + warnings.warn( + "A new implementation of convmtx is provided in v2.2.0 to match " + "MATLAB's convmtx method as stated in the docstring. Prior to v2.2.0," + "The implementation of convmtx provided prior to v2.2.0 was instead " + "not consistent with the documentation. Users are highly encouraged " + "to modify their codes accordingly.", + FutureWarning, + ) + ncp = get_array_module(h) - if len(h) < n: - col_1 = ncp.r_[h[0], ncp.zeros(n - 1, dtype=h.dtype)] + nh = len(h) + if nh < n: + h = np.flipud(h) + col_1 = ncp.r_[h[0], ncp.zeros(n + nh - 2, dtype=h.dtype)] row_1 = ncp.r_[h, ncp.zeros(n - 1, dtype=h.dtype)] + C = get_toeplitz(h)(col_1, row_1) + # apply offset + C = C[:, offset : offset + n] else: + col_1 = ncp.r_[h, ncp.zeros(n + nh - 2, dtype=h.dtype)] row_1 = ncp.r_[h[0], ncp.zeros(n - 1, dtype=h.dtype)] - col_1 = ncp.r_[h, ncp.zeros(n - 1, dtype=h.dtype)] - C = get_toeplitz(h)(col_1, row_1) + C = get_toeplitz(h)(col_1, row_1) + # apply offset + C = C[offset : offset + nh + n - 1] return C diff --git a/pytests/test_signalutils.py b/pytests/test_signalutils.py index 79c6d4e8..d7551a2e 100755 --- a/pytests/test_signalutils.py +++ b/pytests/test_signalutils.py @@ -4,47 +4,89 @@ from pylops.utils.signalprocessing import convmtx, nonstationary_convmtx -par1 = {"nt": 51, "imag": 0, "dtype": "float32"} # real -par2 = {"nt": 51, "imag": 1j, "dtype": "complex64"} # complex +par1 = {"nt": 51, "nh": 7, "imag": 0, "dtype": "float32"} # odd sign, odd filt, real +par1j = { + "nt": 51, + "nh": 7, + "imag": 1j, + "dtype": "complex64", +} # odd sign, odd filt, complex +par2 = {"nt": 50, "nh": 7, "imag": 0, "dtype": "float32"} # even sign, odd filt, real +par2j = { + "nt": 50, + "nh": 7, + "imag": 1j, + "dtype": "complex64", +} # even sign, odd filt, complex +par3 = {"nt": 51, "nh": 6, "imag": 0, "dtype": "float32"} # odd sign, even filt, real +par3j = { + "nt": 51, + "nh": 6, + "imag": 1j, + "dtype": "complex64", +} # odd sign, even filt, complex +par4 = {"nt": 50, "nh": 6, "imag": 0, "dtype": "float32"} # even sign, even filt, real +par4j = { + "nt": 50, + "nh": 6, + "imag": 1j, + "dtype": "complex64", +} # even sign, even filt, complex + np.random.seed(10) -@pytest.mark.parametrize("par", [(par1), (par2)]) +@pytest.mark.parametrize("par", [(par1), (par1j), (par2), (par2j)]) def test_convmtx(par): - """Compare convmtx with np.convolve""" + """Compare convmtx with np.convolve (small filter)""" + x = np.random.normal(0, 1, par["nt"]) + par["imag"] * np.random.normal( + 0, 1, par["nt"] + ) + + h = np.hanning(par["nh"]) + H = convmtx(h, par["nt"], par["nh"] // 2) + + y = np.convolve(x, h, mode="same") + y1 = np.dot(H[: par["nt"]], x) + assert_array_almost_equal(y, y1, decimal=4) + + +@pytest.mark.parametrize("par", [(par1), (par1j), (par2), (par2j)]) +def test_convmtx1(par): + """Compare convmtx with np.convolve (large filter)""" x = np.random.normal(0, 1, par["nt"]) + par["imag"] * np.random.normal( 0, 1, par["nt"] ) - nh = 7 - h = np.hanning(7) - H = convmtx(h, par["nt"]) - H = H[:, nh // 2 : -nh // 2 + 1] + h = np.hanning(par["nh"]) + X = convmtx( + x, par["nh"], par["nh"] // 2 - 1 if par["nh"] % 2 == 0 else par["nh"] // 2 + ) y = np.convolve(x, h, mode="same") - y1 = np.dot(H, x) + y1 = np.dot(X[: par["nt"]], h) assert_array_almost_equal(y, y1, decimal=4) -@pytest.mark.parametrize("par", [(par1), (par2)]) +@pytest.mark.parametrize("par", [(par1), (par1j)]) def test_nonstationary_convmtx(par): """Compare nonstationary_convmtx with convmtx for stationary filter""" x = np.random.normal(0, 1, par["nt"]) + par["imag"] * np.random.normal( 0, 1, par["nt"] ) - nh = 7 - h = np.hanning(7) - H = convmtx(h, par["nt"]) - H = H[:, nh // 2 : -nh // 2 + 1] + h = np.hanning(par["nh"]) + H = convmtx( + h, par["nt"], par["nh"] // 2 - 1 if par["nh"] % 2 == 0 else par["nh"] // 2 + ) H1 = nonstationary_convmtx( np.repeat(h[:, np.newaxis], par["nt"], axis=1).T, par["nt"], - hc=nh // 2, + hc=par["nh"] // 2, pad=(par["nt"], par["nt"]), ) - y = np.dot(H, x) + y = np.dot(H[: par["nt"]], x) y1 = np.dot(H1, x) assert_array_almost_equal(y, y1, decimal=4) From e8ec0f89c9128aa224ca69ce6885f3a182939fc4 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Mon, 7 Aug 2023 22:24:55 +0200 Subject: [PATCH 45/69] fix: correct basic solvers' handling of ndarrays Basic solvers return 1d- or nd-arrays depending on the attribute dims of the operator indipendently on the shape of x0 and data. This change introduces some logic that reshapes the output to dims only if the starting guess was also an nd-array and the operator does not have forceflat=True. --- pylops/linearoperator.py | 2 +- pylops/utils/decorators.py | 6 ++++++ pytests/test_solver.py | 31 ++++++++++++++++++++++++++++++- 3 files changed, 37 insertions(+), 2 deletions(-) diff --git a/pylops/linearoperator.py b/pylops/linearoperator.py index 67ffa59f..0a719cf3 100644 --- a/pylops/linearoperator.py +++ b/pylops/linearoperator.py @@ -97,7 +97,7 @@ class LinearOperator(_LinearOperator): Force an array to be flattened after matvec/rmatvec if the input is ambiguous (i.e., is a 1D array both when operating with ND arrays and with 1D arrays). - Defaults to ``None`` for operators that have no ambiguity or to ``False`` + Defaults to ``None`` for operators that have no ambiguity or to ``True`` for those with ambiguity. name : :obj:`str`, optional .. versionadded:: 2.0.0 diff --git a/pylops/utils/decorators.py b/pylops/utils/decorators.py index f8767b28..887d2d53 100644 --- a/pylops/utils/decorators.py +++ b/pylops/utils/decorators.py @@ -54,10 +54,16 @@ def add_ndarray_support_to_solver(func: Callable) -> Callable: @wraps(func) def wrapper(A, b, *args, **kwargs): # SciPy-type signature + x0flat = False if "x0" in kwargs and kwargs["x0"] is not None: + if kwargs["x0"].ndim == 1: + x0flat = True kwargs["x0"] = kwargs["x0"].ravel() with disabled_ndarray_multiplication(): res = list(func(A, b.ravel(), *args, **kwargs)) + # reshape if x0 was provided unflattened + # (unless the operator has forceflat=True) + if not x0flat and not getattr(A, "forceflat", None): res[0] = res[0].reshape(getattr(A, "dims", (A.shape[1],))) return tuple(res) diff --git a/pytests/test_solver.py b/pytests/test_solver.py index 6030bf99..ba32f79b 100644 --- a/pytests/test_solver.py +++ b/pytests/test_solver.py @@ -94,7 +94,7 @@ def test_cg(par): "par", [(par1), (par2), (par3), (par4), (par1j), (par2j), (par3j), (par3j)] ) def test_cg_ndarray(par): - """CG with linear operator""" + """CG with linear operator (and ndarray as input and output)""" np.random.seed(10) dims = dimsd = (par["nx"], par["ny"]) @@ -119,6 +119,35 @@ def test_cg_ndarray(par): assert_array_almost_equal(x, xinv, decimal=4) +@pytest.mark.parametrize( + "par", [(par1), (par2), (par3), (par4), (par1j), (par2j), (par3j), (par3j)] +) +def test_cg_forceflat(par): + """CG with linear operator (and forced 1darray as input and output)""" + np.random.seed(10) + + dims = dimsd = (par["nx"], par["ny"]) + x = np.ones(dims) + par["imag"] * np.ones(dims) + + A = np.random.normal(0, 10, (x.size, x.size)) + par["imag"] * np.random.normal( + 0, 10, (x.size, x.size) + ) + A = np.conj(A).T @ A # to ensure definite positive matrix + Aop = MatrixMult(A, dtype=par["dtype"], forceflat=True) + Aop.dims = dims + Aop.dimsd = dimsd + + if par["x0"]: + x0 = np.random.normal(0, 10, dims) + par["imag"] * np.random.normal(0, 10, dims) + else: + x0 = None + + y = Aop * x + xinv = cg(Aop, y, x0=x0, niter=2 * x.size, tol=1e-5, show=True)[0] + assert xinv.shape == x.ravel().shape + assert_array_almost_equal(x.ravel(), xinv, decimal=4) + + @pytest.mark.parametrize( "par", [(par1), (par2), (par3), (par4), (par1j), (par2j), (par3j), (par3j)] ) From 390cd98abb7fc0443280029122a6d9ffca79f2c3 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Wed, 23 Aug 2023 16:24:26 -0500 Subject: [PATCH 46/69] bug: fixed convmtx as discussed in https://github.com/PyLops/pylops/issues/517 --- pylops/utils/signalprocessing.py | 29 +++++++++++------------------ 1 file changed, 11 insertions(+), 18 deletions(-) diff --git a/pylops/utils/signalprocessing.py b/pylops/utils/signalprocessing.py index 0aa82e3a..4aa277fa 100644 --- a/pylops/utils/signalprocessing.py +++ b/pylops/utils/signalprocessing.py @@ -24,8 +24,9 @@ def convmtx(h: npt.ArrayLike, n: int, offset: int = 0) -> NDArray: the filter :math:`h` centered on `offset` and the input signal :math:`x`. Equivalent of `MATLAB's convmtx function - `_ for - ``offset=0``. + `_ for: + - ``mode='full'`` when used with ``offset=0``. + - ``mode='same'`` when used with ``offset=len(h)//2`` (after truncating the rows as ``C[:n]``) Parameters ---------- @@ -44,28 +45,20 @@ def convmtx(h: npt.ArrayLike, n: int, offset: int = 0) -> NDArray: """ warnings.warn( "A new implementation of convmtx is provided in v2.2.0 to match " - "MATLAB's convmtx method as stated in the docstring. Prior to v2.2.0," - "The implementation of convmtx provided prior to v2.2.0 was instead " - "not consistent with the documentation. Users are highly encouraged " + "MATLAB's convmtx method as stated in the docstring. The implementation " + "of convmtx provided prior to v2.2.0 was instead not consistent " + "with the documentation. Users are highly encouraged " "to modify their codes accordingly.", FutureWarning, ) ncp = get_array_module(h) nh = len(h) - if nh < n: - h = np.flipud(h) - col_1 = ncp.r_[h[0], ncp.zeros(n + nh - 2, dtype=h.dtype)] - row_1 = ncp.r_[h, ncp.zeros(n - 1, dtype=h.dtype)] - C = get_toeplitz(h)(col_1, row_1) - # apply offset - C = C[:, offset : offset + n] - else: - col_1 = ncp.r_[h, ncp.zeros(n + nh - 2, dtype=h.dtype)] - row_1 = ncp.r_[h[0], ncp.zeros(n - 1, dtype=h.dtype)] - C = get_toeplitz(h)(col_1, row_1) - # apply offset - C = C[offset : offset + nh + n - 1] + col_1 = ncp.r_[h, ncp.zeros(n + nh - 2, dtype=h.dtype)] + row_1 = ncp.r_[h[0], ncp.zeros(n - 1, dtype=h.dtype)] + C = get_toeplitz(h)(col_1, row_1) + # apply offset + C = C[offset : offset + nh + n - 1] return C From 01e2c45b29fd2c25af821ee410c71b6c4f19fad9 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sun, 10 Sep 2023 21:14:28 +0300 Subject: [PATCH 47/69] minor: second implementation to BlendingContinuous A second implementation is provided for the matvec/rmatvec routines of BlendingContinuous, where shifting is performed simultaneously for all sources. This implementation is more favourable when dealing with small number of receivers. --- pylops/waveeqprocessing/blending.py | 109 +++++++++++++++++++++------- 1 file changed, 82 insertions(+), 27 deletions(-) diff --git a/pylops/waveeqprocessing/blending.py b/pylops/waveeqprocessing/blending.py index 70242cb6..ca18ccb3 100644 --- a/pylops/waveeqprocessing/blending.py +++ b/pylops/waveeqprocessing/blending.py @@ -33,8 +33,9 @@ class BlendingContinuous(LinearOperator): Time sampling in seconds times : :obj:`np.ndarray` Absolute ignition times for each source - nproc : :obj:`int`, optional - Number of processors used when applying operator + shiftall : :obj:`bool`, optional + Shift all shots together (``True``) or one at the time (``False``). Defaults to ``shiftall=False`` (original + implementation), however ``shiftall=True`` should be preferred when ``nr`` is small. dtype : :obj:`str`, optional Operator dtype name : :obj:`str`, optional @@ -64,6 +65,7 @@ def __init__( ns: int, dt: float, times: NDArray, + shiftall: bool = False, dtype: DTypeLike = "float64", name: str = "B", ) -> None: @@ -73,40 +75,85 @@ def __init__( self.ns = ns self.dt = dt self.times = times + self.shiftall = shiftall self.nttot = int(np.max(self.times) / self.dt + self.nt + 1) - self.PadOp = Pad((self.nr, self.nt), ((0, 0), (0, 1)), dtype=self.dtype) - # Define shift operators - self.shifts = [] - self.ShiftOps = [] - for i in range(self.ns): - shift = self.times[i] - # This is the part that fits on the grid - shift_int = int(shift // self.dt) - self.shifts.append(shift_int) - # This is the fractional part - diff = (shift / self.dt - shift_int) * self.dt - if diff == 0: - self.ShiftOps.append(None) - else: - self.ShiftOps.append( - Shift( - (self.nr, self.nt + 1), - diff, - axis=1, - sampling=self.dt, - real=False, - dtype=self.dtype, + if not self.shiftall: + # original implementation, where each source is shifted indipendently + self.PadOp = Pad((self.nr, self.nt), ((0, 0), (0, 1)), dtype=self.dtype) + # Define shift operators + self.shifts = [] + self.ShiftOps = [] + for i in range(self.ns): + shift = self.times[i] + # This is the part that fits on the grid + shift_int = int(shift // self.dt) + self.shifts.append(shift_int) + # This is the fractional part + diff = (shift / self.dt - shift_int) * self.dt + if diff == 0: + self.ShiftOps.append(None) + else: + self.ShiftOps.append( + Shift( + (self.nr, self.nt + 1), + diff, + axis=1, + sampling=self.dt, + real=True, + dtype=self.dtype, + ) ) - ) + else: + # alternative implementation, where all sources are shifted at the same time + self.PadOp = Pad( + (self.ns, self.nr, self.nt), ((0, 0), (0, 0), (0, 1)), dtype=self.dtype + ) + # Define shift operator + self.shifts = (times // self.dt).astype(np.int32) + diff = (times / self.dt - self.shifts) * self.dt + diff = np.repeat(diff[:, np.newaxis], self.nr, axis=1) + self.ShiftOp = Shift( + (self.ns, self.nr, self.nt + 1), + diff, + axis=-1, + sampling=self.dt, + real=True, + dtype=self.dtype, + ) + self.diff = diff + super().__init__( dtype=np.dtype(dtype), dims=(self.ns, self.nr, self.nt), dimsd=(self.nr, self.nttot), name=name, ) + self._register_multiplications() + + @reshaped + def _matvec_smallrecs(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + blended_data = ncp.zeros((self.nr, self.nttot), dtype=self.dtype) + shifted_data = self.ShiftOp._matvec(self.PadOp._matvec(x.ravel())).reshape( + self.ns, self.nr, self.nt + 1 + ) + for i, shift_int in enumerate(self.shifts): + blended_data[:, shift_int : shift_int + self.nt + 1] += shifted_data[i] + return blended_data + + @reshaped + def _rmatvec_smallrecs(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + shifted_data = ncp.zeros((self.ns, self.nr, self.nt + 1), dtype=self.dtype) + for i, shift_int in enumerate(self.shifts): + shifted_data[i, :, :] = x[:, shift_int : shift_int + self.nt + 1] + deblended_data = self.PadOp._rmatvec( + self.ShiftOp._rmatvec(shifted_data.ravel()) + ).reshape(self.dims) + return deblended_data @reshaped - def _matvec(self, x: NDArray) -> NDArray: + def _matvec_largerecs(self, x: NDArray) -> NDArray: ncp = get_array_module(x) blended_data = ncp.zeros((self.nr, self.nttot), dtype=self.dtype) for i, shift_int in enumerate(self.shifts): @@ -118,7 +165,7 @@ def _matvec(self, x: NDArray) -> NDArray: return blended_data @reshaped - def _rmatvec(self, x: NDArray) -> NDArray: + def _rmatvec_largerecs(self, x: NDArray) -> NDArray: ncp = get_array_module(x) deblended_data = ncp.zeros((self.ns, self.nr, self.nt), dtype=self.dtype) for i, shift_int in enumerate(self.shifts): @@ -133,6 +180,14 @@ def _rmatvec(self, x: NDArray) -> NDArray: deblended_data[i, :, :] = shifted_data return deblended_data + def _register_multiplications(self) -> None: + if self.shiftall: + self._matvec = self._matvec_smallrecs + self._rmatvec = self._rmatvec_smallrecs + else: + self._matvec = self._matvec_largerecs + self._rmatvec = self._rmatvec_largerecs + def BlendingGroup( nt: int, From e6e0338d36712159fa596fdd35a5405984b7a4d0 Mon Sep 17 00:00:00 2001 From: Matteo Ravasi Date: Sun, 1 Oct 2023 22:03:59 +0300 Subject: [PATCH 48/69] build: switch from setup.py to pyproject.toml (#530) build: switch from setup.py to pyproject.toml --- .github/workflows/build.yaml | 2 +- .../workflows/codacy-coverage-reporter.yaml | 13 +++-- .readthedocs.yaml | 2 +- azure-pipelines.yml | 8 +-- pyproject.toml | 3 ++ setup.cfg | 54 ++++++++++++++++++- setup.py | 54 ------------------- 7 files changed, 71 insertions(+), 65 deletions(-) create mode 100644 pyproject.toml delete mode 100755 setup.py diff --git a/.github/workflows/build.yaml b/.github/workflows/build.yaml index f6852e4c..63b4a4e3 100644 --- a/.github/workflows/build.yaml +++ b/.github/workflows/build.yaml @@ -7,7 +7,7 @@ jobs: strategy: matrix: platform: [ ubuntu-latest, macos-latest ] - python-version: ["3.7", "3.8", "3.9", "3.10"] + python-version: ["3.8", "3.9", "3.10"] runs-on: ${{ matrix.platform }} steps: diff --git a/.github/workflows/codacy-coverage-reporter.yaml b/.github/workflows/codacy-coverage-reporter.yaml index b0a8f90d..5e3ed9fe 100644 --- a/.github/workflows/codacy-coverage-reporter.yaml +++ b/.github/workflows/codacy-coverage-reporter.yaml @@ -6,22 +6,27 @@ on: [push, pull_request_target] jobs: build: + strategy: + matrix: + platform: [ ubuntu-latest, macos-latest ] + python-version: ["3.8", "3.9", "3.10"] - runs-on: ubuntu-latest + runs-on: ${{ matrix.platform }} steps: - uses: actions/checkout@v3 - - name: Set up Python + - name: Set up Python ${{ matrix.python-version }} uses: actions/setup-python@v4 with: - python-version: 3.9 + python-version: ${{ matrix.python-version }} - name: Install dependencies run: | python -m pip install --upgrade pip - pip install flake8 pytest coverage + pip install flake8 pytest if [ -f requirements.txt ]; then pip install -r requirements-dev.txt; fi - name: Install pylops run: | pip install . + pip install coverage - name: Code coverage with coverage run: | coverage run -m pytest diff --git a/.readthedocs.yaml b/.readthedocs.yaml index 146206bb..a74e3480 100644 --- a/.readthedocs.yaml +++ b/.readthedocs.yaml @@ -19,5 +19,5 @@ sphinx: python: install: - requirements: requirements-dev.txt - - method: setuptools + - method: pip path: . diff --git a/azure-pipelines.yml b/azure-pipelines.yml index d778e70a..ef5306bb 100644 --- a/azure-pipelines.yml +++ b/azure-pipelines.yml @@ -55,7 +55,7 @@ jobs: steps: - task: UsePythonVersion@0 inputs: - versionSpec: '3.7' + versionSpec: '3.8' architecture: 'x64' - script: | @@ -65,7 +65,7 @@ jobs: displayName: 'Install prerequisites and library' - script: | - python setup.py test + pytest condition: succeededOrFailed() displayName: 'Run tests' @@ -84,7 +84,7 @@ jobs: steps: - task: UsePythonVersion@0 inputs: - versionSpec: '3.7' + versionSpec: '3.8' architecture: 'x64' - script: | @@ -94,6 +94,6 @@ jobs: displayName: 'Install prerequisites and library' - script: | - python setup.py test + pytest condition: succeededOrFailed() displayName: 'Run tests' diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 00000000..8fd8d67e --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,3 @@ +[build-system] +requires = ["setuptools", "setuptools-scm"] +build-backend = "setuptools.build_meta" diff --git a/setup.cfg b/setup.cfg index 1691b3bc..67636c69 100755 --- a/setup.cfg +++ b/setup.cfg @@ -1,3 +1,56 @@ +[metadata] +name = pylops +fullname = PyLops +description = Python library implementing linear operators to allow solving large-scale optimization problems +long_description = file: README.md +long_description_content_type = text/markdown +author = The PyLops Development Team +author_email = matteoravasi@gmail.com +maintainer = "Matteo Ravasi" +maintainer_email = matteoravasi@gmail.com +license = LGPL-3.0 License +license_file = LICENSE.md +platform = any +keywords = algebra, inverse problems, large-scale optimization +classifiers = + Development Status :: 5 - Production/Stable + Intended Audience :: Developers + Intended Audience :: Science/Research + Intended Audience :: Education + License :: OSI Approved :: GNU Lesser General Public License v3 (LGPLv3), + Natural Language :: English + Operating System :: OS Independent + Programming Language :: Python :: 3 :: Only + Programming Language :: Python :: 3.8 + Programming Language :: Python :: 3.9 + Programming Language :: Python :: 3.10 + Topic :: Scientific/Engineering : Mathematics +url = https://github.com/pylops/pylops + +project_urls = + Documentation = https://pylops.readthedocs.io/ + Release Notes = https://github.com/pylops/pylops/releases + Bug Tracker = https://github.com/pylops/pylops/issues + Source Code = https://github.com/pylops/pylops + +[options] +zip_safe = True +include_package_data = True +packages = find: +python_requires = >=3.8 +install_requires = + numpy >= 1.21.0 + scipy >= 1.4.0 + +[options.extras_require] +advanced = + llvmlite + numba + pyfftw + PyWavelets + scikit-fmm + spgl1 + [aliases] test=pytest @@ -11,7 +64,6 @@ per-file-ignores = __init__.py: F401, F403, F405 max-line-length = 88 - # mypy global options [mypy] plugins = numpy.typing.mypy_plugin diff --git a/setup.py b/setup.py deleted file mode 100755 index 8b82afa2..00000000 --- a/setup.py +++ /dev/null @@ -1,54 +0,0 @@ -import os - -from setuptools import find_packages, setup - - -def src(pth): - return os.path.join(os.path.dirname(__file__), pth) - - -# Project description -descr = ( - "Python library implementing linear operators to allow solving large-scale optimization " - "problems without requiring to explicitly create a dense (or sparse) matrix." -) - -# Setup -setup( - name="pylops", - description=descr, - long_description=open(src("README.md")).read(), - long_description_content_type="text/markdown", - keywords=["algebra", "inverse problems", "large-scale optimization"], - classifiers=[ - "Development Status :: 5 - Production/Stable", - "Intended Audience :: Developers", - "Intended Audience :: Science/Research", - "License :: OSI Approved :: GNU Lesser General Public License v2 or later (LGPLv2+)", - "Natural Language :: English", - "Programming Language :: Python :: 3.8", - "Programming Language :: Python :: 3.9", - "Topic :: Scientific/Engineering :: Mathematics", - ], - author="mrava", - author_email="matteoravasi@gmail.com", - install_requires=["numpy >= 1.21.0", "scipy >= 1.4.0"], - extras_require={ - "advanced": [ - "llvmlite", - "numba", - "pyfftw", - "PyWavelets", - "scikit-fmm", - "spgl1", - ] - }, - packages=find_packages(exclude=["pytests"]), - use_scm_version=dict( - root=".", relative_to=__file__, write_to=src("pylops/version.py") - ), - setup_requires=["pytest-runner", "setuptools_scm"], - test_suite="pytests", - tests_require=["pytest"], - zip_safe=True, -) From 32a40c20ecf450e11937b869a6a3ae3768d3a09e Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sun, 1 Oct 2023 22:08:10 +0300 Subject: [PATCH 49/69] feature: restyling of Kirchhoff operator --- pylops/waveeqprocessing/kirchhoff.py | 607 +++++++-------------------- pytests/test_kirchhoff.py | 18 +- 2 files changed, 149 insertions(+), 476 deletions(-) diff --git a/pylops/waveeqprocessing/kirchhoff.py b/pylops/waveeqprocessing/kirchhoff.py index a743ec6d..3ea6b6f4 100644 --- a/pylops/waveeqprocessing/kirchhoff.py +++ b/pylops/waveeqprocessing/kirchhoff.py @@ -35,10 +35,12 @@ class Kirchhoff(LinearOperator): - r"""Kirchhoff Demigration operator. + r"""Kirchhoff demigration operator. Kirchhoff-based demigration/migration operator. Uses a high-frequency - approximation of Green's function propagators based on ``trav``. + approximation of the Green's function propagators based on traveltimes + and amplitudes that are either computed internally by solving the Eikonal equation, + or passed directly by the user (which can use any other propagation engine of choice). Parameters ---------- @@ -73,23 +75,21 @@ class Kirchhoff(LinearOperator): dynamic : :obj:`bool`, optional .. versionadded:: 2.0.0 - Include dynamic weights in computations (``True``) or not (``False``). - Note that when ``mode=byot``, the user is required to provide such weights - in ``amp``. + Apply dynamic weights in computations (``True``) or not (``False``). This includes both the amplitude + terms of the Green's function and the reflectivity-related scaling term (see equations below). trav : :obj:`numpy.ndarray` or :obj:`tuple`, optional Traveltime table of size :math:`\lbrack (n_y) n_x n_z \times n_s n_r \rbrack` or pair of traveltime tables of size :math:`\lbrack (n_y) n_x n_z \times n_s \rbrack` and :math:`\lbrack (n_y) n_x n_z \times n_r \rbrack` (to be provided if ``mode='byot'``). Note that the latter approach is recommended as less memory demanding - than the former. + than the former. Moreover, only ``mode='dynamic'`` is only possible when traveltimes are provided in + the latter form. amp : :obj:`numpy.ndarray`, optional .. versionadded:: 2.0.0 - Amplitude table of size - :math:`\lbrack (n_y) n_x n_z \times n_s n_r \rbrack` or pair of amplitude tables - of size :math:`\lbrack (n_y) n_x n_z \times n_s \rbrack` and :math:`\lbrack (n_y) n_x n_z \times n_r \rbrack` - (to be provided if ``mode='byot'``). Note that the latter approach is recommended as less memory demanding - than the former. + Pair of amplitude tables of size :math:`\lbrack (n_y) n_x n_z \times n_s \rbrack` and + :math:`\lbrack (n_y) n_x n_z \times n_r \rbrack` (to be provided if ``mode='byot'``). Note that this parameter + is only required when ``mode='dynamic'`` is chosen. aperture : :obj:`float` or :obj:`tuple`, optional .. versionadded:: 2.0.0 @@ -104,18 +104,9 @@ class Kirchhoff(LinearOperator): Maximum allowed angle (either source or receiver side) in degrees. If ``None``, angle aperture limitations are not introduced. See ``aperture`` for implementation details regarding scalar and tuple cases. - - anglerefl : :obj:`np.ndarray`, optional - .. versionadded:: 2.0.0 - - Angle between the normal of the reflectors and the vertical axis in degrees snell : :obj:`float` or :obj:`tuple`, optional - .. versionadded:: 2.0.0 - - Threshold on Snell's law evaluation. If larger, the source-receiver-image - point is discarded. If ``None``, no check on the validity of the Snell's - law is performed. See ``aperture`` for implementation details regarding - scalar and tuple cases. + Deprecated, will be removed in v3.0.0. Simply kept for back-compatibility with previous implementation, + but effectively not affecting the behaviour of the operator. engine : :obj:`str`, optional Engine used for computations (``numpy`` or ``numba``). dtype : :obj:`str`, optional @@ -141,29 +132,34 @@ class Kirchhoff(LinearOperator): Notes ----- The Kirchhoff demigration operator synthesizes seismic data given a - propagation velocity model :math:`v` and a reflectivity model :math:`m`. - In forward mode [1]_, [2]_: + propagation velocity model :math:`v(\mathbf{x})` and a + reflectivity model :math:`m(\mathbf{x})`. In forward mode [1]_, [2]_: .. math:: d(\mathbf{x_r}, \mathbf{x_s}, t) = - \widetilde{w}(t) * \int_V G(\mathbf{x_r}, \mathbf{x}, t) - m(\mathbf{x}) G(\mathbf{x}, \mathbf{x_s}, t)\,\mathrm{d}\mathbf{x} + \widetilde{w}(t) * \int_V \frac{2 cos(\theta)} {v(\mathbf{x})} + G(\mathbf{x_r}, \mathbf{x}, t) G(\mathbf{x}, \mathbf{x_s}, t) + m(\mathbf{x}) \,\mathrm{d}\mathbf{x} + + where :math:`G(\mathbf{x}, \mathbf{x_s}, t)` and :math:`G(\mathbf{x_r}, \mathbf{x}, t)` + are the Green's functions from source-to-subsurface-to-receiver and finally + :math:`\widetilde{w}(t)` is a filtered version of the wavelet :math:`w(t)` + as explained below (or the wavelet itself when ``wavfilter=False``). + Moreover, an angle scaling is included in the modelling operator, + where :math:`\theta=(\theta_s-\theta_r)/2` is half of the opening angle, + where :math:`\theta_s` and :math:`\theta_r` are the angles between the source-side + and receiver-side rays and the vertical at the image point. - where :math:`m(\mathbf{x})` represents the reflectivity - at every location in the subsurface, :math:`G(\mathbf{x}, \mathbf{x_s}, t)` - and :math:`G(\mathbf{x_r}, \mathbf{x}, t)` are the Green's functions - from source-to-subsurface-to-receiver and finally :math:`\widetilde{w}(t)` is - a filtered version of the wavelet :math:`w(t)` [3]_ (or the wavelet itself when - ``wavfilter=False``). In our implementation, the following high-frequency - approximation of the Green's functions is adopted: + In our implementation, the following high-frequency approximation of the + Green's functions is adopted: .. math:: G(\mathbf{x_r}, \mathbf{x}, \omega) = a(\mathbf{x_r}, \mathbf{x}) e^{j \omega t(\mathbf{x_r}, \mathbf{x})} where :math:`a(\mathbf{x_r}, \mathbf{x})` is the amplitude and - :math:`t(\mathbf{x_r}, \mathbf{x})` is the traveltime. When ``dynamic=False`` the - amplitude is disregarded leading to a kinematic-only Kirchhoff operator. + :math:`t(\mathbf{x_r}, \mathbf{x})` is the traveltime. When ``dynamic=False``, the + amplitude correction terms are disregarded leading to a kinematic-only Kirchhoff operator. .. math:: d(\mathbf{x_r}, \mathbf{x_s}, t) = @@ -171,36 +167,32 @@ class Kirchhoff(LinearOperator): t(\mathbf{x}, \mathbf{x_s}))} m(\mathbf{x}) \,\mathrm{d}\mathbf{x} On the other hand, when ``dynamic=True``, the amplitude scaling is defined as - :math:`a(\mathbf{x}, \mathbf{y})=\frac{1}{\|\mathbf{x} - \mathbf{y}\|}`, - that is, the reciprocal of the distance between the two points, - approximating the geometrical spreading of the wavefront. - Moreover an angle scaling is included in the modelling operator - added as follows: - .. math:: - d(\mathbf{x_r}, \mathbf{x_s}, t) = - \tilde{w}(t) * \int_V a(\mathbf{x}, \mathbf{x_s}) a(\mathbf{x}, \mathbf{x_r}) - \frac{|cos \theta_s + cos \theta_r|} {v(\mathbf{x})} e^{j \omega (t(\mathbf{x_r}, \mathbf{x}) + - t(\mathbf{x}, \mathbf{x_s}))} m(\mathbf{x}) \,\mathrm{d}\mathbf{x} + * ``2D``: :math:`a(\mathbf{x}, \mathbf{y})=\frac{1}{\sqrt{dist(\mathbf{x}, \mathbf{y})}}` + * ``3D``: :math:`a(\mathbf{x}, \mathbf{y})=\frac{1}{dist(\mathbf{x}, \mathbf{y})}` - where :math:`\theta_s` and :math:`\theta_r` are the angles between the source-side - and receiver-side rays and the normal to the reflector at the image point (or - the vertical axis at the image point when ``reflslope=None``), respectively. + approximating the geometrical spreading of the wavefront. For ``mode=analytic``, + :math:`dist(\mathbf{x}, \mathbf{y})=\|\mathbf{x} - \mathbf{y}\|`, whilst for + ``mode=eikonal``, this is instead computed internally by the eikonal solver. + + The wavelet filtering is applied as follows: + + * ``2D``: :math:`\tilde{W}(f)=\sqrt{j\omega} \cdot W(f)` + * ``3D``: :math:`\tilde{W}(f)=j\omega \cdot W(f)` Depending on the choice of ``mode`` the traveltime and amplitude of the Green's function will be also computed differently: - * ``mode=analytic`` or ``mode=eikonal``: traveltimes, geometrical spreading, and angles - are computed for every source-image point-receiver triplets and the - Green's functions are implemented from traveltime look-up tables, placing - scaled reflectivity values at corresponding source-to-receiver time in the - data. - * ``byot``: bring your own tables. Traveltime table are provided + * ``mode=analytic`` or ``mode=eikonal``: traveltimes, amplitudes, and angles + are computed for every source-image point-receiver triplets upfront and the + Green's functions are implemented from traveltime and amplitude look-up tables, + placing scaled reflectivity values at corresponding source-to-receiver time + in the data. + * ``mode=byot``: bring your own tables. Traveltime table are provided directly by user using ``trav`` input parameter. Similarly, in this case one - can provide their own amplitude scaling ``amp`` (which should include the angle - scaling too). + can also provide their own amplitude scaling ``amp``. - Three aperture limitations have been also implemented as defined by: + Two aperture limitations have been also implemented as defined by: * ``aperture``: the maximum allowed aperture is expressed as the ratio of offset over depth. This aperture limitation avoid including grazing angles @@ -208,15 +200,10 @@ class Kirchhoff(LinearOperator): edges of the aperture; * ``angleaperture``: the maximum allowed angle aperture is expressed as the difference between the incident or emerging angle at every image point and - the vertical axis (or the normal to the reflector if ``anglerefl`` is provided. - This aperture limitation also avoid including grazing angles whose contributions - can introduce aliasing effects. Note that for a homogenous medium and slowly varying - heterogenous medium the offset and angle aperture limits may work in the same way; - * ``snell``: the maximum allowed snell's angle is expressed as the absolute value of - the sum between incident and emerging angles defined as in the ``angleaperture`` case. - This aperture limitation is introduced to turn a scattering-based Kirchhoff engine into - a reflection-based Kirchhoff engine where each image point is not considered as scatter - but as a local horizontal (or straight) reflector. + the vertical axis. This aperture limitation also avoid including grazing angles + whose contributions can introduce aliasing effects. Note that for a homogenous + medium and slowly varying heterogeneous medium the offset and angle aperture limits + may work in the same way. Finally, the adjoint of the demigration operator is a *migration* operator which projects data in the model domain creating an image of the subsurface @@ -252,7 +239,6 @@ def __init__( amp: Optional[NDArray] = None, aperture: Optional[Tuple[float, float]] = None, angleaperture: Union[float, Tuple[float, float]] = 90.0, - anglerefl: Optional[NDArray] = None, snell: Optional[Tuple[float, float]] = None, engine: str = "numpy", dtype: DTypeLike = "float64", @@ -288,17 +274,20 @@ def __init__( self.dt = t[1] - t[0] self.nt = len(t) - # store ix-iy locations of sources and receivers - dx = x[1] - x[0] - if self.ndims == 2: - self.six = np.tile((srcs[0] - x[0]) // dx, (nr, 1)).T.astype(int).ravel() - self.rix = np.tile((recs[0] - x[0]) // dx, (ns, 1)).astype(int).ravel() - elif self.ndims == 3: - # TODO: 3D normalized distances - pass - - # compute traveltime + # store ix-iy locations of sources and receivers for aperture filter self.dynamic = dynamic + if self.dynamic: + dx = x[1] - x[0] + if self.ndims == 2: + self.six = ( + np.tile((srcs[0] - x[0]) // dx, (nr, 1)).T.astype(int).ravel() + ) + self.rix = np.tile((recs[0] - x[0]) // dx, (ns, 1)).astype(int).ravel() + elif self.ndims == 3: + # TODO: 3D normalized distances + raise NotImplementedError("dynamic=True currently not available in 3D") + + # compute traveltime and distances self.travsrcrec = True # use separate tables for src and rec traveltimes if mode in ["analytic", "eikonal", "byot"]: if mode in ["analytic", "eikonal"]: @@ -306,63 +295,72 @@ def __init__( ( self.trav_srcs, self.trav_recs, - self.dist_srcs, - self.dist_recs, + dist_srcs, + dist_recs, trav_srcs_grad, trav_recs_grad, ) = Kirchhoff._traveltime_table(z, x, srcs, recs, vel, y=y, mode=mode) if self.dynamic: # need to add a scalar in the denominator of amplitude term to avoid - # division by 0, currently set to 1/100 of max distance to avoid having - # to large scaling around the source. This number may change in future + # division by 0, currently set to 1e-4 of max distance to avoid having + # too large scaling around the source. This number may change in future # or left to the user to define epsdist = 1e-2 - self.maxdist = epsdist * ( - np.max(self.dist_srcs) + np.max(self.dist_recs) - ) - # compute angles + self.maxdist = epsdist * (np.max(dist_srcs) + np.max(dist_recs)) if self.ndims == 2: - # 2d with vertical - if anglerefl is None: - self.angle_srcs = np.arctan2( - trav_srcs_grad[0], trav_srcs_grad[1] - ).reshape(np.prod(dims), ns) - self.angle_recs = np.arctan2( - trav_recs_grad[0], trav_recs_grad[1] - ).reshape(np.prod(dims), nr) - self.cosangle_srcs = np.cos(self.angle_srcs) - self.cosangle_recs = np.cos(self.angle_recs) - else: - # TODO: 2D with normal - raise NotImplementedError( - "angle scaling with anglerefl currently not available" - ) + self.amp_srcs, self.amp_recs = 1.0 / np.sqrt( + dist_srcs + self.maxdist + ), 1.0 / np.sqrt(dist_recs + self.maxdist) else: - # TODO: 3D - raise NotImplementedError( - "dynamic=True currently not available in 3D" - ) + self.amp_srcs, self.amp_recs = 1.0 / ( + dist_srcs + self.maxdist + ), 1.0 / (dist_recs + self.maxdist) else: if isinstance(trav, tuple): self.trav_srcs, self.trav_recs = trav else: self.travsrcrec = False self.trav = trav + + if self.dynamic and not self.travsrcrec: + raise NotImplementedError( + "separate traveltime tables must be provided " + "when selecting mode=dynamic" + ) if self.dynamic: if isinstance(amp, tuple): self.amp_srcs, self.amp_recs = amp else: - self.amp = amp - # in byot mode, angleaperture and snell checks are not performed - self.angle_srcs = np.ones( - (self.ny * self.nx * self.nz, ns), dtype=dtype - ) - self.angle_recs = np.ones( - (self.ny * self.nx * self.nz, nr), dtype=dtype - ) + raise NotImplementedError( + "separate amplitude tables must be provided " + ) + + if self.travsrcrec: + # compute traveltime gradients at image points + trav_srcs_grad = np.gradient( + self.trav_srcs.reshape(*dims, ns), + axis=np.arange(self.ndims), + ) + trav_recs_grad = np.gradient( + self.trav_recs.reshape(*dims, nr), + axis=np.arange(self.ndims), + ) else: raise NotImplementedError("method must be analytic, eikonal or byot") + # compute angles with vertical + if self.dynamic: + if self.ndims == 2: + self.angle_srcs = np.arctan2( + trav_srcs_grad[0], trav_srcs_grad[1] + ).reshape(np.prod(dims), ns) + self.angle_recs = np.arctan2( + trav_recs_grad[0], trav_recs_grad[1] + ).reshape(np.prod(dims), nr) + else: + # TODO: 3D + raise NotImplementedError("dynamic=True currently not available in 3D") + # pre-compute traveltime indices if total traveltime is used if not self.travsrcrec: self.itrav = (self.trav / self.dt).astype("int32") @@ -383,31 +381,25 @@ def __init__( self.aperturetap = taper(41, 20, "hanning")[20:] # define aperture + # if aperture=None, we want to ensure the check is always matched (no aperture limits...) if aperture is not None: warnings.warn( "Aperture is currently defined as ratio of offset over depth, " - "and may be not ideal for highly heterogenous media" + "and may be not ideal for highly heterogeneous media" ) self.aperture = ( - (2 * self.nx / self.nz,) - if aperture is None - else _value_or_sized_to_array(aperture) + (self.nx + 1,) if aperture is None else _value_or_sized_to_array(aperture) ) if len(self.aperture) == 1: self.aperture = np.array([0.8 * self.aperture[0], self.aperture[0]]) - # define angle aperture and snell law + # define angle aperture angleaperture = [0.0, 1000.0] if angleaperture is None else angleaperture self.angleaperture = np.deg2rad(_value_or_sized_to_array(angleaperture)) if len(self.angleaperture) == 1: self.angleaperture = np.array( [0.8 * self.angleaperture[0], self.angleaperture[0]] ) - self.snell = ( - (np.pi,) if snell is None else np.deg2rad(_value_or_sized_to_array(snell)) - ) - if len(self.snell) == 1: - self.snell = np.array([0.8 * self.snell[0], self.snell[0]]) # dimensions self.ns, self.nr = ns, nr @@ -635,18 +627,19 @@ def _wavelet_reshaping( dimrec: int, dimv: int, ) -> NDArray: - """Apply wavelet reshaping as from theory in [1]_""" + """Apply wavelet reshaping to account for omega scaling factor + originating from the wave equation""" f = np.fft.rfftfreq(len(wav), dt) W = np.fft.rfft(wav, n=len(wav)) if dimsrc == 2 and dimv == 2: # 2D - Wfilt = W * (2 * np.pi * f) + Wfilt = W * np.sqrt(1j * 2 * np.pi * f) elif (dimsrc == 2 or dimrec == 2) and dimv == 3: # 2.5D raise NotImplementedError("2.D wavelet currently not available") elif dimsrc == 3 and dimrec == 3 and dimv == 3: # 3D - Wfilt = W * (-1j * 2 * np.pi * f) + Wfilt = W * (1j * 2 * np.pi * f) wavfilt = np.fft.irfft(Wfilt, n=len(wav)) return wavfilt @@ -750,225 +743,6 @@ def _travsrcrec_kirch_rmatvec( ) return y - @staticmethod - def _amp_kirch_matvec( - x: NDArray, - y: NDArray, - nsnr: int, - nt: int, - ni: int, - itrav: NDArray, - travd: NDArray, - amp: NDArray, - aperturemin: float, - aperturemax: float, - aperturetap: NDArray, - nz: int, - six: NDArray, - rix: NDArray, - angleaperturemin: float, - angleaperturemax: float, - angles_srcs: NDArray, - angles_recs: NDArray, - snellmin: float, - snellmax: float, - ) -> NDArray: - nr = angles_recs.shape[-1] - daperture = aperturemax - aperturemin - dangleaperture = angleaperturemax - angleaperturemin - dsnell = snellmax - snellmin - for isrcrec in prange(nsnr): - # extract traveltime, amplitude, src/rec coordinates at given src/pair - itravisrcrec = itrav[:, isrcrec] - travdisrcrec = travd[:, isrcrec] - ampisrcrec = amp[:, isrcrec] - sixisrcrec = six[isrcrec] - rixisrcrec = rix[isrcrec] - # extract source and receiver angles - angles_src = angles_srcs[:, isrcrec // nr] - angles_rec = angles_recs[:, isrcrec % nr] - for ii in range(ni): - # extract traveltime, amplitude at given image point - itravii = itravisrcrec[ii] - travdii = travdisrcrec[ii] - damp = ampisrcrec[ii] - # extract source and receiver angle - angle_src = angles_src[ii] - angle_rec = angles_rec[ii] - abs_angle_src = abs(angle_src) - abs_angle_rec = abs(angle_rec) - abs_angle_src_rec = abs(angle_src + angle_rec) - aptap = 1.0 - # angle apertures checks - if ( - abs_angle_src < angleaperturemax - and abs_angle_rec < angleaperturemax - and abs_angle_src_rec < snellmax - ): - if abs_angle_src >= angleaperturemin: - # extract source angle aperture taper value - aptap = ( - aptap - * aperturetap[ - int( - 20 - * (abs_angle_src - angleaperturemin) - // dangleaperture - ) - ] - ) - if abs_angle_rec >= angleaperturemin: - # extract receiver angle aperture taper value - aptap = ( - aptap - * aperturetap[ - int( - 20 - * (abs_angle_rec - angleaperturemin) - // dangleaperture - ) - ] - ) - if abs_angle_src_rec >= snellmin: - # extract snell taper value - aptap = ( - aptap - * aperturetap[ - int(20 * (abs_angle_src_rec - snellmin) // dsnell) - ] - ) - - # identify x-index of image point - iz = ii % nz - # aperture check - aperture = abs(sixisrcrec - rixisrcrec) / iz - if aperture < aperturemax: - if aperture >= aperturemin: - # extract aperture taper value - aptap = ( - aptap - * aperturetap[ - int(20 * ((aperture - aperturemin) // daperture)) - ] - ) - # time limit check - if 0 <= itravii < nt - 1: - # assign values - y[isrcrec, itravii] += x[ii] * (1 - travdii) * damp * aptap - y[isrcrec, itravii + 1] += x[ii] * travdii * damp * aptap - return y - - @staticmethod - def _amp_kirch_rmatvec( - x: NDArray, - y: NDArray, - nsnr: int, - nt: int, - ni: int, - itrav: NDArray, - travd: NDArray, - amp: NDArray, - aperturemin: float, - aperturemax: float, - aperturetap: NDArray, - nz: int, - six: NDArray, - rix: NDArray, - angleaperturemin: float, - angleaperturemax: float, - angles_srcs: NDArray, - angles_recs: NDArray, - snellmin: float, - snellmax: float, - ) -> NDArray: - nr = angles_recs.shape[-1] - daperture = aperturemax - aperturemin - dangleaperture = angleaperturemax - angleaperturemin - dsnell = snellmax - snellmin - for ii in prange(ni): - itravii = itrav[ii] - travdii = travd[ii] - ampii = amp[ii] - # extract source and receiver angles - angle_srcs = angles_srcs[ii] - angle_recs = angles_recs[ii] - # identify x-index of image point - iz = ii % nz - for isrcrec in range(nsnr): - itravisrcrecii = itravii[isrcrec] - travdisrcrecii = travdii[isrcrec] - sixisrcrec = six[isrcrec] - rixisrcrec = rix[isrcrec] - # extract source and receiver angle - angle_src = angle_srcs[isrcrec // nr] - angle_rec = angle_recs[isrcrec % nr] - abs_angle_src = abs(angle_src) - abs_angle_rec = abs(angle_rec) - abs_angle_src_rec = abs(angle_src + angle_rec) - aptap = 1.0 - # angle apertures checks - if ( - abs_angle_src < angleaperturemax - and abs_angle_rec < angleaperturemax - and abs_angle_src_rec < snellmax - ): - if abs_angle_src >= angleaperturemin: - # extract source angle aperture taper value - aptap = ( - aptap - * aperturetap[ - int( - 20 - * (abs_angle_src - angleaperturemin) - // dangleaperture - ) - ] - ) - if abs_angle_rec >= angleaperturemin: - # extract receiver angle aperture taper value - aptap = ( - aptap - * aperturetap[ - int( - 20 - * (abs_angle_rec - angleaperturemin) - // dangleaperture - ) - ] - ) - if abs_angle_src_rec >= snellmin: - # extract snell taper value - aptap = ( - aptap - * aperturetap[ - int(20 * (abs_angle_src_rec - snellmin) // dsnell) - ] - ) - - # aperture check - aperture = abs(sixisrcrec - rixisrcrec) / iz - if aperture < aperturemax: - if aperture >= aperturemin: - # extract aperture taper value - aptap = ( - aptap - * aperturetap[ - int(20 * ((aperture - aperturemin) // daperture)) - ] - ) - # time limit check - if 0 <= itravisrcrecii < nt - 1: - # assign values - y[ii] += ( - ( - x[isrcrec, itravisrcrecii] * (1 - travdisrcrecii) - + x[isrcrec, itravisrcrecii + 1] * travdisrcrecii - ) - * ampii[isrcrec] - * aptap - ) - return y - @staticmethod def _ampsrcrec_kirch_matvec( x: NDArray, @@ -981,8 +755,8 @@ def _ampsrcrec_kirch_matvec( vel: NDArray, trav_srcs: NDArray, trav_recs: NDArray, - dist_srcs: NDArray, - dist_recs: NDArray, + amp_srcs: NDArray, + amp_recs: NDArray, aperturemin: float, aperturemax: float, aperturetap: NDArray, @@ -993,44 +767,39 @@ def _ampsrcrec_kirch_matvec( angleaperturemax: float, angles_srcs: NDArray, angles_recs: NDArray, - snellmin: float, - snellmax: float, - maxdist: float, ) -> NDArray: daperture = aperturemax - aperturemin dangleaperture = angleaperturemax - angleaperturemin - dsnell = snellmax - snellmin for isrc in prange(ns): travisrc = trav_srcs[:, isrc] - distisrc = dist_srcs[:, isrc] + ampisrc = amp_srcs[:, isrc] angleisrc = angles_srcs[:, isrc] for irec in range(nr): travirec = trav_recs[:, irec] trav = travisrc + travirec itrav = (trav / dt).astype("int32") travd = trav / dt - itrav - distirec = dist_recs[:, irec] + ampirec = amp_recs[:, irec] angleirec = angles_recs[:, irec] - dist = distisrc + distirec - amp = np.abs(np.cos(angleisrc) + np.cos(angleirec)) / (dist + maxdist) sixisrcrec = six[isrc * nr + irec] rixisrcrec = rix[isrc * nr + irec] + # compute cosine of half opening angle and total amplitude scaling + cosangle = np.cos((angleisrc - angleirec) / 2.0) + amp = 2.0 * cosangle * ampisrc * ampirec / vel for ii in range(ni): itravii = itrav[ii] travdii = travd[ii] - damp = amp[ii] / vel[ii] + damp = amp[ii] # extract source and receiver angle at given image point angle_src = angleisrc[ii] angle_rec = angleirec[ii] abs_angle_src = abs(angle_src) abs_angle_rec = abs(angle_rec) - abs_angle_src_rec = abs(angle_src + angle_rec) - aptap = 1.0 # angle apertures checks + aptap = 1.0 if ( abs_angle_src < angleaperturemax and abs_angle_rec < angleaperturemax - and abs_angle_src_rec < snellmax ): if abs_angle_src >= angleaperturemin: # extract source angle aperture taper value @@ -1056,19 +825,11 @@ def _ampsrcrec_kirch_matvec( ) ] ) - if abs_angle_src_rec >= snellmin: - # extract snell taper value - aptap = ( - aptap - * aperturetap[ - int(20 * (abs_angle_src_rec - snellmin) // dsnell) - ] - ) # identify x-index of image point iz = ii % nz # aperture check - aperture = abs(sixisrcrec - rixisrcrec) / iz + aperture = abs(sixisrcrec - rixisrcrec) / (iz + 1) if aperture < aperturemax: if aperture >= aperturemin: # extract aperture taper value @@ -1102,8 +863,8 @@ def _ampsrcrec_kirch_rmatvec( vel: NDArray, trav_srcs: NDArray, trav_recs: NDArray, - dist_srcs: NDArray, - dist_recs: NDArray, + amp_srcs: NDArray, + amp_recs: NDArray, aperturemin: float, aperturemax: float, aperturetap: NDArray, @@ -1114,18 +875,14 @@ def _ampsrcrec_kirch_rmatvec( angleaperturemax: float, angles_srcs: NDArray, angles_recs: NDArray, - snellmin: float, - snellmax: float, - maxdist: float, ) -> NDArray: daperture = aperturemax - aperturemin dangleaperture = angleaperturemax - angleaperturemin - dsnell = snellmax - snellmin for ii in prange(ni): trav_srcsii = trav_srcs[ii] trav_recsii = trav_recs[ii] - dist_srcsii = dist_srcs[ii] - dist_recsii = dist_recs[ii] + amp_srcsii = amp_srcs[ii] + amp_recsii = amp_recs[ii] velii = vel[ii] angle_srcsii = angles_srcs[ii] angle_recsii = angles_recs[ii] @@ -1133,31 +890,27 @@ def _ampsrcrec_kirch_rmatvec( iz = ii % nz for isrc in range(ns): trav_srcii = trav_srcsii[isrc] - dist_srcii = dist_srcsii[isrc] + amp_srcii = amp_srcsii[isrc] angle_src = angle_srcsii[isrc] for irec in range(nr): trav_recii = trav_recsii[irec] travii = trav_srcii + trav_recii itravii = int(travii / dt) travdii = travii / dt - itravii - dist_recii = dist_recsii[irec] + amp_recii = amp_recsii[irec] angle_rec = angle_recsii[irec] - dist = dist_srcii + dist_recii - ampii = np.abs(np.cos(angle_src) + np.cos(angle_rec)) / ( - dist + maxdist - ) - damp = ampii / velii sixisrcrec = six[isrc * nr + irec] rixisrcrec = rix[isrc * nr + irec] abs_angle_src = abs(angle_src) abs_angle_rec = abs(angle_rec) - abs_angle_src_rec = abs(angle_src + angle_rec) - aptap = 1.0 + # compute cosine of half opening angle and total amplitude scaling + cosangle = np.cos((angle_src - angle_rec) / 2.0) + damp = 2.0 * cosangle * amp_srcii * amp_recii / velii # angle apertures checks + aptap = 1.0 if ( abs_angle_src < angleaperturemax and abs_angle_rec < angleaperturemax - and abs_angle_src_rec < snellmax ): if abs_angle_src >= angleaperturemin: # extract source angle aperture taper value @@ -1183,17 +936,9 @@ def _ampsrcrec_kirch_rmatvec( ) ] ) - if abs_angle_src_rec >= snellmin: - # extract snell taper value - aptap = ( - aptap - * aperturetap[ - int(20 * (abs_angle_src_rec - snellmin) // dsnell) - ] - ) # aperture check - aperture = abs(sixisrcrec - rixisrcrec) / iz + aperture = abs(sixisrcrec - rixisrcrec) / (iz + 1) if aperture < aperturemax: if aperture >= aperturemin: # extract aperture taper value @@ -1229,9 +974,6 @@ def _register_multiplications(self, engine: str) -> None: if self.dynamic and self.travsrcrec: self._kirch_matvec = jit(**numba_opts)(self._ampsrcrec_kirch_matvec) self._kirch_rmatvec = jit(**numba_opts)(self._ampsrcrec_kirch_rmatvec) - elif self.dynamic and not self.travsrcrec: - self._kirch_matvec = jit(**numba_opts)(self._amp_kirch_matvec) - self._kirch_rmatvec = jit(**numba_opts)(self._amp_kirch_rmatvec) elif self.travsrcrec: self._kirch_matvec = jit(**numba_opts)(self._travsrcrec_kirch_matvec) self._kirch_rmatvec = jit(**numba_opts)(self._travsrcrec_kirch_rmatvec) @@ -1245,9 +987,6 @@ def _register_multiplications(self, engine: str) -> None: if self.dynamic and self.travsrcrec: self._kirch_matvec = self._ampsrcrec_kirch_matvec self._kirch_rmatvec = self._ampsrcrec_kirch_rmatvec - elif self.dynamic and not self.travsrcrec: - self._kirch_matvec = self._amp_kirch_matvec - self._kirch_rmatvec = self._amp_kirch_rmatvec elif self.travsrcrec: self._kirch_matvec = self._travsrcrec_kirch_matvec self._kirch_rmatvec = self._travsrcrec_kirch_rmatvec @@ -1270,32 +1009,8 @@ def _matvec(self, x: NDArray) -> NDArray: self.vel, self.trav_srcs, self.trav_recs, - self.dist_srcs, - self.dist_recs, - self.aperture[0], - self.aperture[1], - self.aperturetap, - self.nz, - self.six, - self.rix, - self.angleaperture[0], - self.angleaperture[1], - self.angle_srcs, - self.angle_recs, - self.snell[0], - self.snell[1], - self.maxdist, - ) - elif self.dynamic and not self.travsrcrec: - inputs = ( - x.ravel(), - y, - self.nsnr, - self.nt, - self.ni, - self.itrav, - self.travd, - self.amp, + self.amp_srcs, + self.amp_recs, self.aperture[0], self.aperture[1], self.aperturetap, @@ -1306,10 +1021,8 @@ def _matvec(self, x: NDArray) -> NDArray: self.angleaperture[1], self.angle_srcs, self.angle_recs, - self.snell[0], - self.snell[1], ) - elif not self.dynamic and self.travsrcrec: + elif self.travsrcrec: inputs = ( x.ravel(), y, @@ -1321,7 +1034,7 @@ def _matvec(self, x: NDArray) -> NDArray: self.trav_srcs, self.trav_recs, ) - elif not self.dynamic and not self.travsrcrec: + elif not self.travsrcrec: inputs = (x.ravel(), y, self.nsnr, self.nt, self.ni, self.itrav, self.travd) y = self._kirch_matvec(*inputs) @@ -1345,32 +1058,8 @@ def _rmatvec(self, x: NDArray) -> NDArray: self.vel, self.trav_srcs, self.trav_recs, - self.dist_srcs, - self.dist_recs, - self.aperture[0], - self.aperture[1], - self.aperturetap, - self.nz, - self.six, - self.rix, - self.angleaperture[0], - self.angleaperture[1], - self.angle_srcs, - self.angle_recs, - self.snell[0], - self.snell[1], - self.maxdist, - ) - elif self.dynamic and not self.travsrcrec: - inputs = ( - x, - y, - self.nsnr, - self.nt, - self.ni, - self.itrav, - self.travd, - self.amp, + self.amp_srcs, + self.amp_recs, self.aperture[0], self.aperture[1], self.aperturetap, @@ -1381,10 +1070,8 @@ def _rmatvec(self, x: NDArray) -> NDArray: self.angleaperture[1], self.angle_srcs, self.angle_recs, - self.snell[0], - self.snell[1], ) - elif not self.dynamic and self.travsrcrec: + elif self.travsrcrec: inputs = ( x, y, @@ -1396,7 +1083,7 @@ def _rmatvec(self, x: NDArray) -> NDArray: self.trav_srcs, self.trav_recs, ) - elif not self.dynamic and not self.travsrcrec: + elif not self.travsrcrec: inputs = (x, y, self.nsnr, self.nt, self.ni, self.itrav, self.travd) y = self._kirch_rmatvec(*inputs) diff --git a/pytests/test_kirchhoff.py b/pytests/test_kirchhoff.py index 317821c3..b210a95b 100755 --- a/pytests/test_kirchhoff.py +++ b/pytests/test_kirchhoff.py @@ -287,7 +287,7 @@ def test_kirchhoff3d(par): ) def test_kirchhoff2d_trav_vs_travsrcrec(par): """Compare 2D Kirchhoff operator forward and adjoint when using trav (original behavior) - or trav_src and trav_rec (new reccomended behaviour)""" + or trav_src and trav_rec (new reccommended behaviour)""" # new behaviour Dop = Kirchhoff( @@ -311,21 +311,7 @@ def test_kirchhoff2d_trav_vs_travsrcrec(par): ) + Dop.trav_recs.reshape(PAR["nx"] * PAR["nz"], 1, PAR["nrx"]) trav = trav.reshape(PAR["nx"] * PAR["nz"], PAR["nsx"] * PAR["nrx"]) if par["dynamic"]: - dist = Dop.dist_srcs.reshape( - PAR["nx"] * PAR["nz"], PAR["nsx"], 1 - ) + Dop.dist_recs.reshape(PAR["nx"] * PAR["nz"], 1, PAR["nrx"]) - dist = dist.reshape(PAR["nx"] * PAR["nz"], PAR["nsx"] * PAR["nrx"]) - - cosangle = np.cos(Dop.angle_srcs).reshape( - PAR["nx"] * PAR["nz"], PAR["nsx"], 1 - ) + np.cos(Dop.angle_recs).reshape(PAR["nx"] * PAR["nz"], 1, PAR["nrx"]) - cosangle = cosangle.reshape(PAR["nx"] * PAR["nz"], PAR["nsx"] * PAR["nrx"]) - - epsdist = 1e-2 - amp = 1 / (dist + epsdist * np.max(dist)) - - amp *= np.abs(cosangle) - amp /= v0 + amp = (Dop.amp_srcs, Dop.amp_recs) D1op = Kirchhoff( z, From 7ed1d7f657ec37d602eb65d0df8d0fcaa7d51fbb Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sun, 1 Oct 2023 22:10:17 +0300 Subject: [PATCH 50/69] minor: coverage action to use only python 3.8 --- .github/workflows/codacy-coverage-reporter.yaml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/codacy-coverage-reporter.yaml b/.github/workflows/codacy-coverage-reporter.yaml index 5e3ed9fe..b16411b4 100644 --- a/.github/workflows/codacy-coverage-reporter.yaml +++ b/.github/workflows/codacy-coverage-reporter.yaml @@ -8,8 +8,8 @@ jobs: build: strategy: matrix: - platform: [ ubuntu-latest, macos-latest ] - python-version: ["3.8", "3.9", "3.10"] + platform: [ ubuntu-latest, ] + python-version: ["3.8", ] runs-on: ${{ matrix.platform }} steps: From 937a24b92e78b2f8790136de7e83f0bda6d6c4f2 Mon Sep 17 00:00:00 2001 From: fedorgoncharov Date: Tue, 3 Oct 2023 11:56:02 +0200 Subject: [PATCH 51/69] remark on classical sinograms added + jacobian correction; added two plots - data with and without jacobian correction; added more explanations on implementation of custom curve --- tutorials/ctscan.py | 42 +++++++++++++++++++++++++++++++++++++----- 1 file changed, 37 insertions(+), 5 deletions(-) diff --git a/tutorials/ctscan.py b/tutorials/ctscan.py index 7869593d..1374de58 100755 --- a/tutorials/ctscan.py +++ b/tutorials/ctscan.py @@ -28,11 +28,17 @@ # such a type: # # .. math:: -# t(r,\theta; x) = \tan(90°-\theta)x + \frac{r}{\sin(\theta)} +# t\sin(\theta) + x\cos(\theta) = r, # # where :math:`\theta` is the angle between the x-axis (:math:`x`) and # the perpendicular to the summation line and :math:`r` is the distance -# from the origin of the summation line. +# from the origin of the summation line. Radon transform in CT +# corresponds to the integral of the input image along the straight line above. +# To implement the integration in PyLops we simply need to express +# :math:`t(r,\theta;x)` which is given by: +# +# .. math:: +# t(r,\theta; x) = \tan\left(\frac{\pi}{2}-\theta\right)x + \frac{r}{\sin(\theta)}. @jit(nopython=True) @@ -43,6 +49,12 @@ def radoncurve(x, r, theta): + ny // 2 ) +############################################################################### +# Note that in the above implementation we added centering :math:`t \mapsto t - n_y/2` and +# :math:`r \mapsto r - n_y/2` so that origin of integration lines is exactly in the +# center of the image (centering for :math:`x` is not needed because we will use +# ``centeredh=True`` in the constructor of ``Radon2D``). + x = np.load("../testdata/optimization/shepp_logan_phantom.npy").T x = x / x.max() @@ -81,14 +93,34 @@ def radoncurve(x, r, theta): axs[1].set_title("Data") axs[1].axis("tight") axs[2].imshow(xrec.T, cmap="gray") -axs[2].set_title("Adjoint model") +axs[2].set_title("Adjoint model in PyLops") axs[2].axis("tight") fig.tight_layout() +############################################################################### +# Note that our raw data ``y`` *does not represent exactly* classical sinograms +# in medical imaging. Integration along curves in the adjoint form of +# :func:`pylops.signalprocessing.Radon2D` is performed with respect to +# :math:`dx`, whereas canonically it is assumed to be with respect to the natural +# parametrization :math:`dl = \sqrt{(dx)^2 + (dt)^2}`. To retrieve back the +# classical sinogram we have to divide data by the jacobian +# :math:`j(x,l) = \left\vert dx/dl \right\vert = |\sin(\theta)|`. + +sinogram = np.divide(y.T, np.abs(np.sin(theta) + 1e-15)) # small shift to avoid zero-division +fig, axs = plt.subplots(1, 2, figsize=(10, 4)) +axs[0].imshow(y.T, cmap="gray") +axs[0].set_title("Data") +axs[0].axis("tight") +axs[1].imshow(sinogram, cmap="gray") +axs[1].set_title("Sinogram in medical imaging") +axs[1].axis("tight") +fig.tight_layout() ############################################################################### -# Finally we take advantage of our different solvers and try to invert the -# modelling operator both in a least-squares sense and using TV-reg. +# We will not pursue further working with the "true sinogram", but will +# reconstruct the original phantom directly from ``y``. For this we take advantage +# of our different solvers and try to invert the modelling operator both in a +# least-squares sense and using TV-reg. Dop = [ pylops.FirstDerivative( (nx, ny), axis=0, edge=True, kind="backward", dtype=np.float64 From 25315b168c8136e93b248925e8b07620fcc81b0a Mon Sep 17 00:00:00 2001 From: mrava87 Date: Tue, 3 Oct 2023 20:42:01 +0300 Subject: [PATCH 52/69] fix: change make rule for tests --- Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Makefile b/Makefile index 0ac4b3d5..2db298d1 100755 --- a/Makefile +++ b/Makefile @@ -31,7 +31,7 @@ dev-install_conda: tests: make pythoncheck - $(PYTHON) setup.py test + pytest doc: cd docs && rm -rf source/api/generated && rm -rf source/gallery &&\ From f57a2a2e963e8c892b75f5c593696aff814fa1cf Mon Sep 17 00:00:00 2001 From: mrava87 Date: Thu, 5 Oct 2023 19:42:31 +0300 Subject: [PATCH 53/69] fix: change handling of version in pyproject.toml/setup.cfg Setup parameters are moved to pyproject.toml as the current setup installs pylops0.0.0 (does not recognize the correct version). --- pyproject.toml | 47 +++++++++++++++++++++++++++++++++++++++++++- setup.cfg | 53 -------------------------------------------------- 2 files changed, 46 insertions(+), 54 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 8fd8d67e..c93b77fb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,3 +1,48 @@ [build-system] -requires = ["setuptools", "setuptools-scm"] +requires = ["setuptools>=45", "setuptools_scm[toml]>=6.2"] build-backend = "setuptools.build_meta" + +[project] +name = "pylops" +description = "Python library implementing linear operators to allow solving large-scale optimization problems" +readme = "README.md" +authors = [ + {name = "Matteo Ravasi", email = "matteoravasi@gmail.com"}, +] +license = {file = "LICENSE.md"} +keywords = ["algebra", "inverse problems", "large-scale optimization"] +classifiers = [ + "Development Status :: 5 - Production/Stable", + "Intended Audience :: Developers", + "Intended Audience :: Science/Research", + "Intended Audience :: Education", + "License :: OSI Approved :: GNU Lesser General Public License v3 (LGPLv3)", + "Natural Language :: English", + "Operating System :: OS Independent", + "Programming Language :: Python :: 3 :: Only", + "Programming Language :: Python :: 3.8", + "Programming Language :: Python :: 3.9", + "Programming Language :: Python :: 3.10", + "Topic :: Scientific/Engineering : Mathematics", +] +dependencies = [ + "numpy >= 1.21.0", + "scipy >= 1.4.0", +] +dynamic = ["version"] + +[project.optional-dependencies] +advanced = [ + "llvmlite", + "numba", + "pyfftw", + "PyWavelets", + "scikit-fmm", + "spgl1", +] + +[tool.setuptools.packages.find] +exclude = ["pytests"] + +[tool.setuptools_scm] +version_file = "pylops/version.py" diff --git a/setup.cfg b/setup.cfg index 67636c69..4c13f09c 100755 --- a/setup.cfg +++ b/setup.cfg @@ -1,56 +1,3 @@ -[metadata] -name = pylops -fullname = PyLops -description = Python library implementing linear operators to allow solving large-scale optimization problems -long_description = file: README.md -long_description_content_type = text/markdown -author = The PyLops Development Team -author_email = matteoravasi@gmail.com -maintainer = "Matteo Ravasi" -maintainer_email = matteoravasi@gmail.com -license = LGPL-3.0 License -license_file = LICENSE.md -platform = any -keywords = algebra, inverse problems, large-scale optimization -classifiers = - Development Status :: 5 - Production/Stable - Intended Audience :: Developers - Intended Audience :: Science/Research - Intended Audience :: Education - License :: OSI Approved :: GNU Lesser General Public License v3 (LGPLv3), - Natural Language :: English - Operating System :: OS Independent - Programming Language :: Python :: 3 :: Only - Programming Language :: Python :: 3.8 - Programming Language :: Python :: 3.9 - Programming Language :: Python :: 3.10 - Topic :: Scientific/Engineering : Mathematics -url = https://github.com/pylops/pylops - -project_urls = - Documentation = https://pylops.readthedocs.io/ - Release Notes = https://github.com/pylops/pylops/releases - Bug Tracker = https://github.com/pylops/pylops/issues - Source Code = https://github.com/pylops/pylops - -[options] -zip_safe = True -include_package_data = True -packages = find: -python_requires = >=3.8 -install_requires = - numpy >= 1.21.0 - scipy >= 1.4.0 - -[options.extras_require] -advanced = - llvmlite - numba - pyfftw - PyWavelets - scikit-fmm - spgl1 - [aliases] test=pytest From 130cc2925bd6d1a990833327d40c923fe6f580ae Mon Sep 17 00:00:00 2001 From: mrava87 Date: Thu, 5 Oct 2023 19:58:11 +0300 Subject: [PATCH 54/69] minor: fix setuptools version --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index c93b77fb..848c54c8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,5 +1,5 @@ [build-system] -requires = ["setuptools>=45", "setuptools_scm[toml]>=6.2"] +requires = ["setuptools>=60", "setuptools-scm>=8.0"] build-backend = "setuptools.build_meta" [project] From 385b10888678b64a2db5dff6993f6d3b39094913 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Thu, 5 Oct 2023 20:04:25 +0300 Subject: [PATCH 55/69] minor: try to fetch versions --- .github/workflows/build.yaml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/build.yaml b/.github/workflows/build.yaml index 63b4a4e3..2e7a8adc 100644 --- a/.github/workflows/build.yaml +++ b/.github/workflows/build.yaml @@ -23,6 +23,8 @@ jobs: if [ -f requirements.txt ]; then pip install -r requirements-dev.txt; fi - name: Install pylops run: | + git fetch --all --tags + python -m setuptools_scm pip install . - name: Test with pytest run: | From 30d421799bb64f58a56839abcaec30fde09dced1 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Thu, 5 Oct 2023 20:12:57 +0300 Subject: [PATCH 56/69] minor: change of versions of setuptools --- .github/workflows/build.yaml | 2 +- pyproject.toml | 6 +++++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/.github/workflows/build.yaml b/.github/workflows/build.yaml index 2e7a8adc..01cf1705 100644 --- a/.github/workflows/build.yaml +++ b/.github/workflows/build.yaml @@ -23,7 +23,7 @@ jobs: if [ -f requirements.txt ]; then pip install -r requirements-dev.txt; fi - name: Install pylops run: | - git fetch --all --tags + git pull --tags python -m setuptools_scm pip install . - name: Test with pytest diff --git a/pyproject.toml b/pyproject.toml index 848c54c8..4deee5a6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,5 +1,9 @@ [build-system] -requires = ["setuptools>=60", "setuptools-scm>=8.0"] +requires = [ + "setuptools >= 65", + "setuptools_scm[toml]", + "wheel", +] build-backend = "setuptools.build_meta" [project] From 2a471d9482bb8bb24c2b9a8dc1a20bd87f2ba123 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Thu, 5 Oct 2023 20:18:47 +0300 Subject: [PATCH 57/69] minor: add fetch of versions --- .github/workflows/build.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/build.yaml b/.github/workflows/build.yaml index 01cf1705..0408ec17 100644 --- a/.github/workflows/build.yaml +++ b/.github/workflows/build.yaml @@ -23,7 +23,7 @@ jobs: if [ -f requirements.txt ]; then pip install -r requirements-dev.txt; fi - name: Install pylops run: | - git pull --tags + git fetch --tags python -m setuptools_scm pip install . - name: Test with pytest From 0ba2276387716c0db4787ad87385e2863a157fdb Mon Sep 17 00:00:00 2001 From: mrava87 Date: Thu, 5 Oct 2023 20:35:24 +0300 Subject: [PATCH 58/69] minor: Update setuptools --- .github/workflows/build.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/build.yaml b/.github/workflows/build.yaml index 0408ec17..d75a4005 100644 --- a/.github/workflows/build.yaml +++ b/.github/workflows/build.yaml @@ -18,7 +18,7 @@ jobs: python-version: ${{ matrix.python-version }} - name: Install dependencies run: | - python -m pip install --upgrade pip + python -m pip install --upgrade pip setuptools pip install flake8 pytest if [ -f requirements.txt ]; then pip install -r requirements-dev.txt; fi - name: Install pylops From de821eb786ba6b2426700edc78c51225d63b9745 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Thu, 5 Oct 2023 20:55:06 +0300 Subject: [PATCH 59/69] minor: Workaround in GA to get right version with SCM --- .github/workflows/build.yaml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.github/workflows/build.yaml b/.github/workflows/build.yaml index d75a4005..027d5bc2 100644 --- a/.github/workflows/build.yaml +++ b/.github/workflows/build.yaml @@ -12,6 +12,10 @@ jobs: runs-on: ${{ matrix.platform }} steps: - uses: actions/checkout@v3 + - name: Get history and tags for SCM versioning to work + run: | + git fetch --prune --unshallow + git fetch --depth=1 origin +refs/tags/*:refs/tags/* - name: Set up Python ${{ matrix.python-version }} uses: actions/setup-python@v4 with: From a98d82666237f04d1993db6355e31bf8dd36ce98 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Thu, 5 Oct 2023 20:56:02 +0300 Subject: [PATCH 60/69] minor: Workaround in GA to get right version with SCM --- .github/workflows/build.yaml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/build.yaml b/.github/workflows/build.yaml index 027d5bc2..abf7ff7e 100644 --- a/.github/workflows/build.yaml +++ b/.github/workflows/build.yaml @@ -27,7 +27,6 @@ jobs: if [ -f requirements.txt ]; then pip install -r requirements-dev.txt; fi - name: Install pylops run: | - git fetch --tags python -m setuptools_scm pip install . - name: Test with pytest From a894e158c90a30cd1b5476be6a93f5588aa9bdf3 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Fri, 6 Oct 2023 21:18:37 +0300 Subject: [PATCH 61/69] minor: small text changes to ctscan tutorial --- tutorials/ctscan.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/tutorials/ctscan.py b/tutorials/ctscan.py index 1374de58..9017499f 100755 --- a/tutorials/ctscan.py +++ b/tutorials/ctscan.py @@ -32,7 +32,7 @@ # # where :math:`\theta` is the angle between the x-axis (:math:`x`) and # the perpendicular to the summation line and :math:`r` is the distance -# from the origin of the summation line. Radon transform in CT +# from the origin of the summation line. Radon transform in CT # corresponds to the integral of the input image along the straight line above. # To implement the integration in PyLops we simply need to express # :math:`t(r,\theta;x)` which is given by: @@ -49,9 +49,10 @@ def radoncurve(x, r, theta): + ny // 2 ) + ############################################################################### # Note that in the above implementation we added centering :math:`t \mapsto t - n_y/2` and -# :math:`r \mapsto r - n_y/2` so that origin of integration lines is exactly in the +# :math:`r \mapsto r - n_y/2` so that the origin of the integration lines is exactly in the # center of the image (centering for :math:`x` is not needed because we will use # ``centeredh=True`` in the constructor of ``Radon2D``). @@ -106,7 +107,9 @@ def radoncurve(x, r, theta): # classical sinogram we have to divide data by the jacobian # :math:`j(x,l) = \left\vert dx/dl \right\vert = |\sin(\theta)|`. -sinogram = np.divide(y.T, np.abs(np.sin(theta) + 1e-15)) # small shift to avoid zero-division +sinogram = np.divide( + y.T, np.abs(np.sin(theta) + 1e-15) +) # small shift to avoid zero-division fig, axs = plt.subplots(1, 2, figsize=(10, 4)) axs[0].imshow(y.T, cmap="gray") axs[0].set_title("Data") @@ -117,8 +120,8 @@ def radoncurve(x, r, theta): fig.tight_layout() ############################################################################### -# We will not pursue further working with the "true sinogram", but will -# reconstruct the original phantom directly from ``y``. For this we take advantage +# From now on, we will not pursue further working with the "true sinogram", instead +# we will reconstruct the original phantom directly from ``y``. For this we take advantage # of our different solvers and try to invert the modelling operator both in a # least-squares sense and using TV-reg. Dop = [ From 41dc437b9ef7e3b6d190c46463438199fe5aa85b Mon Sep 17 00:00:00 2001 From: mrava87 Date: Fri, 6 Oct 2023 21:20:49 +0300 Subject: [PATCH 62/69] minor: added fedor-goncharov to contributors --- README.md | 1 + docs/source/credits.rst | 3 ++- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index a39436be..3d1fc89f 100644 --- a/README.md +++ b/README.md @@ -148,3 +148,4 @@ A list of video tutorials to learn more about PyLops: * Aniket Singh Rawat, dikwickley * Rohan Babbar, rohanbabbar04 * Wei Zhang, ZhangWeiGeo +* Fedor Goncharov, fedor-goncharov diff --git a/docs/source/credits.rst b/docs/source/credits.rst index 036554b8..2a5daf61 100755 --- a/docs/source/credits.rst +++ b/docs/source/credits.rst @@ -19,4 +19,5 @@ Contributors * `Juan Daniel Romero `_, jdromerom * `Aniket Singh Rawat `_, dikwickley * `Rohan Babbar `_, rohanbabbar04 -* `Wei Zhang `_, ZhangWeiGeo \ No newline at end of file +* `Wei Zhang `_, ZhangWeiGeo +* `Fedor Goncharov `_, fedor-goncharov \ No newline at end of file From ce165368c550e54db6102bca3085ca955f80e36c Mon Sep 17 00:00:00 2001 From: cako Date: Mon, 9 Oct 2023 15:59:41 -0700 Subject: [PATCH 63/69] update isort precommit hook --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 14011ef6..a4ca2e16 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -15,7 +15,7 @@ repos: - --line-length=88 - repo: https://github.com/pycqa/isort - rev: 5.10.1 + rev: 5.12.0 hooks: - id: isort name: isort (python) From e782fd3bfe56bd9294db4f23f1d4b23c3bc1b3f1 Mon Sep 17 00:00:00 2001 From: cako Date: Mon, 9 Oct 2023 15:59:56 -0700 Subject: [PATCH 64/69] improve language --- pylops/waveeqprocessing/kirchhoff.py | 22 +++++++++++----------- pytests/test_kirchhoff.py | 2 +- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/pylops/waveeqprocessing/kirchhoff.py b/pylops/waveeqprocessing/kirchhoff.py index 3ea6b6f4..5c7fef59 100644 --- a/pylops/waveeqprocessing/kirchhoff.py +++ b/pylops/waveeqprocessing/kirchhoff.py @@ -137,18 +137,18 @@ class Kirchhoff(LinearOperator): .. math:: d(\mathbf{x_r}, \mathbf{x_s}, t) = - \widetilde{w}(t) * \int_V \frac{2 cos(\theta)} {v(\mathbf{x})} + \widetilde{w}(t) * \int_V \frac{2 \cos\theta} {v(\mathbf{x})} G(\mathbf{x_r}, \mathbf{x}, t) G(\mathbf{x}, \mathbf{x_s}, t) m(\mathbf{x}) \,\mathrm{d}\mathbf{x} where :math:`G(\mathbf{x}, \mathbf{x_s}, t)` and :math:`G(\mathbf{x_r}, \mathbf{x}, t)` are the Green's functions from source-to-subsurface-to-receiver and finally - :math:`\widetilde{w}(t)` is a filtered version of the wavelet :math:`w(t)` - as explained below (or the wavelet itself when ``wavfilter=False``). + :math:`\widetilde{w}(t)` is either a filtered version of the wavelet :math:`w(t)` + as explained below (``wavfilter=True``) or the wavelet itself when (``wavfilter=False``). Moreover, an angle scaling is included in the modelling operator, - where :math:`\theta=(\theta_s-\theta_r)/2` is half of the opening angle, - where :math:`\theta_s` and :math:`\theta_r` are the angles between the source-side - and receiver-side rays and the vertical at the image point. + where the reflection angle :math:`\theta=(\theta_s-\theta_r)/2` is half of the opening angle, + with :math:`\theta_s` and :math:`\theta_r` representing the angles between the source-side + and receiver-side rays and the vertical at the image point, respectively. In our implementation, the following high-frequency approximation of the Green's functions is adopted: @@ -168,14 +168,14 @@ class Kirchhoff(LinearOperator): On the other hand, when ``dynamic=True``, the amplitude scaling is defined as - * ``2D``: :math:`a(\mathbf{x}, \mathbf{y})=\frac{1}{\sqrt{dist(\mathbf{x}, \mathbf{y})}}` - * ``3D``: :math:`a(\mathbf{x}, \mathbf{y})=\frac{1}{dist(\mathbf{x}, \mathbf{y})}` + * ``2D``: :math:`a(\mathbf{x}, \mathbf{y})=\frac{1}{\sqrt{\text{dist}(\mathbf{x}, \mathbf{y})}}` + * ``3D``: :math:`a(\mathbf{x}, \mathbf{y})=\frac{1}{\text{dist}(\mathbf{x}, \mathbf{y})}` approximating the geometrical spreading of the wavefront. For ``mode=analytic``, - :math:`dist(\mathbf{x}, \mathbf{y})=\|\mathbf{x} - \mathbf{y}\|`, whilst for - ``mode=eikonal``, this is instead computed internally by the eikonal solver. + :math:`\text{dist}(\mathbf{x}, \mathbf{y})=\|\mathbf{x} - \mathbf{y}\|`, whilst for + ``mode=eikonal``, this is computed internally by the Eikonal solver. - The wavelet filtering is applied as follows: + The wavelet filtering is applied as follows [3]_: * ``2D``: :math:`\tilde{W}(f)=\sqrt{j\omega} \cdot W(f)` * ``3D``: :math:`\tilde{W}(f)=j\omega \cdot W(f)` diff --git a/pytests/test_kirchhoff.py b/pytests/test_kirchhoff.py index b210a95b..64139544 100755 --- a/pytests/test_kirchhoff.py +++ b/pytests/test_kirchhoff.py @@ -287,7 +287,7 @@ def test_kirchhoff3d(par): ) def test_kirchhoff2d_trav_vs_travsrcrec(par): """Compare 2D Kirchhoff operator forward and adjoint when using trav (original behavior) - or trav_src and trav_rec (new reccommended behaviour)""" + or trav_src and trav_rec (new recommended behaviour)""" # new behaviour Dop = Kirchhoff( From 81466584955fb08d59086f84f8ab569592a7ba6a Mon Sep 17 00:00:00 2001 From: mrava87 Date: Tue, 10 Oct 2023 15:59:25 +0300 Subject: [PATCH 65/69] minor: fixed epsdist and sign of 3d wavelet --- pylops/waveeqprocessing/kirchhoff.py | 16 ++++++++++------ pytests/test_kirchhoff.py | 2 +- 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/pylops/waveeqprocessing/kirchhoff.py b/pylops/waveeqprocessing/kirchhoff.py index 5c7fef59..7cd4a7da 100644 --- a/pylops/waveeqprocessing/kirchhoff.py +++ b/pylops/waveeqprocessing/kirchhoff.py @@ -133,7 +133,7 @@ class Kirchhoff(LinearOperator): ----- The Kirchhoff demigration operator synthesizes seismic data given a propagation velocity model :math:`v(\mathbf{x})` and a - reflectivity model :math:`m(\mathbf{x})`. In forward mode [1]_, [2]_: + reflectivity model :math:`m(\mathbf{x})`. In forward mode [1]_, [2]_, [3]_: .. math:: d(\mathbf{x_r}, \mathbf{x_s}, t) = @@ -175,7 +175,7 @@ class Kirchhoff(LinearOperator): :math:`\text{dist}(\mathbf{x}, \mathbf{y})=\|\mathbf{x} - \mathbf{y}\|`, whilst for ``mode=eikonal``, this is computed internally by the Eikonal solver. - The wavelet filtering is applied as follows [3]_: + The wavelet filtering is applied as follows [4]_: * ``2D``: :math:`\tilde{W}(f)=\sqrt{j\omega} \cdot W(f)` * ``3D``: :math:`\tilde{W}(f)=j\omega \cdot W(f)` @@ -209,14 +209,18 @@ class Kirchhoff(LinearOperator): projects data in the model domain creating an image of the subsurface reflectivity. - .. [1] Bleistein, N., Cohen, J.K., and Stockwell, J.W.. + .. [1] Bleistein, N., Cohen, J.K., and Stockwell, J.W. "Mathematics of Multidimensional Seismic Imaging, Migration and Inversion", 2000. .. [2] Santos, L.T., Schleicher, J., Tygel, M., and Hubral, P. "Seismic modeling by demigration", Geophysics, 65(4), pp. 1281-1289, 2000. - .. [3] Safron, L. "Multicomponent least-squares Kirchhoff depth migration", + .. [3] Yang, K., and Zhang, J. "Comparison between Born and Kirchhoff operators for + least-squares reverse time migration and the constraint of the propagation of the + background wavefield", Geophysics, 84(5), pp. R725-R739, 2019. + + .. [4] Safron, L. "Multicomponent least-squares Kirchhoff depth migration", MSc Thesis, 2018. """ @@ -302,7 +306,7 @@ def __init__( ) = Kirchhoff._traveltime_table(z, x, srcs, recs, vel, y=y, mode=mode) if self.dynamic: # need to add a scalar in the denominator of amplitude term to avoid - # division by 0, currently set to 1e-4 of max distance to avoid having + # division by 0, currently set to 1e-2 of max distance to avoid having # too large scaling around the source. This number may change in future # or left to the user to define epsdist = 1e-2 @@ -639,7 +643,7 @@ def _wavelet_reshaping( raise NotImplementedError("2.D wavelet currently not available") elif dimsrc == 3 and dimrec == 3 and dimv == 3: # 3D - Wfilt = W * (1j * 2 * np.pi * f) + Wfilt = W * (-1j * 2 * np.pi * f) wavfilt = np.fft.irfft(Wfilt, n=len(wav)) return wavfilt diff --git a/pytests/test_kirchhoff.py b/pytests/test_kirchhoff.py index 64139544..d601d91b 100755 --- a/pytests/test_kirchhoff.py +++ b/pytests/test_kirchhoff.py @@ -347,7 +347,7 @@ def test_kirchhoff2d_trav_vs_travsrcrec(par): ) def test_kirchhoff3d_trav_vs_travsrcrec(par): """Compare 3D Kirchhoff operator forward and adjoint when using trav (original behavior) - or trav_src and trav_rec (new reccomended behaviour)""" + or trav_src and trav_rec (new recommended behaviour)""" # new behaviour Dop = Kirchhoff( From 00f7034def7a0d26f3e17d02ba234d9d43cbef52 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sat, 14 Oct 2023 22:54:36 +0300 Subject: [PATCH 66/69] minor: fix doc to reflect wavelet sign --- pylops/waveeqprocessing/kirchhoff.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pylops/waveeqprocessing/kirchhoff.py b/pylops/waveeqprocessing/kirchhoff.py index 7cd4a7da..06c29251 100644 --- a/pylops/waveeqprocessing/kirchhoff.py +++ b/pylops/waveeqprocessing/kirchhoff.py @@ -178,7 +178,7 @@ class Kirchhoff(LinearOperator): The wavelet filtering is applied as follows [4]_: * ``2D``: :math:`\tilde{W}(f)=\sqrt{j\omega} \cdot W(f)` - * ``3D``: :math:`\tilde{W}(f)=j\omega \cdot W(f)` + * ``3D``: :math:`\tilde{W}(f)=-j\omega \cdot W(f)` Depending on the choice of ``mode`` the traveltime and amplitude of the Green's function will be also computed differently: From 35aeae521d7a0d26aae12ad911572b032d006dc2 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Thu, 26 Oct 2023 21:56:56 +0300 Subject: [PATCH 67/69] doc: added pylops-mpi to extensions page --- docs/source/extensions.rst | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/docs/source/extensions.rst b/docs/source/extensions.rst index 57c5c8d0..25220df5 100755 --- a/docs/source/extensions.rst +++ b/docs/source/extensions.rst @@ -15,7 +15,8 @@ for academic purposes. Spin-off projects that aim at extending the capabilities of PyLops are: -* `PyLops-GPU `_ : PyLops for GPU arrays (incorporated into PyLops). -* `PyLops-Distributed `_: PyLops for distributed systems with many computing nodes. +* `PyLops-MPI `_: PyLops for distributed systems with many computing nodes using MPI * `PyProximal `_: Proximal solvers which integrate with PyLops Linear Operators. -* `Curvelops `_: Python wrapper for the Curvelab 2D and 3D digital curvelet transforms. \ No newline at end of file +* `Curvelops `_: Python wrapper for the Curvelab 2D and 3D digital curvelet transforms. +* `PyLops-GPU `_ : PyLops for GPU arrays (unmantained! the core features are now incorporated into PyLops). +* `PyLops-Distributed `_ : PyLops for distributed systems with many computing nodes using Dask (unmantained!). From 47b672d807ebb290d7da54162be65572a8a898ff Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sun, 29 Oct 2023 11:37:20 +0300 Subject: [PATCH 68/69] minor: improved signal definition in sliding example --- examples/plot_sliding.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/examples/plot_sliding.py b/examples/plot_sliding.py index 2d4be7f8..ad553487 100644 --- a/examples/plot_sliding.py +++ b/examples/plot_sliding.py @@ -29,12 +29,14 @@ ############################################################################### # Let's start by creating a 1-dimensional array of size :math:`n_t` and create # a sliding operator to compute its transformed representation. -nwins = 4 + +# sliding window parameters nwin = 26 nover = 3 nop = 64 -dimd = nwin * nwins - 3 * nover +# length of input signal (chosen to ensure perfect match with sliding windows) +dimd = 95 t = np.arange(dimd) * 0.004 data = np.sin(2 * np.pi * 20 * t) From c8f54fc844ae13ecb811565f73acd87fbde32248 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Fri, 10 Nov 2023 22:01:34 +0300 Subject: [PATCH 69/69] minor: prepare for v2.2.0 --- CHANGELOG.md | 15 +++++++++++++++ docs/source/changelog.rst | 18 ++++++++++++++++++ pyproject.toml | 2 +- 3 files changed, 34 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d632bd4a..93ce3982 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,18 @@ +# 2.2.0 + +* Added `pylops.signalprocessing.NonStationaryConvolve3D` operator +* Added nd-array capabilities to `pylops.basicoperators.Identity` and `pylops.basicoperators.Zero` +* Added second implementation in `pylops.waveeqprocessing.BlendingContinuous` which is more + performant when dealing with small number of receivers +* Added `forceflat` property to operators with ambiguous `rmatvec` (`pylops.basicoperators.Block`, + `pylops.basicoperators.Bilinear`, `pylops.basicoperators.BlockDiag`, `pylops.basicoperators.HStack`, + `pylops.basicoperators.MatrixMult`, `pylops.basicoperators.VStack`, and `pylops.basicoperators.Zero`) +* Improved `dynamic` mode of `pylops.waveeqprocessing.Kirchhoff` operator +* Modified `pylops.signalprocessing.Convolve1D` to allow both filters that are both shorter and longer of the + input vector +* Modified all solvers to use `matvec/rmatvec` instead of `@/.H @` to improve performance + + # 2.1.0 * Added `pylops.signalprocessing.DCT`, `pylops.signalprocessing.NonStationaryConvolve1D`, `pylops.signalprocessing.NonStationaryConvolve2D`, `pylops.signalprocessing.NonStationaryFilters1D`, and diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index a56dd9aa..ad3d94b1 100644 --- a/docs/source/changelog.rst +++ b/docs/source/changelog.rst @@ -3,6 +3,24 @@ Changelog ========= +Version 2.2.0 +------------- + +*Released on: 11/11/2023* + +* Added :class:`pylops.signalprocessing.NonStationaryConvolve3D` operator +* Added nd-array capabilities to :class:`pylops.basicoperators.Identity` and :class:`pylops.basicoperators.Zero` +* Added second implementation in :class:`pylops.waveeqprocessing.BlendingContinuous` which is more + performant when dealing with small number of receivers +* Added `forceflat` property to operators with ambiguous `rmatvec` (:class:`pylops.basicoperators.Block`, + :class:`pylops.basicoperators.Bilinear`, :class:`pylops.basicoperators.BlockDiag`, :class:`pylops.basicoperators.HStack`, + :class:`pylops.basicoperators.MatrixMult`, :class:`pylops.basicoperators.VStack`, and :class:`pylops.basicoperators.Zero`) +* Improved `dynamic` mode of :class:`pylops.waveeqprocessing.Kirchhoff` operator +* Modified :class:`pylops.signalprocessing.Convolve1D` to allow both filters that are both shorter and longer of the + input vector +* Modified all solvers to use `matvec/rmatvec` instead of `@/.H @` to improve performance + + Version 2.1.0 ------------- diff --git a/pyproject.toml b/pyproject.toml index 4deee5a6..d2f6c854 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,7 +27,7 @@ classifiers = [ "Programming Language :: Python :: 3.8", "Programming Language :: Python :: 3.9", "Programming Language :: Python :: 3.10", - "Topic :: Scientific/Engineering : Mathematics", + "Topic :: Scientific/Engineering :: Mathematics", ] dependencies = [ "numpy >= 1.21.0",