Merge branch 'dev' into patch-N-dims

Author: cako

.pre-commit-config.yaml

@@ -19,4 +19,12 @@ repos:
     hooks:
       - id: isort
         name: isort (python)
-        args: ["--profile", "black", "--filter-files", "--line-length=88"]
+        args:
+          [
+            "--profile",
+            "black",
+            "--skip",
+            "__init__.py",
+            "--filter-files",
+            "--line-length=88",
+          ]

MIGRATION_V1_V2.md

@@ -7,3 +7,4 @@ should be used as a checklist when converting a piece of code using PyLops from
 
 - Several operators have deprecated `N` as a keyword. To migrate, pass only `dims` if both `N` and `dims` are currently being passed.
 If only `N` is being passed, ensure it is being passed as a value and not a keyword argument (e.g., change `Flip(N=100)` to `Flip(100)`).
+- `utils.dottest`: The relative tolerance is new set via `rtol` (before `tol`), and absolute tolerance is new 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`.

examples/plot_nmo.py

@@ -296,7 +296,7 @@ def _rmatvec(self, y):
 # and adjoint transforms truly are adjoints of each other.
 
 NMOOp = NMO(t, x, vel_t)
-dottest(NMOOp, *NMOOp.shape, tol=1e-4)
+dottest(NMOOp, *NMOOp.shape, rtol=1e-4)
 
 ###############################################################################
 # NMO using :py:class:`pylops.Spread`
@@ -370,7 +370,7 @@ def create_tables(taxis, haxis, vels_rms):
     dtable=nmo_dtable,  # Table of weights for linear interpolation
     engine="numba",  # numba or numpy
 ).H  # To perform NMO *correction*, we need the adjoint
-dottest(SpreadNMO, *SpreadNMO.shape, tol=1e-4)
+dottest(SpreadNMO, *SpreadNMO.shape, rtol=1e-4)
 
 ###############################################################################
 # We see it passes the dot test, but are the results right? Let's find out.

examples/plot_slopeest.py

@@ -12,24 +12,34 @@
 compressed and the sparse nature of the Seislet transform can also be used to
 precondition sparsity-promoting inverse problems.
 
+We will show examples of a variety of different settings, including a comparison
+with the original implementation in [1].
+
+.. [1] van Vliet, L. J.,  Verbeek, P. W., "Estimators for orientation and
+    anisotropy in digitized images", Journal ASCI Imaging Workshop. 1995.
+
 """
 
 import matplotlib.pyplot as plt
 import numpy as np
+from matplotlib.image import imread
 from mpl_toolkits.axes_grid1 import make_axes_locatable
 
 import pylops
 from pylops.signalprocessing.Seislet import _predict_trace
+from pylops.utils.signalprocessing import slope_estimate
 
 plt.close("all")
 np.random.seed(10)
 
 ###############################################################################
+# Python logo
+# -----------
 # To start we import a 2d image and estimate the local slopes of the image.
 im = np.load("../testdata/python.npy")[..., 0]
 im = im / 255.0 - 0.5
 
-slopes, anisotropy = pylops.utils.signalprocessing.slope_estimate(im, smooth=7)
+slopes, anisotropy = slope_estimate(im, smooth=7)
 angles = -np.rad2deg(np.arctan(slopes))
 
 ###############################################################################
@@ -51,6 +61,8 @@
 fig.tight_layout()
 
 ###############################################################################
+# Seismic data
+# ------------
 # We can now repeat the same using some seismic data. We will first define
 # a single trace and a slope field, apply such slope field to the trace
 # recursively to create the other traces of the data and finally try to recover
@@ -87,7 +99,7 @@
     d[:, ix] = tr
 
 # Estimate slopes
-slope_est, _ = pylops.utils.signalprocessing.slope_estimate(d, dt, dx, smooth=10)
+slope_est, _ = slope_estimate(d, dt, dx, smooth=10)
 slope_est *= -1
 
 ###############################################################################
@@ -123,6 +135,120 @@
 fig.tight_layout()
 
 ###############################################################################
+# Concentric circles
+# ------------------
+# The original paper by van Vliet and Verbeek [1] has an example with concentric
+# circles. We recover their original images and compare our implementation with
+# theirs.
+def rgb2gray(rgb):
+    return np.dot(rgb[..., :3], [0.2989, 0.5870, 0.1140])
+
+
+circles_input = rgb2gray(imread("../testdata/slope_estimate/concentric.png"))
+circles_angles = rgb2gray(imread("../testdata/slope_estimate/concentric_angles.png"))
+
+slope, anisos_sm0 = slope_estimate(circles_input, smooth=0)
+angles_sm0 = np.rad2deg(np.arctan(slope))
+
+slope, anisos_sm4 = slope_estimate(circles_input, smooth=4)
+angles_sm4 = np.rad2deg(np.arctan(slope))
+
+###############################################################################
+fig, axs = plt.subplots(2, 3, figsize=(6, 4), sharex=True, sharey=True)
+axs[0, 0].imshow(circles_input, cmap="gray", aspect="equal")
+axs[0, 0].set(title="Original Image")
+cax = make_axes_locatable(axs[0, 0]).append_axes("right", size="5%", pad=0.05)
+cax.axis("off")
+
+axs[1, 0].imshow(-circles_angles, cmap="RdBu_r")
+axs[1, 0].set(title="Original Angles")
+cax = make_axes_locatable(axs[1, 0]).append_axes("right", size="5%", pad=0.05)
+cax.axis("off")
+
+im = axs[0, 1].imshow(angles_sm0, cmap="RdBu_r", vmin=-90, vmax=90)
+cax = make_axes_locatable(axs[0, 1]).append_axes("right", size="5%", pad=0.05)
+cb = fig.colorbar(im, cax=cax, orientation="vertical")
+axs[0, 1].set(title="Angles (smooth=0)")
+
+im = axs[1, 1].imshow(angles_sm4, cmap="RdBu_r", vmin=-90, vmax=90)
+cax = make_axes_locatable(axs[1, 1]).append_axes("right", size="5%", pad=0.05)
+cb = fig.colorbar(im, cax=cax, orientation="vertical")
+axs[1, 1].set(title="Angles (smooth=4)")
+
+im = axs[0, 2].imshow(anisos_sm0, cmap="Reds", vmin=0, vmax=1)
+cax = make_axes_locatable(axs[0, 2]).append_axes("right", size="5%", pad=0.05)
+cb = fig.colorbar(im, cax=cax, orientation="vertical")
+axs[0, 2].set(title="Anisotropy (smooth=0)")
+
+im = axs[1, 2].imshow(anisos_sm4, cmap="Reds", vmin=0, vmax=1)
+cax = make_axes_locatable(axs[1, 2]).append_axes("right", size="5%", pad=0.05)
+cb = fig.colorbar(im, cax=cax, orientation="vertical")
+axs[1, 2].set(title="Anisotropy (smooth=4)")
+
+for ax in axs.ravel():
+    ax.axis("off")
+fig.tight_layout()
+
+
+###############################################################################
+# Core samples
+# ------------------
+# The original paper by van Vliet and Verbeek [1] also has an example with images
+# of core samples. Since the original paper does not have a scale with which to
+# plot the angles, we have chosen ours it to match their image as closely as
+# possible.
+
+core_input = rgb2gray(imread("../testdata/slope_estimate/core_sample.png"))
+core_angles = rgb2gray(imread("../testdata/slope_estimate/core_sample_orientation.png"))
+core_aniso = rgb2gray(imread("../testdata/slope_estimate/core_sample_anisotropy.png"))
+
+
+slope, anisos_sm4 = slope_estimate(core_input, smooth=4)
+angles_sm4 = np.rad2deg(np.arctan(slope))
+
+slope, anisos_sm8 = slope_estimate(core_input, smooth=8)
+angles_sm8 = np.rad2deg(np.arctan(slope))
+
+###############################################################################
+fig, axs = plt.subplots(1, 6, figsize=(10, 6))
+
+axs[0].imshow(core_input, cmap="gray_r", aspect="equal")
+axs[0].set(title="Original\nImage")
+cax = make_axes_locatable(axs[0]).append_axes("right", size="20%", pad=0.05)
+cax.axis("off")
+
+axs[1].imshow(-core_angles, cmap="YlGnBu_r")
+axs[1].set(title="Original\nAngles")
+cax = make_axes_locatable(axs[1]).append_axes("right", size="20%", pad=0.05)
+cax.axis("off")
+
+im = axs[2].imshow(angles_sm8, cmap="YlGnBu_r", vmin=-49, vmax=-11)
+cax = make_axes_locatable(axs[2]).append_axes("right", size="20%", pad=0.05)
+cb = fig.colorbar(im, cax=cax, orientation="vertical")
+axs[2].set(title="Angles\n(smooth=8)")
+
+im = axs[3].imshow(angles_sm4, cmap="YlGnBu_r", vmin=-49, vmax=-11)
+cax = make_axes_locatable(axs[3]).append_axes("right", size="20%", pad=0.05)
+cb = fig.colorbar(im, cax=cax, orientation="vertical")
+axs[3].set(title="Angles\n(smooth=4)")
+
+im = axs[4].imshow(anisos_sm8, cmap="Reds", vmin=0, vmax=1)
+cax = make_axes_locatable(axs[4]).append_axes("right", size="20%", pad=0.05)
+cb = fig.colorbar(im, cax=cax, orientation="vertical")
+axs[4].set(title="Anisotropy\n(smooth=8)")
+
+im = axs[5].imshow(anisos_sm4, cmap="Reds", vmin=0, vmax=1)
+cax = make_axes_locatable(axs[5]).append_axes("right", size="20%", pad=0.05)
+cb = fig.colorbar(im, cax=cax, orientation="vertical")
+axs[5].set(title="Anisotropy\n(smooth=4)")
+
+for ax in axs.ravel():
+    ax.axis("off")
+fig.tight_layout()
+
+###############################################################################
+# Final considerations
+# --------------------
 # As you can see the Structure Tensor algorithm is a very fast, general purpose
 # algorithm that can be used to estimate local slopes to input datasets of
-# very different nature.
+# very different natures.

pylops/__init__.py

@@ -1,61 +1,70 @@
-# isort: skip_file
-from .LinearOperator import LinearOperator
-from .basicoperators import Regression
-from .basicoperators import LinearRegression
-from .basicoperators import MatrixMult
-from .basicoperators import Identity
-from .basicoperators import Zero
-from .basicoperators import Diagonal
-from .basicoperators import Flip
-from .basicoperators import Symmetrize
-from .basicoperators import Spread
-from .basicoperators import Transpose
-from .basicoperators import Roll
-from .basicoperators import Pad
-from .basicoperators import Sum
-from .basicoperators import FunctionOperator
-from .basicoperators import MemoizeOperator
+"""
+PyLops
+======
 
-from .basicoperators import VStack
-from .basicoperators import HStack
-from .basicoperators import Block
-from .basicoperators import BlockDiag
-from .basicoperators import Kronecker
-from .basicoperators import Real
-from .basicoperators import Imag
-from .basicoperators import Conj
+Linear operators and inverse problems are at the core of many of the most used
+algorithms in signal processing, image processing, and remote sensing.
+When dealing with small-scale problems, the Python numerical scientific
+libraries `numpy <http://www.numpy.org>`_
+and `scipy <http://www.scipy.org/scipylib/index.html>`_  allow to perform most
+of the underlying matrix operations (e.g., computation of matrix-vector
+products and manipulation of matrices) in a simple and expressive way.
 
