Add TV inversion to postack and prestack

docs/source/adding.rst

@@ -85,7 +85,7 @@ And that's it, we have implemented our first linear operator!
 
 Testing the operator
 --------------------
-Being able to write an operator is not yet a guarantee of the fact the operator is correct, or in other words
+Being able to write an operator is not yet a guarantee of the fact that the operator is correct, or in other words
 that the adjoint code is actually the *adjoint* of the forward code. Luckily for us, a simple test can be performed
 to check the validity of forward and adjoint operators, the so called *dot-test*.
 

docs/source/faq.rst

@@ -1,4 +1,5 @@
 .. _faq:
 
 Frequenty Asked Questions
-=========================
+=========================
+None so far...

docs/source/roadmap.rst

@@ -30,8 +30,9 @@ Code cleaning
 
 * Change all ``np.flatten()`` into ``np.ravel()`` -
   `Issue #24 <https://github.com/Statoil/pylops/issues/24>`_.
-* Fix all ``if: return ... else: ...`` statements using the one-liner
-  ``return ... else ...`` - `Issue #26 <https://github.com/Statoil/pylops/issues/26>`_.
+* Fix all ``if: return ... else: ...`` statements to enforce a single return
+  with the same number of outputs
+  - `Issue #26 <https://github.com/Statoil/pylops/issues/26>`_.
 * Protected attributes and @property attributes in linear operator classes?
   - `Issue #27 <https://github.com/Statoil/pylops/issues/27>`_.
 
@@ -84,8 +85,11 @@ signalprocessing
   `pywavelets <https://pywavelets.readthedocs.io/en/latest/>`_ can be used as back-end -
   `Issue #21 <https://github.com/Statoil/pylops/issues/21>`_.
 
-* ``Fredholm1`` and ``Fredholm2`` operators applying Fredholm integrals
-  of first and second kind  - `Issue #31 <https://github.com/Statoil/pylops/issues/31>`_.
+* :strike:`Fredholm1 operator applying Fredholm integrals
+  of first kind`  - `Issue #31 <https://github.com/Statoil/pylops/issues/31>`_.
+
+* ``Fredholm2`` operators applying Fredholm integrals
+  of second kind  - `Issue #31 <https://github.com/Statoil/pylops/issues/31>`_.
 
 utils
 ~~~~~
@@ -96,8 +100,8 @@ Nothing so far
 waveeqprocessing
 ~~~~~~~~~~~~~~~~
 
-* Use ``numpy.matmul`` as a way to speed up integral computation (i.e., inner for loop)
-  in ``MDC`` operator - `Issue #32 <https://github.com/Statoil/pylops/issues/32>`_.
+* :strike:`numpy.matmul as a way to speed up integral computation (i.e., inner for loop)
+  in ``MDC`` operator` - `Issue #32 <https://github.com/Statoil/pylops/issues/32>`_.
 
 * ``NMO`` operator performing NMO modelling -
   `Issue #29 <https://github.com/Statoil/pylops/issues/29>`_.

examples/plot_tvreg.py

@@ -21,7 +21,8 @@
 the solution:
 
 .. math::
