Merge pull request #347 from cako/patch-dimsd

Implement `dimsd` for all operators with `dims`

Author: mrava87

MIGRATION_V1_V2.md

@@ -12,3 +12,5 @@ should be used as a checklist when converting a piece of code using PyLops from
 - `utils.dottest`: Change `tol` into `rtol`. Absolute tolerance is now also supported via the keyword `atol`.
   When calling it with purely positional arguments, note that after `rtol` comes now first `atol` before `complexflag`.
   When using `raiseerror=True` it now emits an `AttributeError` instead of a `ValueError`.
+- `dims_fft` in the FFT operators is deprecated in favor of `dimsd`.
+- `dims_d` in `Sum` is deprecated in favor or `dimsd`

examples/plot_matrixmult.py

@@ -196,11 +196,11 @@
 #
 # and apply it using the same implementation of the
 # :py:class:`pylops.MatrixMult` operator by simply telling the operator how we
-# want the model to be organized through the ``dims`` input parameter.
+# want the model to be organized through the ``otherdims`` input parameter.
 A = np.array([[1.0, 2.0], [4.0, 5.0]])
 x = np.array([[1.0, 1.0], [2.0, 2.0], [3.0, 3.0]])
 
-Aop = pylops.MatrixMult(A, dims=(3,), dtype="float64")
+Aop = pylops.MatrixMult(A, otherdims=(3,), dtype="float64")
 y = Aop * x.ravel()
 
 xest, istop, itn, r1norm, r2norm = lsqr(Aop, y, damp=1e-10, iter_lim=10, show=0)[0:5]

pylops/avo/poststack.py