-from .basicoperators import CausalIntegration
-from .basicoperators import FirstDerivative
-from .basicoperators import SecondDerivative
-from .basicoperators import Laplacian
-from .basicoperators import Gradient
-from .basicoperators import FirstDirectionalDerivative
-from .basicoperators import SecondDirectionalDerivative
-from .basicoperators import Restriction
-from .basicoperators import Smoothing1D
-from .basicoperators import Smoothing2D
+Many useful operators, however, do not lend themselves to an explicit matrix
+representation when used to solve large-scale problems. PyLops operators,
+on the other hand, still represent a matrix and can be treated in a similar
+way, but do not rely on the explicit creation of a dense (or sparse) matrix
+itself. Conversely, the forward and adjoint operators are represented by small
+pieces of codes that mimic the effect of the matrix on a vector or
+another matrix.
 
-from .avo.poststack import PoststackLinearModelling
-from .avo.prestack import PrestackWaveletModelling, PrestackLinearModelling
+Luckily, many iterative methods (e.g. cg, lsqr) do not need to know the
+individual entries of a matrix to solve a linear system. Such solvers only
+require the computation of forward and adjoint matrix-vector products as
+done for any of the PyLops operators.
 
-from .optimization.solver import cg, cgls
-from .optimization.leastsquares import NormalEquationsInversion, RegularizedInversion
-from .optimization.leastsquares import PreconditionedInversion
-from .optimization.sparsity import IRLS, OMP, ISTA, FISTA, SPGL1, SplitBregman
+PyLops provides
+  1. A general construct for creating Linear Operators
+  2. An extensive set of commonly used linear operators
+  3. A set of least-squares and sparse solvers for linear operators.
 
-from .utils.seismicevents import makeaxis, linear2d, parabolic2d
-from .utils.tapers import hanningtaper, cosinetaper, taper2d, taper3d
-from .utils.wavelets import ricker, gaussian
-from .utils.utils import Report
-from .utils.deps import *
+Available subpackages
+---------------------
+basicoperators
+    Basic Linear Operators
+signalprocessing
+    Linear Operators for Signal Processing operations
+avo
+    Linear Operators for Seismic Reservoir Characterization
+waveeqprocessing
+    Linear Operators for Wave Equation oriented processing
+optimization
+    Solvers
+utils
+    Utility routines
+
+"""
 
-from . import avo
-from . import basicoperators
-from . import optimization
-from . import signalprocessing
-from . import utils
-from . import waveeqprocessing
+from .LinearOperator import LinearOperator
+from .basicoperators import *
+
+from . import (
+    avo,
+    basicoperators,
+    optimization,
+    signalprocessing,
+    utils,
+    waveeqprocessing,
+)
+from .avo.poststack import *
+from .avo.prestack import *
+from .optimization.solver import *
+from .optimization.leastsquares import *
+from .optimization.sparsity import *
+from .utils.seismicevents import *
+from .utils.tapers import *
+from .utils.utils import Report
+from .utils.wavelets import *
 
 try:
     from .version import version as __version__

pylops/avo/__init__.py

@@ -1 +1,34 @@
-# isort: skip_file
+"""
+AVO Operators
+=============
+
+The subpackage avo provides linear operators and applications aimed at
+solving various inverse problems in the area of Seismic Reservoir
+Characterization.
+
+A list of available operators present in pylops.avo:
+
+    AVOLinearModelling                      AVO modelling.
+    PoststackLinearModelling                Post-stack seismic modelling.
+    PrestackLinearModelling                 Pre-stack seismic modelling.
+    PrestackWaveletModelling                Pre-stack modelling operator for wavelet.
+
+and a list of applications:
+
+    PoststackInversion                      Post-stack seismic inversion.
+    PrestackInversion                       Pre-stack seismic inversion.
+
+"""
+
+from .poststack import *
+from .prestack import *
+
+
+__all__ = [
+    "AVOLinearModelling",
+    "PoststackLinearModelling",
+    "PrestackWaveletModelling",
+    "PrestackLinearModelling",
+    "PoststackInversion",
+    "PrestackInversion",
+]