-        J = ||\mathbf{y} - \mathbf{I} \mathbf{x}||_2 + || \nabla \mathbf{x}||_1
+        J = \mu/2  ||\mathbf{y} - \mathbf{I} \mathbf{x}||_2 +
+        || \nabla \mathbf{x}||_1
 
 """
 
@@ -95,7 +96,7 @@
 ###############################################################################
 # Finally, we repeat the same exercise on a 2-dimensional image. In this case
 # we however consider the MRI imaging problem where the Split Bregman
-# solver shines. The data is this time created by first appling a 2D Fourier
+# solver shines. The data is created by first appling a 2D Fourier
 # Transform of the input model and by randomly sampling 60% of its values.
 x = np.load('../testdata/optimization/shepp_logan_phantom.npy')
 x = x/x.max()
@@ -130,15 +131,24 @@
 
 ###############################################################################
 # Let's attempt now to reconstruct the model using the Split Bregman
-# solver and an isotropic Besov regularization (aka L1 norm of the sum of
-# first derivatives over x and y).
+# with anisotropic TV regularization (aka sum of L1 norms of the
+# first derivatives over x and y):
+#
+# .. math::
+#         J = \mu/2 ||\mathbf{y} - \mathbf{R} \mathbf{F} \mathbf{x}||_2
+#         + || \nabla_x \mathbf{x}||_1 + || \nabla_y \mathbf{x}||_1
+
+
+#Dop = \
+#    [pylops.FirstDerivative(ny * nx, dims=(ny, nx), dir=0, edge=False, dtype=np.complex) + \
+#     pylops.FirstDerivative(ny * nx, dims=(ny, nx), dir=1, edge=False, dtype=np.complex)]
 Dop = \
-    [pylops.FirstDerivative(ny * nx, dims=(ny, nx), dir=0, edge=False, dtype=np.complex) + \
+    [pylops.FirstDerivative(ny * nx, dims=(ny, nx), dir=0, edge=False, dtype=np.complex),
      pylops.FirstDerivative(ny * nx, dims=(ny, nx), dir=1, edge=False, dtype=np.complex)]
 
 # TV
 mu = 1.5
-lamda = [0.1]
+lamda = [0.1, 0.1]
 niter = 20
 niterinner = 10
 
@@ -159,12 +169,12 @@
 axs[1].axis('tight')
 
 fig, axs = plt.subplots(2, 1, figsize=(10, 5))
-axs[0].plot(x[ny//2], 'k', lw=3, label='x')
-axs[0].plot(xinv[ny//2], 'r', lw=5, label='xinv TV')
+axs[0].plot(x[ny//2], 'k', lw=5, label='x')
+axs[0].plot(xinv[ny//2], 'r', lw=3, label='xinv TV')
 axs[0].set_title('Horizontal section')
 axs[0].legend()
-axs[1].plot(x[:, nx//2], 'k', lw=3, label='x')
-axs[1].plot(xinv[:, nx//2], 'r', lw=5, label='xinv TV')
+axs[1].plot(x[:, nx//2], 'k', lw=5, label='x')
+axs[1].plot(xinv[:, nx//2], 'r', lw=3, label='xinv TV')
 axs[1].set_title('Vertical section')
 axs[1].legend()
 

pylops/avo/poststack.py

@@ -8,6 +8,7 @@
 from pylops.utils import dottest as Dottest
 from pylops import MatrixMult, FirstDerivative, SecondDerivative, Laplacian
 from pylops.optimization.leastsquares import RegularizedInversion
+from pylops.optimization.sparsity import SplitBregman
 
 logging.basicConfig(format='%(levelname)s: %(message)s', level=logging.WARNING)
 
@@ -94,7 +95,7 @@ def PoststackLinearModelling(wav, nt0, spatdims=None, explicit=False):
 
 
 def PoststackInversion(data, wav, m0=None, explicit=False, simultaneous=False,
-                       epsI=None, epsR=None, dottest=False,
+                       epsI=None, epsR=None, dottest=False, epsRL1=None,
                        **kwargs_solver):
     r"""Post-stack linearized seismic inversion.
 
@@ -129,11 +130,13 @@ def PoststackInversion(data, wav, m0=None, explicit=False, simultaneous=False,
         Damping factor for additional Laplacian regularization term
     dottest : :obj:`bool`, optional
         Apply dot-test
+    epsRL1 : :obj:`float`, optional
+        Damping factor for additional blockiness regularization term
     **kwargs_solver
         Arbitrary keyword arguments for :py:func:`scipy.linalg.lstsq`
         solver (if ``explicit=True`` and  ``epsR=None``)
         or :py:func:`scipy.sparse.linalg.lsqr` solver (if ``explicit=False``
-        and/or ``epsR`` is not ``None``))
+        and/or ``epsR`` is not ``None``)
 
     Returns
     -------
@@ -161,9 +164,13 @@ def PoststackInversion(data, wav, m0=None, explicit=False, simultaneous=False,
       if ``simultaneous=True``)
     * ``explicit=False`` and ``epsR=None``: the iterative solver
       :py:func:`scipy.sparse.linalg.lsqr` is used