@@ -87,7 +87,7 @@ def _PoststackLinearModelling(
         M = ncp.dot(C, D)
         if sparse:
             M = get_csc_matrix(wav)(M)
-        Pop = _MatrixMult(M, dims=spatdims, dtype=dtype, **args_MatrixMult)
+        Pop = _MatrixMult(M, otherdims=spatdims, dtype=dtype, **args_MatrixMult)
     else:
         # Create wavelet operator
         if len(wav.shape) == 1:
@@ -102,7 +102,7 @@ def _PoststackLinearModelling(
         else:
             Cop = _MatrixMult(
                 nonstationary_convmtx(wav, nt0, hc=wav.shape[1] // 2, pad=(nt0, nt0)),
-                dims=spatdims,
+                otherdims=spatdims,
                 dtype=dtype,
                 **args_MatrixMult
             )
@@ -350,7 +350,7 @@ def PoststackInversion(
                     minv = get_lstsq(data)(PP, datarn, **kwargs_solver)[0]
                 else:
                     # solve regularized normal equations simultaneously
-                    PPop_reg = MatrixMult(PP, dims=nspatprod)
+                    PPop_reg = MatrixMult(PP, otherdims=nspatprod)
                     if ncp == np:
                         minv = lsqr(PPop_reg, datar.ravel(), **kwargs_solver)[0]
                     else:
@@ -364,7 +364,7 @@ def PoststackInversion(
                 # create regularized normal eqs. and solve them simultaneously
                 PP = ncp.dot(PPop.A.T, PPop.A) + epsI * ncp.eye(nt0, dtype=PPop.A.dtype)
                 datarn = PPop.A.T * datar.reshape(nt0, nspatprod)
-                PPop_reg = MatrixMult(PP, dims=nspatprod)
+                PPop_reg = MatrixMult(PP, otherdims=nspatprod)
                 minv = get_lstsq(data)(PPop_reg.A, datarn.ravel(), **kwargs_solver)[0]
         else:
             # solve unregularized normal equations simultaneously with lop

pylops/avo/prestack.py

@@ -182,7 +182,7 @@ def PrestackLinearModelling(
 
         # Combine operators
         M = ncp.dot(C, ncp.dot(G, D))
-        Preop = MatrixMult(M, dims=spatdims, dtype=dtype)
+        Preop = MatrixMult(M, otherdims=spatdims, dtype=dtype)
 
     else:
         # Create wavelet operator
@@ -504,7 +504,7 @@ def PrestackInversion(
         ]
         if explicit:
             PPop = MatrixMult(
-                np.vstack([Op.A for Op in PPop]), dims=nspat, dtype=PPop[0].A.dtype
+                np.vstack([Op.A for Op in PPop]), otherdims=nspat, dtype=PPop[0].A.dtype
             )
         else:
             PPop = VStack(PPop)
@@ -560,7 +560,7 @@ def PrestackInversion(
                     minv = get_lstsq(data)(PP, datarn, **kwargs_solver)[0]
                 else:
                     # solve regularized normal equations simultaneously
-                    PPop_reg = MatrixMult(PP, dims=nspatprod)
+                    PPop_reg = MatrixMult(PP, otherdims=nspatprod)
                     if ncp == np:
                         minv = lsqr(PPop_reg, datarn.ravel(), **kwargs_solver)[0]
                     else:
@@ -574,7 +574,7 @@ def PrestackInversion(
             #    # create regularized normal eqs. and solve them simultaneously
             #    PP = np.dot(PPop.A.T, PPop.A) + epsI * np.eye(nt0*nm)
             #    datarn = PPop.A.T * datar.reshape(nt0*ntheta, nspatprod)
-            #    PPop_reg = MatrixMult(PP, dims=ntheta*nspatprod)
+            #    PPop_reg = MatrixMult(PP, otherdims=ntheta*nspatprod)
             #    minv = lstsq(PPop_reg, datarn.ravel(), **kwargs_solver)[0]
         else:
             # solve unregularized normal equations simultaneously with lop

pylops/basicoperators/CausalIntegration.py

@@ -3,7 +3,7 @@
 import numpy as np
 
 from pylops import LinearOperator
-from pylops.utils._internal import _value_or_list_like_to_array
+from pylops.utils._internal import _value_or_list_like_to_tuple
 
 
 class CausalIntegration(LinearOperator):
@@ -97,16 +97,18 @@ def __init__(
         kind="full",
         removefirst=False,
     ):
-        self.dims = _value_or_list_like_to_array(dims)
+        self.dims = _value_or_list_like_to_tuple(dims)
         self.axis = axis
         self.sampling = sampling
         self.kind = kind
         if kind == "full" and halfcurrent:  # ensure backcompatibility
             self.kind = "half"
         self.removefirst = removefirst
-        self.dimsd = self.dims.copy()
+        dimsd = list(self.dims)
         if self.removefirst:
-            self.dimsd[self.axis] -= 1
+            dimsd[self.axis] -= 1
+        self.dimsd = tuple(dimsd)
+
         self.shape = (np.prod(self.dimsd), np.prod(self.dims))
         self.dtype = np.dtype(dtype)
         self.explicit = False

pylops/basicoperators/Conj.py

@@ -1,6 +1,7 @@
 import numpy as np
 
 from pylops import LinearOperator
+from pylops.utils._internal import _value_or_list_like_to_tuple
 from pylops.utils.backend import get_array_module
 
 
@@ -41,7 +42,9 @@ class Conj(LinearOperator):
     """
 
     def __init__(self, dims, dtype="complex128"):
-        self.shape = (np.prod(np.array(dims)), np.prod(np.array(dims)))
+        self.dims = self.dimsd = _value_or_list_like_to_tuple(dims)
+
+        self.shape = (np.prod(self.dimsd), np.prod(self.dims))
         self.dtype = np.dtype(dtype)
         self.explicit = False
         self.clinear = False

pylops/basicoperators/Diagonal.py

@@ -3,6 +3,7 @@
 import numpy as np
 
 from pylops import LinearOperator
+from pylops.utils._internal import _value_or_list_like_to_tuple
 from pylops.utils.backend import get_array_module, to_cupy_conditional
 
 
@@ -64,28 +65,26 @@ def __init__(self, diag, dims=None, axis=-1, dtype="float64"):
         self.diag = diag.ravel()
         self.complex = True if ncp.iscomplexobj(self.diag) else False
 
-        if dims is None:
-            self.shape = (len(self.diag), len(self.diag))
-            self.dims = None
-            self.reshape = False
-        else:
-            diagdims = [1] * len(dims)
-            diagdims[axis] = dims[axis]
-            self.diag = self.diag.reshape(diagdims)
-            self.shape = (np.prod(dims), np.prod(dims))
-            self.dims = dims
-            self.reshape = True
+        ncp = get_array_module(diag)
+        self.diag = diag.ravel()
+        self.complex = True if ncp.iscomplexobj(self.diag) else False
+        self.dims = self.dimsd = (
+            (len(self.diag),) if dims is None else _value_or_list_like_to_tuple(dims)
+        )
+
+        diagdims = np.ones_like(self.dims)
+        diagdims[axis] = self.dims[axis]
+        self.diag = self.diag.reshape(diagdims)
+
+        self.shape = (np.prod(self.dimsd), np.prod(self.dims))
         self.dtype = np.dtype(dtype)
         self.explicit = False
 
     def _matvec(self, x):
         if type(self.diag) != type(x):
             self.diag = to_cupy_conditional(x, self.diag)
-        if not self.reshape:
-            y = self.diag * x.ravel()
-        else:
-            x = x.reshape(self.dims)
-            y = self.diag * x
+        x = x.reshape(self.dims)
+        y = self.diag * x
         return y.ravel()
 
     def _rmatvec(self, x):
@@ -95,11 +94,8 @@ def _rmatvec(self, x):
             diagadj = self.diag.conj()
         else:
             diagadj = self.diag
-        if not self.reshape:
-            y = diagadj * x.ravel()
-        else:
-            x = x.reshape(self.dims)
-            y = diagadj * x
+        x = x.reshape(self.dims)
+        y = diagadj * x
         return y.ravel()
 
     def matrix(self):

pylops/basicoperators/DirectionalDerivative.py

@@ -14,7 +14,6 @@ def FirstDirectionalDerivative(
     .. note:: At least 2 dimensions are required, consider using
       :py:func:`pylops.FirstDerivative` for 1d arrays.
 
-
     Parameters
     ----------
     dims : :obj:`tuple`
@@ -67,8 +66,10 @@ def FirstDirectionalDerivative(
     else:
         Dop = Diagonal(v.ravel(), dtype=dtype)
     Sop = Sum(dims=[len(dims)] + list(dims), axis=0, dtype=dtype)
-    ddop = Sop * Dop * Gop
-    return LinearOperator(ddop)
+    ddop = LinearOperator(Sop * Dop * Gop)
+    ddop.dims = ddop.dimsd = dims
+    ddop.sampling = sampling
+    return ddop
 
 
 def SecondDirectionalDerivative(dims, v, sampling=1, edge=False, dtype="float64"):
@@ -118,5 +119,7 @@ def SecondDirectionalDerivative(dims, v, sampling=1, edge=False, dtype="float64"
     in the literature.
     """
     Dop = FirstDirectionalDerivative(dims, v, sampling=sampling, edge=edge, dtype=dtype)
-    ddop = -Dop.H * Dop
-    return LinearOperator(ddop)
+    ddop = LinearOperator(-Dop.H * Dop)
+    ddop.dims = ddop.dimsd = dims
+    ddop.sampling = sampling
+    return ddop

pylops/basicoperators/FirstDerivative.py

@@ -4,7 +4,7 @@
 from numpy.core.multiarray import normalize_axis_index
 
 from pylops import LinearOperator
-from pylops.utils._internal import _value_or_list_like_to_array
+from pylops.utils._internal import _value_or_list_like_to_tuple
 from pylops.utils.backend import get_array_module
 
 
@@ -73,13 +73,13 @@ def __init__(
         dtype="float64",
         kind="centered",
     ):
-        self.dims = _value_or_list_like_to_array(dims)
+        self.dims = self.dimsd = _value_or_list_like_to_tuple(dims)
         self.axis = normalize_axis_index(axis, len(self.dims))
         self.sampling = sampling
         self.edge = edge
         self.kind = kind
-        N = np.prod(self.dims)
-        self.shape = (N, N)
+
+        self.shape = (np.prod(self.dimsd), np.prod(self.dims))
         self.dtype = np.dtype(dtype)
         self.explicit = False
 

pylops/basicoperators/Flip.py

@@ -3,7 +3,7 @@
 import numpy as np
 
 from pylops import LinearOperator
-from pylops.utils._internal import _value_or_list_like_to_array
+from pylops.utils._internal import _value_or_list_like_to_tuple
 
 
 class Flip(LinearOperator):
@@ -46,10 +46,10 @@ class Flip(LinearOperator):
     """
 
     def __init__(self, dims, axis=-1, dtype="float64"):
-        self.dims = _value_or_list_like_to_array(dims)
+        self.dims = self.dimsd = _value_or_list_like_to_tuple(dims)
         self.axis = axis
-        N = np.prod(self.dims)
-        self.shape = (N, N)
+
+        self.shape = (np.prod(self.dimsd), np.prod(self.dims))
         self.dtype = np.dtype(dtype)
         self.explicit = False