-    * ``explicit=False`` with ``epsR``: the iterative solver
-      :py:func:`pylops.optimization.leastsquares.RegularizedInversion` is used
-      to solve the spatially regularized problem.
+    * ``explicit=False`` with ``epsR`` and ``epsRL1=None``: the iterative
+      solver :py:func:`pylops.optimization.leastsquares.RegularizedInversion`
+      is used to solve the spatially regularized problem.
+    * ``explicit=False`` with ``epsR`` and ``epsRL1``: the iterative
+      solver :py:func:`pylops.optimization.sparsity.SplitBregman`
+      is used to solve the blockiness-promoting (in vertical direction)
+      and spatially regularized (in additional horizontal directions) problem.
 
     Note that the convergence of iterative solvers such as
     :py:func:`scipy.sparse.linalg.lsqr` can be very slow for this type of
@@ -237,18 +244,61 @@ def PoststackInversion(data, wav, m0=None, explicit=False, simultaneous=False,
             # solve unregularized normal equations simultaneously with lop
             minv = lsqr(PPop, datar, **kwargs_solver)[0]
     else:
-        # inversion with spatial regularization
-        if dims==1:
-            Regop = SecondDerivative(nt0, dtype=PPop.dtype)
-        elif dims==2:
-            Regop = Laplacian((nt0, nx), dtype=PPop.dtype)
-        else:
-            Regop = Laplacian((nt0, nx, ny), dirs=(1, 2), dtype=PPop.dtype)
+        if epsRL1 is None:
+            # L2 inversion with spatial regularization
+            if dims==1:
+                Regop = SecondDerivative(nt0, dtype=PPop.dtype)
+            elif dims==2:
+                Regop = Laplacian((nt0, nx), dtype=PPop.dtype)
+            else:
+                Regop = Laplacian((nt0, nx, ny), dirs=(1, 2), dtype=PPop.dtype)
 
-        minv = RegularizedInversion(PPop, [Regop], data.flatten(),
-                                    x0=None if m0 is None else m0.flatten(),
-                                    epsRs=[epsR], returninfo=False,
-                                    **kwargs_solver)
+            minv = RegularizedInversion(PPop, [Regop], data.flatten(),
+                                        x0=None if m0 is None else m0.flatten(),
+                                        epsRs=[epsR], returninfo=False,
+                                        **kwargs_solver)
+        else:
+            # Blockiness-promoting inversion with spatial regularization
+            if dims == 1:
+                RegL1op = FirstDerivative(nt0, dtype=PPop.dtype)
+                RegL2op = None
+            elif dims == 2:
+                RegL1op = FirstDerivative(nt0*nx, dims=(nt0, nx),
+                                          dir=0, dtype=PPop.dtype)
+                RegL2op = SecondDerivative(nt0*nx, dims=(nt0, nx),
+                                           dir=1, dtype=PPop.dtype)
+            else:
+                RegL1op = FirstDerivative(nt0*nx*ny, dims=(nt0, nx, ny),
+                                          dir=0, dtype=PPop.dtype)
+                RegL2op = Laplacian((nt0, nx, ny), dirs=(1, 2),
+                                    dtype=PPop.dtype)
+
+            if 'mu' in kwargs_solver.keys():
+                mu = kwargs_solver['mu']
+                kwargs_solver.pop('mu')
+            else:
+                mu = 1.
+            if 'niter_outer' in kwargs_solver.keys():
+                niter_outer = kwargs_solver['niter_outer']
+                kwargs_solver.pop('niter_outer')
+            else:
+                niter_outer = 3
+            if 'niter_inner' in kwargs_solver.keys():
+                niter_inner = kwargs_solver['niter_inner']
+                kwargs_solver.pop('niter_inner')
+            else:
+                niter_inner = 5
+            if not isinstance(epsRL1, (list, tuple)):
+                epsRL1 = list([epsRL1])
+            if not isinstance(epsR, (list, tuple)):
+                epsR = list([epsR])
+            minv = SplitBregman(PPop, [RegL1op], data.ravel(),
+                                RegsL2=[RegL2op], epsRL1s=epsRL1,
+                                epsRL2s=epsR, mu=mu,
+                                niter_outer=niter_outer,
+                                niter_inner=niter_inner,
+                                x0=None if m0 is None else m0.flatten(),
+                                **kwargs_solver)[0]
 
     # compute residual
     if epsR is None: