Overview

PyLops is an open-source Python library focused on providing a backend-agnostic, idiomatic, matrix-free library of linear operators and related computations. It is inspired by the iconic MATLAB Spot – A Linear-Operator Toolbox project.

Linear operators and inverse problems are at the core of many of the most used algorithms in signal processing, image processing, and remote sensing. For small-scale problems, matrices can be explicitly computed and manipulated with Python numerical scientific libraries such as NumPy and SciPy.

Large-scale problems often feature matrices that are prohibitive in size—but whose operations can be described by simple functions. PyLops operators exploit this to represent a linear operator not as array of numbers, but by functions which describe matrix-vector products in forward and adjoint modes. Moreover, many iterative methods (e.g. cg, lsqr) are designed to not rely on the elements of the matrix, only these matrix-vector products. PyLops offers such solvers for many different types of problems, in particular least-squares and sparsity-promoting inversions.

Get started by installing PyLops and following our quick tour.

Terminology

A common terminology is used within the entire documentation of PyLops. Every linear operator and its application to a model will be referred to as forward model (or operation)

\[\mathbf{y} = \mathbf{A} \mathbf{x}\]

while its application to a data is referred to as adjoint model (or operation)

\[\mathbf{x} = \mathbf{A}^H \mathbf{y}\]

where \(\mathbf{x}\) is called model and \(\mathbf{y}\) is called data. The operator \(\mathbf{A}:\mathbb{F}^m \to \mathbb{F}^n\) effectively maps a vector of size \(m\) in the model space to a vector of size \(n\) in the data space, conversely the adjoint operator \(\mathbf{A}^H:\mathbb{F}^n \to \mathbb{F}^m\) maps a vector of size \(n\) in the data space to a vector of size \(m\) in the model space. As linear operators mimics the effect a matrix on a vector we can also loosely refer to \(m\) as the number of columns and \(n\) as the number of rows of the operator.

Ultimately, solving an inverse problems accounts to removing the effect of \(\mathbf{A}\) from the data \(\mathbf{y}\) to retrieve the model \(\mathbf{x}\).

For a more detailed description of the concepts of linear operators, adjoints and inverse problems in general, you can head over to one of Jon Claerbout’s books such as Basic Earth Imaging.

Implementation

PyLops is build on top of the scipy class scipy.sparse.linalg.LinearOperator.

This class allows in fact for the creation of objects (or interfaces) for matrix-vector and matrix-matrix products that can ultimately be used to solve any inverse problem of the form \(\mathbf{y}=\mathbf{A}\mathbf{x}\).

As explained in the scipy LinearOperator official documentation, to construct a scipy.sparse.linalg.LinearOperator, a user is required to pass appropriate callables to the constructor of this class, or subclass it. More specifically one of the methods _matvec and _matmat must be implemented for the forward operator and one of the methods _rmatvec or _adjoint may be implemented to apply the Hermitian adjoint. The attributes/properties shape (pair of integers) and dtype (may be None) must also be provided during __init__ of this class.

Any linear operator developed within the PyLops library follows this philosophy. As explained more in details in Implementing new operators section, a linear operator is created by subclassing the scipy.sparse.linalg.LinearOperator class and _matvec and _rmatvec are implemented.

History

PyLops was initially written by Equinor It is a flexible and scalable python library for large-scale optimization with linear operators that can be tailored to our needs, and as contribution to the free software community. Since June 2021, PyLops is a NUMFOCUS Affiliated Project.

Installation

Dependencies

The PyLops project strives to create a library that is easy to install in any environment and has a very limited number of dependencies. Required dependencies are limited to:

We highly encourage using the Anaconda Python distribution or its standalone package manager Conda. Especially for Intel processors, this ensures a higher performance with no configuration. If you are interested in getting the best code performance, read carefully Advanced installation. For learning, however, the standard installation is often good enough.

Some operators have additional, optional “engines” to improve their performance. These often rely on third-party libraries which are added to the list of our optional dependencies. Optional dependencies therefore refer to those dependencies that are not strictly needed nor installed directly as part of a standard installation. For details more details, see Optional dependencies.

Step-by-step installation for users

Pip

If you are using pip, and simply type the following command in your terminal to install the PyPI distribution:

>> pip install pylops

Note that when installing via pip, only required dependencies are installed.

Docker

If you want to try PyLops but do not have Python in your local machine, you can use our Docker image instead.

After installing Docker in your computer, type the following command in your terminal (note that this will take some time the first time you type it as you will download and install the Docker image):

>> docker run -it -v /path/to/local/folder:/home/jupyter/notebook -p 8888:8888 mrava87/pylops:notebook

This will give you an address that you can put in your browser and will open a Jupyter notebook environment with PyLops and other basic Python libraries installed. Here, /path/to/local/folder is the absolute path of a local folder on your computer where you will create a notebook (or containing notebooks that you want to continue working on). Note that anything you do to the notebook(s) will be saved in your local folder. A larger image with conda a distribution is also available:

>> docker run -it -v /path/to/local/folder:/home/jupyter/notebook -p 8888:8888 mrava87/pylops:conda_notebook

Step-by-step installation for developers

Fork PyLops

Fork the PyLops repository and clone it by executing the following in your terminal:

>> git clone https://github.com/YOUR-USERNAME/pylops.git

We recommend installing dependencies into a separate environment. For that end, we provide a Makefile with useful commands for setting up the environment.

Install dependencies
Conda (recommended)

For a conda environment, run

>> make dev-install_conda

This will create and activate an environment called pylops, with all required and optional dependencies.

Pip

If you prefer a pip installation, we provide the following command

>> make dev-install

Note that, differently from the conda command, the above will not create a virtual environment. Make sure you create and activate your environment previously.

Run tests

To ensure that everything has been setup correctly, run tests:

>> make tests

Make sure no tests fail, this guarantees that the installation has been successful.

Add remote (optional)

To keep up-to-date on the latest changes while you are developing, you may optionally add the PyLops repository as a remote. Run the following command to add the PyLops repo as a remote named upstream:

>> git remote add upstream https://github.com/PyLops/pylops

From then on, you can pull changes (for example, in the dev branch) with:

>> git pull upstream dev
Install pre-commit hooks

To ensure consistency in the coding style of our developers we rely on pre-commit to perform a series of checks when you are ready to commit and push some changes. This is accomplished by means of git hooks that have been configured in the .pre-commit-config.yaml file.

In order to setup such hooks in your local repository, run:

>> pre-commit install

Once this is set up, when committing changes, pre-commit will reject and “fix” your code by running the proper hooks. At this point, the user must check the changes and then stage them before trying to commit again.

Final steps

PyLops does not enforce the use of a linter as a pre-commit hook, but we do highly encourage using one before submitting a Pull Request. A properly configured linter (flake8) can be run with:

>> make lint

In addition, it is highly encouraged to build the docs prior to submitting a Pull Request. Apart from ensuring that docstrings are properly formatted, they can aid in catching bugs during development. Build (or update) the docs with:

>> make doc

or

>> make docupdate

Advanced installation

In this section we discuss some important details regarding code performance when using PyLops.

To get the most out of PyLops operators in terms of speed you will need to follow these guidelines as much as possible or ensure that the Python libraries used by PyLops are efficiently installed in your system.

BLAS

PyLops relies on the NumPy and SciPy, and being able to link these to the most performant BLAS library will ensure optimal performance of PyLops when using only required dependencies.

We strongly encourage using the Anaconda Python distribution as NumPy and SciPy will, when available, be automatically linked to Intel MKL, the most performant library for basic linear algebra operations to date (see Markus Beuckelmann’s benchmarks). The PyPI version installed with pip, however, will default to OpenBLAS. For more information, see NumPy’s section on BLAS.

To check which BLAS NumPy and SciPy were compiled against, run the following commands in a Python interpreter:

import numpy as np
import scipy as sp
print(np.__config__.show())
print(sp.__config__.show())

Intel also provides NumPy and SciPy replacement packages in PyPI intel-numpy and intel-scipy, respectively, which link to Intel MKL. These are an option for an environment without conda that needs Intel MKL without requiring manual compilation.

Warning

intel-numpy and intel-scipy not only link against Intel MKL, but also substitute NumPy and SciPy FFTs for Intel MKL FFT. MKL FFT is not supported and may break PyLops.

Multithreading

It is important to ensure that your environment variable which sets threads is correctly assigned to the maximum number of cores you would like to use in your code. Multiprocessing parallelism in NumPy and SciPy can be controlled in different ways depending on where it comes from.

Environment variable

Library

OMP_NUM_THREADS

OpenMP

NUMEXPR_NUM_THREADS

NumExpr

OPENBLAS_NUM_THREADS

OpenBLAS

MKL_NUM_THREADS

Intel MKL

VECLIB_MAXIMUM_THREADS

Apple Accelerate (vecLib)

For example, try setting one processor to be used with (if using OpenBlas)

>> export OMP_NUM_THREADS=1
>> export NUMEXPR_NUM_THREADS=1
>> export OPENBLAS_NUM_THREADS=1

and run the following code in Python:

import os
import numpy as np
from timeit import timeit

size = 1024
A = np.random.random((size, size)),
B = np.random.random((size, size))
print("Time with %s threads: %f s" \
      %(os.environ.get("OMP_NUM_THREADS"),
        timeit(lambda: np.dot(A, B), number=4)))

Subsequently set the environment variables to 2 or any higher number of threads available in your hardware (multi-threaded), and run the same code. By looking at both the load on your processors (e.g., using top), and at the Python print statement you should see a speed-up in the second case.

Alternatively, you could set the OMP_NUM_THREADS variable directly inside your script using os.environ["OMP_NUM_THREADS"]="2", but ensure that this is done before loading NumPy.

Note

Always remember to set OMP_NUM_THREADS and other relevant variables in your environment when using PyLops

Optional dependencies

To avoid increasing the number of required dependencies, which may lead to conflicts with other libraries that you have in your system, we have decided to build some of the additional features of PyLops in such a way that if an optional dependency is not present in your Python environment, a safe fallback to one of the required dependencies will be enforced.

When available in your system, we recommend using the Conda package manager and install all the required and optional dependencies of PyLops at once using the command:

>> conda install --channel conda-forge pylops

in this case all dependencies will be installed from their Conda distributions.

Alternatively, from version 1.4.0 optional dependencies can also be installed as part of the pip installation via:

>> pip install pylops[advanced]

Dependencies are however installed from their PyPI wheels. An exception is however represented by CuPy. This library is not installed automatically. Users interested to accelerate their computations with the aid of GPUs should install it prior to installing PyLops as described in Optional Dependencies for GPU.

Note

If you are a developer, all the optional dependencies below (except GPU) can be installed automatically by cloning the repository and installing PyLops via make dev-install_conda (conda) or make dev-install (pip).

In alphabetic order:

Devito

Devito is library used to solve PDEs via the finite-difference method. It is used in PyLops to compute wavefields pylops.waveeqprocessing.AcousticWave2D

Install it via pip with

>> pip install devito
FFTW

Three different “engines” are provided by the pylops.signalprocessing.FFT operator: engine="numpy" (default), engine="scipy" and engine="fftw".

The first two engines are part of the required PyLops dependencies. The latter implements the well-known FFTW via the Python wrapper pyfftw.FFTW. While this optimized FFT tends to outperform the other two in many cases, it is not included by default. To use this library, install it manually either via conda:

>> conda install --channel conda-forge pyfftw

or via pip:

>> pip install pyfftw

Warning

Intel MKL FFT is not supported.

Numba

Although we always strive to write code for forward and adjoint operators that takes advantage of the perks of NumPy and SciPy (e.g., broadcasting, ufunc), in some case we may end up using for loops that may lead to poor performance. In those cases we may decide to implement alternative (optional) back-ends in Numba, a Just-In-Time compiler that translates a subset of Python and NumPy code into fast machine code.

A user can simply switch from the native, always available implementation to the Numba implementation by simply providing the following additional input parameter to the operator engine="numba". This is for example the case in the pylops.signalprocessing.Radon2D.

If interested to use Numba backend from conda, you will need to manually install it:

>> conda install numba

It is also advised to install the additional package icc_rt to use optimised transcendental functions as compiler intrinsics.

>> conda install --channel numba icc_rt

Through pip the equivalent would be:

>> pip install numba
>> pip install icc_rt

However, it is important to note that icc_rt will only be identified by Numba if LD_LIBRARY_PATH is properly set. If you are using a virtual environment, you can ensure this with:

>> export LD_LIBRARY_PATH=/path/to/venv/lib/:$LD_LIBRARY_PATH

To ensure that icc_rt is being recognized, run

>> numba -s | grep SVML
__SVML Information__
SVML State, config.USING_SVML                 : True
SVML Library Loaded                           : True
llvmlite Using SVML Patched LLVM              : True
SVML Operational                              : True

Numba also offers threading parallelism through a variety of Threading Layers. You may need to set the environment variable NUMBA_NUM_THREADS define how many threads to use out of the available ones (numba -s | grep "CPU Count"). It can also be checked dynamically with numba.config.NUMBA_DEFAULT_NUM_THREADS.

PyWavelets

PyWavelets is used to implement the wavelet operators. Install it via conda with:

>> conda install pywavelets

or via pip with

>> pip install PyWavelets
scikit-fmm

scikit-fmm is a library which implements the fast marching method. It is used in PyLops to compute traveltime tables in the initialization of pylops.waveeqprocessing.Kirchhoff when choosing mode="eikonal". As this may not be of interest for many users, this library has not been added to the mandatory requirements of PyLops. With conda, install it via

>> conda install --channel conda-forge scikit-fmm

or with pip via

>> pip install scikit-fmm
SPGL1

SPGL1 is used to solve sparsity-promoting basis pursuit, basis pursuit denoise, and Lasso problems in pylops.optimization.sparsity.SPGL1 solver.

Install it via pip with:

>> pip install spgl1
Sympy

This library is used to implement the describe method, which transforms PyLops operators into their mathematical expression.

Install it via conda with:

>> conda install sympy

or via pip with

>> pip install sympy
Torch

Torch used to allow seamless integration between PyLops and PyTorch operators.

Install it via conda with:

>> conda install -c pytorch pytorch

or via pip with

>> pip install torch
Optional Dependencies for GPU

PyLops will automatically check if the libraries below are installed and, in that case, use them any time the input vector passed to an operator is of compatible type. Users can, however, disable this option. For more details of GPU-accelerated PyLops read GPU Support.

CuPy

CuPy is a library used as a drop-in replacement to NumPy for GPU-accelerated computations. Since many different versions of CuPy exist (based on the CUDA drivers of the GPU), users must install CuPy prior to installing PyLops. To do so, follow their installation instructions.

cuSignal

cuSignal is a library is used as a drop-in replacement to SciPy Signal for GPU-accelerated computations. Similar to CuPy, users must install cuSignal prior to installing PyLops. To do so, follow their installation instructions.

GPU Support

Overview

From v1.12.0, PyLops supports computations on GPUs powered by CuPy (cupy-cudaXX>=8.1.0) and cuSignal (cusignal>=0.16.0). They must be installed before PyLops is installed.

Note

Set environment variables CUPY_PYLOPS=0 and/or CUSIGNAL_PYLOPS=0 to force PyLops to ignore cupy and cusignal backends. This can be also used if a previous version of cupy or cusignal is installed in your system, otherwise you will get an error when importing PyLops.

Apart from a few exceptions, all operators and solvers in PyLops can seamlessly work with numpy arrays on CPU as well as with cupy arrays on GPU. Users do simply need to consistently create operators and provide data vectors to the solvers, e.g., when using pylops.MatrixMult the input matrix must be a cupy array if the data provided to a solver is also cupy array.

Warning

Some pylops.LinearOperator methods are currently on GPU:

  • pylops.LinearOperator.eigs

  • pylops.LinearOperator.cond

  • pylops.LinearOperator.tosparse

  • pylops.LinearOperator.estimate_spectral_norm

Warning

Some solvers are currently not available on GPU:

  • pylops.optimization.sparsity.SPGL1

Example

Finally, let’s briefly look at an example. First we write a code snippet using numpy arrays which PyLops will run on your CPU:

ny, nx = 400, 400
G = np.random.normal(0, 1, (ny, nx)).astype(np.float32)
x = np.ones(nx, dtype=np.float32)

Gop = MatrixMult(G, dtype='float32')
y = Gop * x
xest = Gop / y

Now we write a code snippet using cupy arrays which PyLops will run on your GPU:

ny, nx = 400, 400
G = cp.random.normal(0, 1, (ny, nx)).astype(np.float32)
x = cp.ones(nx, dtype=np.float32)

Gop = MatrixMult(G, dtype='float32')
y = Gop * x
xest = Gop / y

The code is almost unchanged apart from the fact that we now use cupy arrays, PyLops will figure this out!

Note

The CuPy backend is in active development, with many examples not yet in the docs. You can find many other examples from the PyLops Notebooks repository.

Extensions

PyLops brings to users the power of linear operators in a simple and easy to use programming interface.

While very powerful on its own, this library is further extended to take advantage of more advanced computational resources, either in terms of multiple-node clusters or GPUs. Moreover, some independent libraries are created to use third party software that cannot be included as dependencies to our main library for licensing issues but may be useful 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.

  • PyProximal: Proximal solvers which integrate with PyLops Linear Operators.

  • Curvelops: Python wrapper for the Curvelab 2D and 3D digital curvelet transforms.

Tutorials

01. The LinearOpeator

01. The LinearOpeator

01. The LinearOpeator
02. The Dot-Test

02. The Dot-Test

02. The Dot-Test
03. Solvers

03. Solvers

03. Solvers
03. Solvers (Advanced)

03. Solvers (Advanced)

03. Solvers (Advanced)
04. Bayesian Inversion

04. Bayesian Inversion

04. Bayesian Inversion
05. Image deblurring

05. Image deblurring

05. Image deblurring
06. 2D Interpolation

06. 2D Interpolation

06. 2D Interpolation
07. Post-stack inversion

07. Post-stack inversion

07. Post-stack inversion
08. Pre-stack (AVO) inversion

08. Pre-stack (AVO) inversion

08. Pre-stack (AVO) inversion
09. Multi-Dimensional Deconvolution

09. Multi-Dimensional Deconvolution

09. Multi-Dimensional Deconvolution
10. Marchenko redatuming by inversion

10. Marchenko redatuming by inversion

10. Marchenko redatuming by inversion
11. Radon filtering

11. Radon filtering

11. Radon filtering
12. Seismic regularization

12. Seismic regularization

12. Seismic regularization
13. Deghosting

13. Deghosting

13. Deghosting
14. Seismic wavefield decomposition

14. Seismic wavefield decomposition

14. Seismic wavefield decomposition
15. Least-squares migration

15. Least-squares migration

15. Least-squares migration
16. CT Scan Imaging

16. CT Scan Imaging

16. CT Scan Imaging
17. Real/Complex Inversion

17. Real/Complex Inversion

17. Real/Complex Inversion
18. Deblending

18. Deblending

18. Deblending
19. Automatic Differentiation

19. Automatic Differentiation

19. Automatic Differentiation

01. The LinearOpeator

This first tutorial is aimed at easing the use of the PyLops library for both new users and developers.

Since PyLops heavily relies on the use of the scipy.sparse.linalg.LinearOperator class of SciPy, 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 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, chaining, etc. multiple operators in very easy and expressive way.

Let’s start by defining a simple operator that applies element-wise multiplication of the model with a vector d in forward mode and element-wise multiplication of the data with the same vector d in adjoint mode. This operator is present in PyLops under the name of pylops.Diagonal and its implementation is discussed in more details in the Implementing new operators page.

import timeit

import matplotlib.pyplot as plt
import numpy as np

import pylops

n = 10
d = np.arange(n) + 1.0
x = np.ones(n)
Dop = pylops.Diagonal(d)

First of all we apply the operator in the forward mode. This can be done in four different ways:

  • _matvec: directly applies the method implemented for forward mode

  • matvec: performs some checks before and after applying _matvec

  • *: operator used to map the special method __matmul__ which checks whether the input x is a vector or matrix and applies _matvec or _matmul accordingly.

  • @: operator used to map the special method __mul__ which performs like the * opetator

We will time these 4 different executions and see how using _matvec (or matvec) will result in the faster computation. It is thus advised to use * (or @) in examples when expressivity has priority but prefer _matvec (or matvec) for efficient implementations.

# setup command
cmd_setup = """\
import numpy as np
import pylops
n = 10
d = np.arange(n) + 1.
x = np.ones(n)
Dop = pylops.Diagonal(d)
DopH = Dop.H
"""

# _matvec
cmd1 = "Dop._matvec(x)"

# matvec
cmd2 = "Dop.matvec(x)"

# @
cmd3 = "Dop@x"

# *
cmd4 = "Dop*x"

# timing
t1 = 1.0e3 * np.array(timeit.repeat(cmd1, setup=cmd_setup, number=500, repeat=5))
t2 = 1.0e3 * np.array(timeit.repeat(cmd2, setup=cmd_setup, number=500, repeat=5))
t3 = 1.0e3 * np.array(timeit.repeat(cmd3, setup=cmd_setup, number=500, repeat=5))
t4 = 1.0e3 * np.array(timeit.repeat(cmd4, setup=cmd_setup, number=500, repeat=5))

plt.figure(figsize=(7, 2))
plt.plot(t1, "k", label=" _matvec")
plt.plot(t2, "r", label="matvec")
plt.plot(t3, "g", label="@")
plt.plot(t4, "b", label="*")
plt.axis("tight")
plt.legend()
plt.tight_layout()
linearoperator

Similarly we now consider the adjoint mode. This can be done in three different ways:

  • _rmatvec: directly applies the method implemented for adjoint mode

  • rmatvec: performs some checks before and after applying _rmatvec

  • .H*: first applies the adjoint .H which creates a new scipy.sparse.linalg._CustomLinearOperator` where _matvec and _rmatvec are swapped and then applies the new _matvec.

Once again, after timing these 3 different executions we can see see how using _rmatvec (or rmatvec) will result in the faster computation while .H* is very unefficient and slow. Note that if the adjoint has to be applied multiple times it is at least advised to create the adjoint operator by applying .H only once upfront. Not surprisingly, the linear solvers in scipy as well as in PyLops actually use matvec and rmatvec when dealing with linear operators.

# _rmatvec
cmd1 = "Dop._rmatvec(x)"

# rmatvec
cmd2 = "Dop.rmatvec(x)"

# .H* (pre-computed H)
cmd3 = "DopH*x"

# .H*
cmd4 = "Dop.H*x"

# timing
t1 = 1.0e3 * np.array(timeit.repeat(cmd1, setup=cmd_setup, number=500, repeat=5))
t2 = 1.0e3 * np.array(timeit.repeat(cmd2, setup=cmd_setup, number=500, repeat=5))
t3 = 1.0e3 * np.array(timeit.repeat(cmd3, setup=cmd_setup, number=500, repeat=5))
t4 = 1.0e3 * np.array(timeit.repeat(cmd4, setup=cmd_setup, number=500, repeat=5))

plt.figure(figsize=(7, 2))
plt.plot(t1, "k", label=" _rmatvec")
plt.plot(t2, "r", label="rmatvec")
plt.plot(t3, "g", label=".H* (pre-computed H)")
plt.plot(t4, "b", label=".H*")
plt.axis("tight")
plt.legend()
plt.tight_layout()
linearoperator

Just to reiterate once again, it is advised to call matvec and rmatvec unless PyLops linear operators are used for teaching purposes.

We now go through some other methods and special methods that are implemented in scipy.sparse.linalg.LinearOperator (and pylops.LinearOperator):

  • Op1+Op2: maps the special method __add__ and performs summation between two operators and returns a pylops.LinearOperator

  • -Op: maps the special method __neg__ and performs negation of an operators and returns a pylops.LinearOperator

  • Op1-Op2: maps the special method __sub__ and performs summation between two operators and returns a pylops.LinearOperator

  • Op1**N: maps the special method __pow__ and performs exponentiation of an operator and returns a pylops.LinearOperator

  • Op/y (and Op.div(y)): maps the special method __truediv__ and performs inversion of an operator

  • Op.eigs(): estimates the eigenvalues of the operator

  • Op.cond(): estimates the condition number of the operator

  • Op.conj(): create complex conjugate operator

Dop = pylops.Diagonal(d)

# +
print(Dop + Dop)

# -
print(-Dop)
print(Dop - 0.5 * Dop)

# **
print(Dop**3)

# * and /
y = Dop * x
print(Dop / y)

# eigs
print(Dop.eigs(neigs=3))

# cond
print(Dop.cond())

# conj
print(Dop.conj())
<10x10 LinearOperator with dtype=float64>
<10x10 LinearOperator with dtype=float64>
<10x10 LinearOperator with dtype=float64>
<10x10 LinearOperator with dtype=float64>
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[10.+0.j  9.+0.j  8.+0.j]
(10.00000000000003+0j)
<10x10 _ConjLinearOperator with dtype=float64>

To understand the effect of conj we need to look into a problem with an operator in the complex domain. Let’s create again our pylops.Diagonal operator but this time we populate it with complex numbers. We will see that the action of the operator and its complex conjugate is different even if the model is real.

n = 5
d = 1j * (np.arange(n) + 1.0)
x = np.ones(n)
Dop = pylops.Diagonal(d)

print(f"y = Dx = {Dop * x}")
print(f"y = conj(D)x = {Dop.conj() * x}")
y = Dx = [0.+1.j 0.+2.j 0.+3.j 0.+4.j 0.+5.j]
y = conj(D)x = [0.-1.j 0.-2.j 0.-3.j 0.-4.j 0.-5.j]

At this point, the concept of linear operator may sound abstract. The convinience method 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))
plt.imshow(np.abs(D))
plt.title("Dense representation of Diagonal operator")
plt.axis("tight")
plt.colorbar()
plt.tight_layout()
Dense representation of Diagonal operator

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 pylops.LinearOperator operator. This means that we can still apply any of the methods implemented in the original scipy class definition like *, as well as those in our class definition like /

Dop1 = Dop - Dop.conj()

y = Dop1 * x
print(f"x = (Dop - conj(Dop))/y = {Dop1 / y}")

D1 = Dop1.todense()

plt.figure(figsize=(5, 5))
plt.imshow(np.abs(D1))
plt.title(r"Dense representation of $|D - D^*|$")
plt.axis("tight")
plt.colorbar()
plt.tight_layout()
Dense representation of $|D - D^*|$
x = (Dop - conj(Dop))/y = [1.+0.j 1.+0.j 1.+0.j 1.+0.j 1.+0.j]

Finally, 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 inside the solver.

Dop = pylops.Diagonal(d)

y = Dop.matvec(x)
y = Dop.matvec(x)
y = Dop.rmatvec(y)

print(f"Forward evaluations: {Dop.matvec_count}")
print(f"Adjoint evaluations: {Dop.rmatvec_count}")

# Reset
Dop.reset_count()
print(f"Forward evaluations: {Dop.matvec_count}")
print(f"Adjoint evaluations: {Dop.rmatvec_count}")
Forward evaluations: 2
Adjoint evaluations: 1
Forward evaluations: 0
Adjoint evaluations: 0

This first tutorial is completed. You have seen the basic operations that can be performed using scipy.sparse.linalg.LinearOperator and our overload of such a class pylops.LinearOperator and you should be able to get started combining various PyLops operators and solving your own inverse problems.

Total running time of the script: ( 0 minutes 1.049 seconds)

Gallery generated by Sphinx-Gallery

02. The Dot-Test

One of the most important aspect of writing a Linear operator is to be able to verify that the code implemented in forward mode and the code implemented in adjoint mode are effectively adjoint to each other. If this is the case, your Linear operator will successfully pass the so-called dot-test. Refer to the Notes section of pylops.utils.dottest) for a more detailed description.

In this example, I will show you how to use the dot-test for a variety of operator when model and data are either real or complex numbers.

import matplotlib.gridspec as pltgs
import matplotlib.pyplot as plt

# pylint: disable=C0103
import numpy as np

import pylops
from pylops.utils import dottest

plt.close("all")

Let’s start with something very simple. We will make a pylops.MatrixMult operator and verify that its implementation passes the dot-test. For this time, we will do this step-by-step, replicating what happens in the pylops.utils.dottest routine.

N, M = 5, 3
Mat = np.arange(N * M).reshape(N, M)
Op = pylops.MatrixMult(Mat)

v = np.random.randn(N)
u = np.random.randn(M)

# Op * u
y = Op.matvec(u)
# Op'* v
x = Op.rmatvec(v)

yy = np.dot(y, v)  # (Op  * u)' * v
xx = np.dot(u, x)  # u' * (Op' * v)

print(f"Dot-test {np.abs((yy - xx) / ((yy + xx + 1e-15) / 2)):.2e}")
Dot-test 1.37e-16

And here is a visual intepretation of what a dot-test is

gs = pltgs.GridSpec(1, 9)
fig = plt.figure(figsize=(7, 3))
ax = plt.subplot(gs[0, 0:2])
ax.imshow(Op.A, cmap="rainbow")
ax.set_title(r"$(Op*$", size=20, fontweight="bold")
ax.set_xticks(np.arange(M - 1) + 0.5)
ax.set_yticks(np.arange(N - 1) + 0.5)
ax.grid(linewidth=3, color="white")
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])
ax.axis("tight")
ax = plt.subplot(gs[0, 2])
ax.imshow(u[:, np.newaxis], cmap="rainbow")
ax.set_title(r"$u)^T$", size=20, fontweight="bold")
ax.set_xticks([])
ax.set_yticks(np.arange(M - 1) + 0.5)
ax.grid(linewidth=3, color="white")
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])
ax.axis("tight")
ax = plt.subplot(gs[0, 3])
ax.imshow(v[:, np.newaxis], cmap="rainbow")
ax.set_title(r"$v$", size=20, fontweight="bold")
ax.set_xticks([])
ax.set_yticks(np.arange(N - 1) + 0.5)
ax.grid(linewidth=3, color="white")
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])
ax = plt.subplot(gs[0, 4])
ax.text(
    0.35,
    0.5,
    "=",
    horizontalalignment="center",
    verticalalignment="center",
    size=40,
    fontweight="bold",
)
ax.axis("off")
ax = plt.subplot(gs[0, 5])
ax.imshow(u[:, np.newaxis].T, cmap="rainbow")
ax.set_title(r"$u^T$", size=20, fontweight="bold")
ax.set_xticks(np.arange(M - 1) + 0.5)
ax.set_yticks([])
ax.grid(linewidth=3, color="white")
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])
ax = plt.subplot(gs[0, 6:8])
ax.imshow(Op.A.T, cmap="rainbow")
ax.set_title(r"$(Op^T*$", size=20, fontweight="bold")
ax.set_xticks(np.arange(N - 1) + 0.5)
ax.set_yticks(np.arange(M - 1) + 0.5)
ax.grid(linewidth=3, color="white")
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])
ax.axis("tight")
ax = plt.subplot(gs[0, 8])
ax.imshow(v[:, np.newaxis], cmap="rainbow")
ax.set_title(r"$v)$", size=20, fontweight="bold")
ax.set_xticks([])
ax.set_yticks(np.arange(N - 1) + 0.5)
ax.grid(linewidth=3, color="white")
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])
plt.tight_layout()
$(Op*$, $u)^T$, $v$, $u^T$, $(Op^T*$, $v)$

From now on, we can simply use the pylops.utils.dottest implementation of the dot-test and pass the operator we would like to validate, its size in the model and data spaces and optionally the tolerance we will be accepting for the dot-test to be considered succesfull. Finally we need to specify if our data or/and model vectors contain complex numbers using the complexflag parameter. While the dot-test will return True when succesfull and False otherwise, we can also ask to print its outcome putting the verb parameters to True.

N = 10
d = np.arange(N)
Dop = pylops.Diagonal(d)

_ = dottest(Dop, N, N, rtol=1e-6, complexflag=0, verb=True)
Dot test passed, v^H(Opu)=-6.91965436475742 - u^H(Op^Hv)=-6.91965436475742

We move now to a more complicated operator, the pylops.signalprocessing.FFT operator. We use once again the pylops.utils.dottest to verify its implementation and since we are dealing with a transform that can be applied to both real and complex array, we try different combinations using the complexflag input.

dt = 0.005
nt = 100
nfft = 2**10

FFTop = pylops.signalprocessing.FFT(
    dims=(nt,), nfft=nfft, sampling=dt, dtype=np.complex128
)
dottest(FFTop, nfft, nt, complexflag=2, verb=True)
_ = dottest(FFTop, nfft, nt, complexflag=3, verb=True)
Dot test passed, v^H(Opu)=(-17.13910350422714+3.8288339959112507j) - u^H(Op^Hv)=(-17.139103504227144+3.8288339959112503j)
Dot test passed, v^H(Opu)=(11.42080013111478-3.688151329246325j) - u^H(Op^Hv)=(11.420800131114781-3.6881513292463266j)

Total running time of the script: ( 0 minutes 0.266 seconds)

Gallery generated by Sphinx-Gallery

03. Solvers

This tutorial will guide you through the pylops.optimization module and show how to use various solvers that are included in the PyLops library.

The main idea here is to provide the user of PyLops with very high-level functionalities to quickly and easily set up and solve complex systems of linear equations as well as include regularization and/or preconditioning terms (all of those constructed by means of PyLops linear operators).

To make this tutorial more interesting, we will present a real life problem and show how the choice of the solver and regularization/preconditioning terms is vital in many circumstances to successfully retrieve an estimate of the model. The problem that we are going to consider is generally referred to as the data reconstruction problem and aims at reconstructing a regularly sampled signal of size \(M\) from \(N\) randomly selected samples:

\[\mathbf{y} = \mathbf{R} \mathbf{x}\]

where the restriction operator \(\mathbf{R}\) that selects the \(M\) elements from \(\mathbf{x}\) at random locations is implemented using pylops.Restriction, and

\[\mathbf{y}= [y_1, y_2,\ldots,y_N]^T, \qquad \mathbf{x}= [x_1, x_2,\ldots,x_M]^T, \qquad\]

with \(M \gg N\).

import matplotlib.pyplot as plt

# pylint: disable=C0103
import numpy as np

import pylops

plt.close("all")
np.random.seed(10)

Let’s first create the data in the frequency domain. The data is composed by the superposition of 3 sinusoids with different frequencies.

# Signal creation in frequency domain
ifreqs = [41, 25, 66]
amps = [1.0, 1.0, 1.0]
N = 200
nfft = 2**11
dt = 0.004
t = np.arange(N) * dt
f = np.fft.rfftfreq(nfft, dt)

FFTop = 10 * pylops.signalprocessing.FFT(N, nfft=nfft, real=True)

X = np.zeros(nfft // 2 + 1, dtype="complex128")
X[ifreqs] = amps
x = FFTop.H * X

fig, axs = plt.subplots(2, 1, figsize=(12, 8))
axs[0].plot(f, np.abs(X), "k", lw=2)
axs[0].set_xlim(0, 30)
axs[0].set_title("Data(frequency domain)")
axs[1].plot(t, x, "k", lw=2)
axs[1].set_title("Data(time domain)")
axs[1].axis("tight")
plt.tight_layout()
Data(frequency domain), Data(time domain)

We now define the locations at which the signal will be sampled.

# subsampling locations
perc_subsampling = 0.2
Nsub = int(np.round(N * perc_subsampling))

iava = np.sort(np.random.permutation(np.arange(N))[:Nsub])

# Create restriction operator
Rop = pylops.Restriction(N, iava, dtype="float64")

y = Rop * x
ymask = Rop.mask(x)

# Visualize data
fig = plt.figure(figsize=(12, 4))
plt.plot(t, x, "k", lw=3)
plt.plot(t, x, ".k", ms=20, label="all samples")
plt.plot(t, ymask, ".g", ms=15, label="available samples")
plt.legend()
plt.title("Data restriction")
plt.tight_layout()
Data restriction

To start let’s consider the simplest ‘solver’, i.e., least-square inversion without regularization. We aim here to minimize the following cost function:

\[J = \|\mathbf{y} - \mathbf{R} \mathbf{x}\|_2^2\]

Depending on the choice of the operator \(\mathbf{R}\), such problem can be solved using explicit matrix solvers as well as iterative solvers. In this case we will be using the latter approach (more specifically the scipy implementation of the LSQR solver - i.e., scipy.sparse.linalg.lsqr) as we do not want to explicitly create and invert a matrix. In most cases this will be the only viable approach as most of the large-scale optimization problems that we are interested to solve using PyLops do not lend naturally to the creation and inversion of explicit matrices.

This first solver can be very easily implemented using the / for PyLops operators, which will automatically call the scipy.sparse.linalg.lsqr with some default parameters.

xinv = Rop / y

We can also use pylops.optimization.leastsquares.regularized_inversion (without regularization term for now) and customize our solvers using kwargs.

xinv = pylops.optimization.leastsquares.regularized_inversion(
    Rop, y, [], **dict(damp=0, iter_lim=10, show=True)
)[0]
RegularizedInversion
-----------------------------------------------------------------
The Operator Op has 40 rows and 200 cols
Regs=[]
epsRs=[]
-----------------------------------------------------------------

LSQR            Least-squares solution of  Ax = b
The matrix A has 40 rows and 200 columns
damp = 0.00000000000000e+00   calc_var =        0
atol = 1.00e-06                 conlim = 1.00e+08
btol = 1.00e-06               iter_lim =       10

   Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
     0  0.00000e+00   2.658e+00  2.658e+00    1.0e+00  3.8e-01
     1  0.00000e+00   0.000e+00  0.000e+00    0.0e+00  0.0e+00   0.0e+00  0.0e+00

LSQR finished
Ax - b is small enough, given atol, btol

istop =       1   r1norm = 0.0e+00   anorm = 0.0e+00   arnorm = 0.0e+00
itn   =       1   r2norm = 0.0e+00   acond = 0.0e+00   xnorm  = 2.7e+00

Finally we can select a different starting guess from the null vector

xinv_fromx0 = pylops.optimization.leastsquares.regularized_inversion(
    Rop, y, [], x0=np.ones(N), **dict(damp=0, iter_lim=10, show=True)
)[0]
RegularizedInversion
-----------------------------------------------------------------
The Operator Op has 40 rows and 200 cols
Regs=[]
epsRs=[]
-----------------------------------------------------------------

LSQR            Least-squares solution of  Ax = b
The matrix A has 40 rows and 200 columns
damp = 0.00000000000000e+00   calc_var =        0
atol = 1.00e-06                 conlim = 1.00e+08
btol = 1.00e-06               iter_lim =       10

   Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
     0  0.00000e+00   6.737e+00  6.737e+00    1.0e+00  1.5e-01
     1  0.00000e+00   0.000e+00  0.000e+00    0.0e+00  0.0e+00   0.0e+00  0.0e+00

LSQR finished
Ax - b is small enough, given atol, btol

istop =       1   r1norm = 0.0e+00   anorm = 0.0e+00   arnorm = 0.0e+00
itn   =       1   r2norm = 0.0e+00   acond = 0.0e+00   xnorm  = 6.7e+00

The cost function above can be also expanded in terms of its normal equations

\[\mathbf{x}_{ne}= (\mathbf{R}^T \mathbf{R})^{-1} \mathbf{R}^T \mathbf{y}\]

The method pylops.optimization.leastsquares.normal_equations_inversion implements such system of equations explicitly and solves them using an iterative scheme suitable for square matrices (i.e., \(M=N\)).

While this approach may seem not very useful, we will soon see how regularization terms could be easily added to the normal equations using this method.

xne = pylops.optimization.leastsquares.normal_equations_inversion(Rop, y, [])[0]

Let’s now visualize the different inversion results

fig = plt.figure(figsize=(12, 4))
plt.plot(t, x, "k", lw=2, label="original")
plt.plot(t, xinv, "b", ms=10, label="inversion")
plt.plot(t, xinv_fromx0, "--r", ms=10, label="inversion from x0")
plt.plot(t, xne, "--g", ms=10, label="normal equations")
plt.legend()
plt.title("Data reconstruction without regularization")
plt.tight_layout()
Data reconstruction without regularization
Regularization

You may have noticed that none of the inversion has been successfull in recovering the original signal. This is a clear indication that the problem we are trying to solve is highly ill-posed and requires some prior knowledge from the user.

We will now see how to add prior information to the inverse process in the form of regularization (or preconditioning). This can be done in two different ways

Let’s start by regularizing the normal equations using a second derivative operator

\[\mathbf{x} = (\mathbf{R^TR}+\epsilon_\nabla^2\nabla^T\nabla)^{-1} \mathbf{R^Ty}\]
# Create regularization operator
D2op = pylops.SecondDerivative(N, dtype="float64")

# Regularized inversion
epsR = np.sqrt(0.1)
epsI = np.sqrt(1e-4)

xne = pylops.optimization.leastsquares.normal_equations_inversion(
    Rop, y, [D2op], epsI=epsI, epsRs=[epsR], **dict(maxiter=50)
)[0]

Note that in case we have access to a fast implementation for the chain of forward and adjoint for the regularization operator (i.e., \(\nabla^T\nabla\)), we can modify our call to pylops.optimization.leastsquares.normal_equations_inversion as follows:

ND2op = pylops.MatrixMult((D2op.H * D2op).tosparse())  # mimic fast D^T D

xne1 = pylops.optimization.leastsquares.normal_equations_inversion(
    Rop, y, [], NRegs=[ND2op], epsI=epsI, epsNRs=[epsR], **dict(maxiter=50)
)[0]

We can do the same while using pylops.optimization.leastsquares.regularized_inversion which solves the following augmented problem

\[\begin{split}\begin{bmatrix} \mathbf{R} \\ \epsilon_\nabla \nabla \end{bmatrix} \mathbf{x} = \begin{bmatrix} \mathbf{y} \\ 0 \end{bmatrix}\end{split}\]
xreg = pylops.optimization.leastsquares.regularized_inversion(
    Rop,
    y,
    [D2op],
    epsRs=[np.sqrt(0.1)],
    **dict(damp=np.sqrt(1e-4), iter_lim=50, show=0)
)[0]

We can also write a preconditioned problem, whose cost function is

\[J= \|\mathbf{y} - \mathbf{R} \mathbf{P} \mathbf{p}\|_2^2\]

where \(\mathbf{P}\) is the precondioned operator, \(\mathbf{p}\) is the projected model in the preconditioned space, and \(\mathbf{x}=\mathbf{P}\mathbf{p}\) is the model in the original model space we want to solve for. Note that a preconditioned problem converges much faster to its solution than its corresponding regularized problem. This can be done using the routine pylops.optimization.leastsquares.preconditioned_inversion.

# Create regularization operator
Sop = pylops.Smoothing1D(nsmooth=11, dims=[N], dtype="float64")

# Invert for interpolated signal
xprec = pylops.optimization.leastsquares.preconditioned_inversion(
    Rop, y, Sop, **dict(damp=np.sqrt(1e-9), iter_lim=20, show=0)
)[0]

Let’s finally visualize these solutions

# sphinx_gallery_thumbnail_number=4
fig = plt.figure(figsize=(12, 4))
plt.plot(t[iava], y, ".k", ms=20, label="available samples")
plt.plot(t, x, "k", lw=3, label="original")
plt.plot(t, xne, "b", lw=3, label="normal equations")
plt.plot(t, xne1, "--c", lw=3, label="normal equations (with direct D^T D)")
plt.plot(t, xreg, "-.r", lw=3, label="regularized")
plt.plot(t, xprec, "--g", lw=3, label="preconditioned equations")
plt.legend()
plt.title("Data reconstruction with regularization")

subax = fig.add_axes([0.7, 0.2, 0.15, 0.6])
subax.plot(t[iava], y, ".k", ms=20)
subax.plot(t, x, "k", lw=3)
subax.plot(t, xne, "b", lw=3)
subax.plot(t, xne1, "--c", lw=3)
subax.plot(t, xreg, "-.r", lw=3)
subax.plot(t, xprec, "--g", lw=3)
subax.set_xlim(0.05, 0.3)
plt.tight_layout()
Data reconstruction with regularization

Much better estimates! We have seen here how regularization and/or preconditioning can be vital to succesfully solve some ill-posed inverse problems.

We have however so far only considered solvers that can include additional norm-2 regularization terms. A very active area of research is that of sparsity-promoting solvers (also sometimes referred to as compressive sensing): the regularization term added to the cost function to minimize has norm-p (\(p \le 1\)) and the problem is generally recasted by considering the model to be sparse in some domain. We can follow this philosophy as our signal to invert was actually created as superposition of 3 sinusoids (i.e., three spikes in the Fourier domain). Our new cost function is:

\[J_1 = \|\mathbf{y} - \mathbf{R} \mathbf{F} \mathbf{p}\|_2^2 + \epsilon \|\mathbf{p}\|_1\]

where \(\mathbf{F}\) is the FFT operator. We will thus use the pylops.optimization.sparsity.ista and pylops.optimization.sparsity.fista solvers to estimate our input signal.

pista, niteri, costi = pylops.optimization.sparsity.ista(
    Rop * FFTop.H,
    y,
    niter=1000,
    eps=0.1,
    tol=1e-7,
)
xista = FFTop.H * pista

pfista, niterf, costf = pylops.optimization.sparsity.fista(
    Rop * FFTop.H,
    y,
    niter=1000,
    eps=0.1,
    tol=1e-7,
)
xfista = FFTop.H * pfista

fig, axs = plt.subplots(2, 1, figsize=(12, 8))
fig.suptitle("Data reconstruction with sparsity", fontsize=14, fontweight="bold", y=0.9)
axs[0].plot(f, np.abs(X), "k", lw=3)
axs[0].plot(f, np.abs(pista), "--r", lw=3)
axs[0].plot(f, np.abs(pfista), "--g", lw=3)
axs[0].set_xlim(0, 30)
axs[0].set_title("Frequency domain")
axs[1].plot(t[iava], y, ".k", ms=20, label="available samples")
axs[1].plot(t, x, "k", lw=3, label="original")
axs[1].plot(t, xista, "--r", lw=3, label="ISTA")
axs[1].plot(t, xfista, "--g", lw=3, label="FISTA")
axs[1].set_title("Time domain")
axs[1].axis("tight")
axs[1].legend()
plt.tight_layout()
plt.subplots_adjust(top=0.8)

fig, ax = plt.subplots(1, 1, figsize=(12, 3))
ax.semilogy(costi, "r", lw=2, label="ISTA")
ax.semilogy(costf, "g", lw=2, label="FISTA")
ax.set_title("Cost functions", size=15, fontweight="bold")
ax.set_xlabel("Iteration")
ax.legend()
ax.grid(True)
plt.tight_layout()
  • Data reconstruction with sparsity, Frequency domain, Time domain
  • Cost functions

As you can see, changing parametrization of the model and imposing sparsity in the Fourier domain has given an extra improvement to our ability of recovering the underlying densely sampled input signal. Moreover, FISTA converges much faster than ISTA as expected and should be preferred when using sparse solvers.

Finally we consider a slightly different cost function (note that in this case we try to solve a constrained problem):

\[J_1 = \|\mathbf{p}\|_1 \quad \text{subject to} \quad \|\mathbf{y} - \mathbf{R} \mathbf{F} \mathbf{p}\|\]

A very popular solver to solve such kind of cost function is called spgl1 and can be accessed via pylops.optimization.sparsity.spgl1.

xspgl1, pspgl1, info = pylops.optimization.sparsity.spgl1(
    Rop, y, SOp=FFTop, tau=3, iter_lim=200
)

fig, axs = plt.subplots(2, 1, figsize=(12, 8))
fig.suptitle("Data reconstruction with SPGL1", fontsize=14, fontweight="bold", y=0.9)
axs[0].plot(f, np.abs(X), "k", lw=3)
axs[0].plot(f, np.abs(pspgl1), "--m", lw=3)
axs[0].set_xlim(0, 30)
axs[0].set_title("Frequency domain")
axs[1].plot(t[iava], y, ".k", ms=20, label="available samples")
axs[1].plot(t, x, "k", lw=3, label="original")
axs[1].plot(t, xspgl1, "--m", lw=3, label="SPGL1")
axs[1].set_title("Time domain")
axs[1].axis("tight")
axs[1].legend()
plt.tight_layout()
plt.subplots_adjust(top=0.8)

fig, ax = plt.subplots(1, 1, figsize=(12, 3))
ax.semilogy(info["rnorm2"], "k", lw=2, label="ISTA")
ax.set_title("Cost functions", size=15, fontweight="bold")
ax.set_xlabel("Iteration")
ax.legend()
ax.grid(True)
plt.tight_layout()
  • Data reconstruction with SPGL1, Frequency domain, Time domain
  • Cost functions

Total running time of the script: ( 0 minutes 3.996 seconds)

Gallery generated by Sphinx-Gallery

03. Solvers (Advanced)

This is a follow up tutorial to the 03. Solvers tutorial. The same example will be considered, however we will showcase how to use the class-based version of our solvers (introduced in PyLops v2).

First of all, when shall you use class-based solvers over function-based ones? The answer is simple, every time you feel you would have like to have more flexibility when using one PyLops function-based solvers.

In fact, a function-based solver in PyLops v2 is nothing more than a thin wrapper over its class-based equivalent, which generally performs the following steps:

  • solver initialization

  • setup

  • run (by calling multiple times step)

  • finalize

The nice thing about class-based solvers is that i) a user can manually orchestrate these steps and do anything in between them; ii) a user can create a class-based pylops.optimization.callback.Callbacks object and define a set of callbacks that will be run pre and post setup, step and run. One example of how such callbacks can be handy to track evolving variables in the solver can be found in Linear Regression.

In the following we will leverage the very same mechanism to keep track of a number of metrics using the predefined pylops.optimization.callback.MetricsCallback callback. Finally we show how to create a customized callback that can track the percentage change of the solution and residual. This is of course just an example, we expect users will find different use cases based on the problem at hand.

import matplotlib.pyplot as plt

# pylint: disable=C0103
import numpy as np

import pylops

plt.close("all")
np.random.seed(10)

Let’s first create the data in the frequency domain. The data is composed by the superposition of 3 sinusoids with different frequencies.

# Signal creation in frequency domain
ifreqs = [41, 25, 66]
amps = [1.0, 1.0, 1.0]
N = 200
nfft = 2**11
dt = 0.004
t = np.arange(N) * dt
f = np.fft.rfftfreq(nfft, dt)

FFTop = 10 * pylops.signalprocessing.FFT(N, nfft=nfft, real=True)

X = np.zeros(nfft // 2 + 1, dtype="complex128")
X[ifreqs] = amps
x = FFTop.H * X

We now define the locations at which the signal will be sampled.

# subsampling locations
perc_subsampling = 0.2
Nsub = int(np.round(N * perc_subsampling))

iava = np.sort(np.random.permutation(np.arange(N))[:Nsub])

# Create restriction operator
Rop = pylops.Restriction(N, iava, dtype="float64")

y = Rop * x
ymask = Rop.mask(x)

Let’s now solve the interpolation problem using the pylops.optimization.sparsity.ISTA class-based solver.

cb = pylops.optimization.callback.MetricsCallback(x, FFTop.H)

istasolve = pylops.optimization.sparsity.ISTA(
    Rop * FFTop.H,
    callbacks=[
        cb,
    ],
)
pista, niteri, costi = istasolve.solve(y, niter=1000, eps=0.1, tol=1e-7)
xista = FFTop.H * pista

fig, axs = plt.subplots(1, 4, figsize=(16, 3))
for i, metric in enumerate(["mae", "mse", "snr", "psnr"]):
    axs[i].plot(cb.metrics[metric], "k", lw=2)
    axs[i].set_title(metric)
plt.tight_layout()
mae, mse, snr, psnr

Finally, we show how we can also define customized callbacks. What we are really interested in here is to store the first residual norm once the setup of the solver is over, and repeat the same after each step (using the previous estimate to compute the percentage change). And, we do the same for the solution norm.

class CallbackISTA(pylops.optimization.callback.Callbacks):
    def __init__(self):
        self.res_perc = []
        self.x_perc = []

    def on_setup_end(self, solver, x):
        self.x = x
        if x is not None:
            self.rec = solver.Op @ x - solver.y
        else:
            self.rec = None

    def on_step_end(self, solver, x):
        self.xold = self.x
        self.x = x
        self.recold = self.rec
        self.rec = solver.Op @ x - solver.y
        if self.xold is not None:
            self.x_perc.append(
                100 * np.linalg.norm(self.x - self.xold) / np.linalg.norm(self.xold)
            )
            self.res_perc.append(
                100
                * np.linalg.norm(self.rec - self.recold)
                / np.linalg.norm(self.recold)
            )

    def on_run_end(self, solver, x):
        # remove first percentage
        self.x_perc = np.array(self.x_perc[1:])
        self.res_perc = np.array(self.res_perc[1:])


cb = CallbackISTA()
istasolve = pylops.optimization.sparsity.ISTA(
    Rop * FFTop.H,
    callbacks=[
        cb,
    ],
)
pista, niteri, costi = istasolve.solve(y, niter=1000, eps=0.1, tol=1e-7)
xista = FFTop.H * pista

cbf = CallbackISTA()
fistasolve = pylops.optimization.sparsity.FISTA(
    Rop * FFTop.H,
    callbacks=[
        cbf,
    ],
)
pfista, niterf, costf = fistasolve.solve(y, niter=1000, eps=0.1, tol=1e-7)
xfista = FFTop.H * pfista

fig, axs = plt.subplots(2, 1, figsize=(12, 8))
fig.suptitle("Data reconstruction with sparsity", fontsize=14, fontweight="bold", y=0.9)
axs[0].plot(f, np.abs(X), "k", lw=3)
axs[0].plot(f, np.abs(pista), "--r", lw=3)
axs[0].plot(f, np.abs(pfista), "--g", lw=3)
axs[0].set_xlim(0, 30)
axs[0].set_title("Frequency domain")
axs[1].plot(t[iava], y, ".k", ms=20, label="available samples")
axs[1].plot(t, x, "k", lw=3, label="original")
axs[1].plot(t, xista, "--r", lw=3, label="ISTA")
axs[1].plot(t, xfista, "--g", lw=3, label="FISTA")
axs[1].set_title("Time domain")
axs[1].axis("tight")
axs[1].legend()
plt.tight_layout()
plt.subplots_adjust(top=0.8)

fig, axs = plt.subplots(2, 1, figsize=(12, 8))
fig.suptitle("Norms history", fontsize=14, fontweight="bold", y=0.9)
axs[0].semilogy(cb.res_perc, "r", lw=3)
axs[0].semilogy(cbf.res_perc, "g", lw=3)
axs[0].set_title("Residual percentage change")
axs[1].semilogy(cb.x_perc, "r", lw=3, label="ISTA")
axs[1].semilogy(cbf.x_perc, "g", lw=3, label="FISTA")
axs[1].set_title("Solution percentage change")
axs[1].legend()
plt.tight_layout()
plt.subplots_adjust(top=0.8)
  • Data reconstruction with sparsity, Frequency domain, Time domain
  • Norms history, Residual percentage change, Solution percentage change

Total running time of the script: ( 0 minutes 3.739 seconds)

Gallery generated by Sphinx-Gallery

04. Bayesian Inversion

This tutorial focuses on Bayesian inversion, a special type of inverse problem that aims at incorporating prior information in terms of model and data probabilities in the inversion process.

In this case we will be dealing with the same problem that we discussed in 03. Solvers, but instead of defining ad-hoc regularization or preconditioning terms we parametrize and model our input signal in the frequency domain in a probabilistic fashion: the central frequency, amplitude and phase of the three sinusoids have gaussian distributions as follows:

\[X(f) = \sum_{i=1}^3 a_i e^{j \phi_i} \delta(f - f_i)\]

where \(f_i \sim N(f_{0,i}, \sigma_{f,i})\), \(a_i \sim N(a_{0,i}, \sigma_{a,i})\), and \(\phi_i \sim N(\phi_{0,i}, \sigma_{\phi,i})\).

Based on the above definition, we construct some prior models in the frequency domain, convert each of them to the time domain and use such an ensemble to estimate the prior mean \(\mu_\mathbf{x}\) and model covariance \(\mathbf{C_x}\).

We then create our data by sampling the true signal at certain locations and solve the resconstruction problem within a Bayesian framework. Since we are assuming gaussianity in our priors, the equation to obtain the posterion mean can be derived analytically:

\[\mathbf{x} = \mathbf{x_0} + \mathbf{C}_x \mathbf{R}^T (\mathbf{R} \mathbf{C}_x \mathbf{R}^T + \mathbf{C}_y)^{-1} (\mathbf{y} - \mathbf{R} \mathbf{x_0})\]
import matplotlib.pyplot as plt

# sphinx_gallery_thumbnail_number = 2
import numpy as np
from scipy.sparse.linalg import lsqr

import pylops

plt.close("all")
np.random.seed(10)

Let’s start by creating our true model and prior realizations

def prior_realization(f0, a0, phi0, sigmaf, sigmaa, sigmaphi, dt, nt, nfft):
    """Create realization from prior mean and std for amplitude, frequency and
    phase
    """
    f = np.fft.rfftfreq(nfft, dt)
    df = f[1] - f[0]
    ifreqs = [int(np.random.normal(f, sigma) / df) for f, sigma in zip(f0, sigmaf)]
    amps = [np.random.normal(a, sigma) for a, sigma in zip(a0, sigmaa)]
    phis = [np.random.normal(phi, sigma) for phi, sigma in zip(phi0, sigmaphi)]

    # input signal in frequency domain
    X = np.zeros(nfft // 2 + 1, dtype="complex128")
    X[ifreqs] = (
        np.array(amps).squeeze() * np.exp(1j * np.deg2rad(np.array(phis))).squeeze()
    )

    # input signal in time domain
    FFTop = pylops.signalprocessing.FFT(nt, nfft=nfft, real=True)
    x = FFTop.H * X
    return x


# Priors
nreals = 100
f0 = [5, 3, 8]
sigmaf = [0.5, 1.0, 0.6]
a0 = [1.0, 1.0, 1.0]
sigmaa = [0.1, 0.5, 0.6]
phi0 = [-90.0, 0.0, 0.0]
sigmaphi = [0.1, 0.2, 0.4]
sigmad = 1e-2

# Prior models
nt = 200
nfft = 2**11
dt = 0.004
t = np.arange(nt) * dt
xs = np.array(
    [
        prior_realization(f0, a0, phi0, sigmaf, sigmaa, sigmaphi, dt, nt, nfft)
        for _ in range(nreals)
    ]
)

# True model (taken as one possible realization)
x = prior_realization(f0, a0, phi0, [0, 0, 0], [0, 0, 0], [0, 0, 0], dt, nt, nfft)

We have now a set of prior models in time domain. We can easily use sample statistics to estimate the prior mean and covariance. For the covariance, we perform a second step where we average values around the main diagonal for each row and find a smooth, compact filter that we use to define a convolution linear operator that mimics the action of the covariance matrix on a vector

x0 = np.average(xs, axis=0)
Cm = ((xs - x0).T @ (xs - x0)) / nreals

N = 30  # lenght of decorrelation
diags = np.array([Cm[i, i - N : i + N + 1] for i in range(N, nt - N)])
diag_ave = np.average(diags, axis=0)
# add a taper at the end to avoid edge effects
diag_ave *= np.hamming(2 * N + 1)

fig, ax = plt.subplots(1, 1, figsize=(12, 4))
ax.plot(t, xs.T, "r", lw=1)
ax.plot(t, x0, "g", lw=4)
ax.plot(t, x, "k", lw=4)
ax.set_title("Prior realizations and mean")
ax.set_xlim(0, 0.8)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
im = ax1.imshow(
    Cm, interpolation="nearest", cmap="seismic", extent=(t[0], t[-1], t[-1], t[0])
)
ax1.set_title(r"$\mathbf{C}_m^{prior}$")
ax1.axis("tight")
ax2.plot(np.arange(-N, N + 1) * dt, diags.T, "--r", lw=1)
ax2.plot(np.arange(-N, N + 1) * dt, diag_ave, "k", lw=4)
ax2.set_title("Averaged covariance 'filter'")
plt.tight_layout()
  • Prior realizations and mean
  • $\mathbf{C}_m^{prior}$, Averaged covariance 'filter'

Let’s define now the sampling operator as well as create our covariance matrices in terms of linear operators. This may not be strictly necessary here but shows how even Bayesian-type of inversion can very easily scale to large model and data spaces.

# Sampling operator
perc_subsampling = 0.2
ntsub = int(np.round(nt * perc_subsampling))
iava = np.sort(np.random.permutation(np.arange(nt))[:ntsub])
iava[-1] = nt - 1  # assume we have the last sample to avoid instability
Rop = pylops.Restriction(nt, iava, dtype="float64")

# Covariance operators
Cm_op = pylops.signalprocessing.Convolve1D(nt, diag_ave, offset=N)
Cd_op = sigmad**2 * pylops.Identity(ntsub)

We model now our data and add noise that respects our prior definition

n = np.random.normal(0, sigmad, nt)
y = Rop * x
yn = Rop * (x + n)
ymask = Rop.mask(x)
ynmask = Rop.mask(x + n)

First we apply the Bayesian inversion equation

xbayes = x0 + Cm_op * Rop.H * (
    lsqr(Rop * Cm_op * Rop.H + Cd_op, yn - Rop * x0, iter_lim=400)[0]
)

# Visualize
fig, ax = plt.subplots(1, 1, figsize=(12, 5))
ax.plot(t, x, "k", lw=6, label="true")
ax.plot(t, ymask, ".k", ms=25, label="available samples")
ax.plot(t, ynmask, ".r", ms=25, label="available noisy samples")
ax.plot(t, xbayes, "r", lw=3, label="bayesian inverse")
ax.legend()
ax.set_title("Signal")
ax.set_xlim(0, 0.8)
plt.tight_layout()
Signal

So far we have been able to estimate our posterion mean. What about its uncertainties (i.e., posterion covariance)?

In real-life applications it is very difficult (if not impossible) to directly compute the posterior covariance matrix. It is much more useful to create a set of models that sample the posterion probability. We can do that by solving our problem several times using different prior realizations as starting guesses:

xpost = [
    x0
    + Cm_op
    * Rop.H
    * (lsqr(Rop * Cm_op * Rop.H + Cd_op, yn - Rop * x0, iter_lim=400)[0])
    for x0 in xs[:30]
]
xpost = np.array(xpost)

x0post = np.average(xpost, axis=0)
Cm_post = ((xpost - x0post).T @ (xpost - x0post)) / nreals

# Visualize
fig, ax = plt.subplots(1, 1, figsize=(12, 5))
ax.plot(t, x, "k", lw=6, label="true")
ax.plot(t, xpost.T, "--r", lw=1)
ax.plot(t, x0post, "r", lw=3, label="bayesian inverse")
ax.plot(t, ymask, ".k", ms=25, label="available samples")
ax.plot(t, ynmask, ".r", ms=25, label="available noisy samples")
ax.legend()
ax.set_title("Signal")
ax.set_xlim(0, 0.8)

fig, ax = plt.subplots(1, 1, figsize=(5, 4))
im = ax.imshow(
    Cm_post, interpolation="nearest", cmap="seismic", extent=(t[0], t[-1], t[-1], t[0])
)
ax.set_title(r"$\mathbf{C}_m^{posterior}$")
ax.axis("tight")
plt.tight_layout()
  • Signal
  • $\mathbf{C}_m^{posterior}$

Note that here we have been able to compute a sample posterior covariance from its estimated samples. By displaying it we can see how both the overall variances and the correlation between different parameters have become narrower compared to their prior counterparts.

Total running time of the script: ( 0 minutes 1.659 seconds)

Gallery generated by Sphinx-Gallery

05. Image deblurring

Deblurring is the process of removing blurring effects from images, caused for example by defocus aberration or motion blur.

In forward mode, such blurring effect is typically modelled as a 2-dimensional convolution between the so-called point spread function and a target sharp input image, where the sharp input image (which has to be recovered) is unknown and the point-spread function can be either known or unknown.

In this tutorial, an example of 2d blurring and deblurring will be shown using the pylops.signalprocessing.Convolve2D operator assuming knowledge of the point-spread function.

import matplotlib.pyplot as plt
import numpy as np

import pylops

Let’s start by importing a 2d image and defining the blurring operator

im = np.load("../testdata/python.npy")[::5, ::5, 0]

Nz, Nx = im.shape

# Blurring guassian operator
nh = [15, 25]
hz = np.exp(-0.1 * np.linspace(-(nh[0] // 2), nh[0] // 2, nh[0]) ** 2)
hx = np.exp(-0.03 * np.linspace(-(nh[1] // 2), nh[1] // 2, nh[1]) ** 2)
hz /= np.trapz(hz)  # normalize the integral to 1
hx /= np.trapz(hx)  # normalize the integral to 1
h = hz[:, np.newaxis] * hx[np.newaxis, :]

fig, ax = plt.subplots(1, 1, figsize=(5, 3))
him = ax.imshow(h)
ax.set_title("Blurring operator")
fig.colorbar(him, ax=ax)
ax.axis("tight")

Cop = pylops.signalprocessing.Convolve2D(
    (Nz, Nx), h=h, offset=(nh[0] // 2, nh[1] // 2), dtype="float32"
)
Blurring operator

We first apply the blurring operator to the sharp image. We then try to recover the sharp input image by inverting the convolution operator from the blurred image. Note that when we perform inversion without any regularization, the deblurred image will show some ringing due to the instabilities of the inverse process. Using a L1 solver with a DWT preconditioner or TV regularization allows to recover sharper contrasts.

imblur = Cop * im

imdeblur = pylops.optimization.leastsquares.normal_equations_inversion(
    Cop, imblur.ravel(), None, maxiter=50  # solvers need 1D arrays
)[0]
imdeblur = imdeblur.reshape(Cop.dims)

Wop = pylops.signalprocessing.DWT2D((Nz, Nx), wavelet="haar", level=3)
Dop = [
    pylops.FirstDerivative((Nz, Nx), axis=0, edge=False),
    pylops.FirstDerivative((Nz, Nx), axis=1, edge=False),
]
DWop = Dop + [Wop]

imdeblurfista = pylops.optimization.sparsity.fista(
    Cop * Wop.H, imblur.ravel(), eps=1e-1, niter=100
)[0]
imdeblurfista = imdeblurfista.reshape((Cop * Wop.H).dims)
imdeblurfista = Wop.H * imdeblurfista

imdeblurtv = pylops.optimization.sparsity.splitbregman(
    Cop,
    imblur.ravel(),
    Dop,
    niter_outer=10,
    niter_inner=5,
    mu=1.5,
    epsRL1s=[2e0, 2e0],
    tol=1e-4,
    tau=1.0,
    show=False,
    **dict(iter_lim=5, damp=1e-4)
)[0]
imdeblurtv = imdeblurtv.reshape(Cop.dims)

imdeblurtv1 = pylops.optimization.sparsity.splitbregman(
    Cop,
    imblur.ravel(),
    DWop,
    niter_outer=10,
    niter_inner=5,
    mu=1.5,
    epsRL1s=[1e0, 1e0, 1e0],
    tol=1e-4,
    tau=1.0,
    show=False,
    **dict(iter_lim=5, damp=1e-4)
)[0]
imdeblurtv1 = imdeblurtv1.reshape(Cop.dims)

Finally we visualize the original, blurred, and recovered images.

# sphinx_gallery_thumbnail_number = 2
fig = plt.figure(figsize=(12, 6))
fig.suptitle("Deblurring", fontsize=14, fontweight="bold", y=0.95)
ax1 = plt.subplot2grid((2, 5), (0, 0))
ax2 = plt.subplot2grid((2, 5), (0, 1))
ax3 = plt.subplot2grid((2, 5), (0, 2))
ax4 = plt.subplot2grid((2, 5), (1, 0))
ax5 = plt.subplot2grid((2, 5), (1, 1))
ax6 = plt.subplot2grid((2, 5), (1, 2))
ax7 = plt.subplot2grid((2, 5), (0, 3), colspan=2)
ax8 = plt.subplot2grid((2, 5), (1, 3), colspan=2)
ax1.imshow(im, cmap="viridis", vmin=0, vmax=250)
ax1.axis("tight")
ax1.set_title("Original")
ax2.imshow(imblur, cmap="viridis", vmin=0, vmax=250)
ax2.axis("tight")
ax2.set_title("Blurred")
ax3.imshow(imdeblur, cmap="viridis", vmin=0, vmax=250)
ax3.axis("tight")
ax3.set_title("Deblurred")
ax4.imshow(imdeblurfista, cmap="viridis", vmin=0, vmax=250)
ax4.axis("tight")
ax4.set_title("FISTA deblurred")
ax5.imshow(imdeblurtv, cmap="viridis", vmin=0, vmax=250)
ax5.axis("tight")
ax5.set_title("TV deblurred")
ax6.imshow(imdeblurtv1, cmap="viridis", vmin=0, vmax=250)
ax6.axis("tight")
ax6.set_title("TV+Haar deblurred")
ax7.plot(im[Nz // 2], "k")
ax7.plot(imblur[Nz // 2], "--r")
ax7.plot(imdeblur[Nz // 2], "--b")
ax7.plot(imdeblurfista[Nz // 2], "--g")
ax7.plot(imdeblurtv[Nz // 2], "--m")
ax7.plot(imdeblurtv1[Nz // 2], "--y")
ax7.axis("tight")
ax7.set_title("Horizontal section")
ax8.plot(im[:, Nx // 2], "k", label="Original")
ax8.plot(imblur[:, Nx // 2], "--r", label="Blurred")
ax8.plot(imdeblur[:, Nx // 2], "--b", label="Deblurred")
ax8.plot(imdeblurfista[:, Nx // 2], "--g", label="FISTA deblurred")
ax8.plot(imdeblurtv[:, Nx // 2], "--m", label="TV deblurred")
ax8.plot(imdeblurtv1[:, Nx // 2], "--y", label="TV+Haar deblurred")
ax8.axis("tight")
ax8.set_title("Vertical section")
ax8.legend(loc=5, fontsize="small")
plt.tight_layout()
plt.subplots_adjust(top=0.8)
Deblurring, Original, Blurred, Deblurred, FISTA deblurred, TV deblurred, TV+Haar deblurred, Horizontal section, Vertical section

Total running time of the script: ( 0 minutes 6.555 seconds)

Gallery generated by Sphinx-Gallery

06. 2D Interpolation

In the mathematical field of numerical analysis, interpolation is the problem of constructing new data points within the range of a discrete set of known data points. In signal and image processing, the data may be recorded at irregular locations and it is often required to regularize the data into a regular grid.

In this tutorial, an example of 2d interpolation of an image is carried out using a combination of PyLops operators (pylops.Restriction and pylops.Laplacian) and the pylops.optimization module.

Mathematically speaking, if we want to interpolate a signal using the theory of inverse problems, we can define the following forward problem:

\[\mathbf{y} = \mathbf{R} \mathbf{x}\]

where the restriction operator \(\mathbf{R}\) selects \(M\) elements from the regularly sampled signal \(\mathbf{x}\) at random locations. The input and output signals are:

\[\mathbf{y}= [y_1, y_2,\ldots,y_N]^T, \qquad \mathbf{x}= [x_1, x_2,\ldots,x_M]^T, \qquad\]

with \(M \gg N\).

import matplotlib.pyplot as plt
import numpy as np

import pylops

plt.close("all")
np.random.seed(0)

To start we import a 2d image and define our restriction operator to irregularly and randomly sample the image for 30% of the entire grid

im = np.load("../testdata/python.npy")[:, :, 0]

Nz, Nx = im.shape
N = Nz * Nx

# Subsample signal
perc_subsampling = 0.2

Nsub2d = int(np.round(N * perc_subsampling))
iava = np.sort(np.random.permutation(np.arange(N))[:Nsub2d])

# Create operators and data
Rop = pylops.Restriction(N, iava, dtype="float64")
D2op = pylops.Laplacian((Nz, Nx), weights=(1, 1), dtype="float64")

x = im.ravel()
y = Rop * x
y1 = Rop.mask(x)

We will now use two different routines from our optimization toolbox to estimate our original image in the regular grid.

xcg_reg_lop = pylops.optimization.leastsquares.normal_equations_inversion(
    Rop, y, [D2op], epsRs=[np.sqrt(0.1)], **dict(maxiter=200)
)[0]

# Invert for interpolated signal, lsqrt
xlsqr_reg_lop = pylops.optimization.leastsquares.regularized_inversion(
    Rop,
    y,
    [D2op],
    epsRs=[np.sqrt(0.1)],
    **dict(damp=0, iter_lim=200, show=0),
)[0]

# Reshape estimated images
im_sampled = y1.reshape((Nz, Nx))
im_rec_lap_cg = xcg_reg_lop.reshape((Nz, Nx))
im_rec_lap_lsqr = xlsqr_reg_lop.reshape((Nz, Nx))

Finally we visualize the original image, the reconstructed images and their error

fig, axs = plt.subplots(1, 4, figsize=(12, 4))
fig.suptitle("Data reconstruction - normal eqs", fontsize=14, fontweight="bold", y=0.95)
axs[0].imshow(im, cmap="viridis", vmin=0, vmax=250)
axs[0].axis("tight")
axs[0].set_title("Original")
axs[1].imshow(im_sampled.data, cmap="viridis", vmin=0, vmax=250)
axs[1].axis("tight")
axs[1].set_title("Sampled")
axs[2].imshow(im_rec_lap_cg, cmap="viridis", vmin=0, vmax=250)
axs[2].axis("tight")
axs[2].set_title("2D Regularization")
axs[3].imshow(im - im_rec_lap_cg, cmap="gray", vmin=-80, vmax=80)
axs[3].axis("tight")
axs[3].set_title("2D Regularization Error")
plt.tight_layout()
plt.subplots_adjust(top=0.8)

fig, axs = plt.subplots(1, 4, figsize=(12, 4))
fig.suptitle(
    "Data reconstruction - regularized eqs", fontsize=14, fontweight="bold", y=0.95
)
axs[0].imshow(im, cmap="viridis", vmin=0, vmax=250)
axs[0].axis("tight")
axs[0].set_title("Original")
axs[1].imshow(im_sampled.data, cmap="viridis", vmin=0, vmax=250)
axs[1].axis("tight")
axs[1].set_title("Sampled")
axs[2].imshow(im_rec_lap_lsqr, cmap="viridis", vmin=0, vmax=250)
axs[2].axis("tight")
axs[2].set_title("2D Regularization")
axs[3].imshow(im - im_rec_lap_lsqr, cmap="gray", vmin=-80, vmax=80)
axs[3].axis("tight")
axs[3].set_title("2D Regularization Error")
plt.tight_layout()
plt.subplots_adjust(top=0.8)
  • Data reconstruction - normal eqs, Original, Sampled, 2D Regularization, 2D Regularization Error
  • Data reconstruction - regularized eqs, Original, Sampled, 2D Regularization, 2D Regularization Error

Total running time of the script: ( 0 minutes 17.005 seconds)

Gallery generated by Sphinx-Gallery

07. Post-stack inversion

Estimating subsurface properties from band-limited seismic data represents an important task for geophysical subsurface characterization.

In this tutorial, the pylops.avo.poststack.PoststackLinearModelling operator is used for modelling of both 1d and 2d synthetic post-stack seismic data from a profile or 2d model of the subsurface acoustic impedence.

\[d(t, \theta=0) = \frac{1}{2} w(t) * \frac{\mathrm{d}\ln \text{AI}(t)}{\mathrm{d}t}\]

where \(\text{AI}(t)\) is the acoustic impedance profile and \(w(t)\) is the time domain seismic wavelet. In compact form:

\[\mathbf{d}= \mathbf{W} \mathbf{D} \mathbf{ai}\]

where \(\mathbf{W}\) is a convolution operator, \(\mathbf{D}\) is a first derivative operator, and \(\mathbf{ai}\) is the input model. Subsequently the acoustic impedance model is estimated via the pylops.avo.poststack.PoststackInversion module. A two-steps inversion strategy is finally presented to deal with the case of noisy data.

import matplotlib.pyplot as plt

# sphinx_gallery_thumbnail_number = 4
import numpy as np
from scipy.signal import filtfilt

import pylops
from pylops.utils.wavelets import ricker

plt.close("all")
np.random.seed(10)

Let’s start with a 1d example. A synthetic profile of acoustic impedance is created and data is modelled using both the dense and linear operator version of pylops.avo.poststack.PoststackLinearModelling operator.

# model
nt0 = 301
dt0 = 0.004
t0 = np.arange(nt0) * dt0
vp = 1200 + np.arange(nt0) + filtfilt(np.ones(5) / 5.0, 1, np.random.normal(0, 80, nt0))
rho = 1000 + vp + filtfilt(np.ones(5) / 5.0, 1, np.random.normal(0, 30, nt0))
vp[131:] += 500
rho[131:] += 100
m = np.log(vp * rho)

# smooth model
nsmooth = 100
mback = filtfilt(np.ones(nsmooth) / float(nsmooth), 1, m)

# wavelet
ntwav = 41
wav, twav, wavc = ricker(t0[: ntwav // 2 + 1], 20)

# dense operator
PPop_dense = pylops.avo.poststack.PoststackLinearModelling(
    wav / 2, nt0=nt0, explicit=True
)

# lop operator
PPop = pylops.avo.poststack.PoststackLinearModelling(wav / 2, nt0=nt0)

# data
d_dense = PPop_dense * m.ravel()
d = PPop * m

# add noise
dn_dense = d_dense + np.random.normal(0, 2e-2, d_dense.shape)

We can now estimate the acoustic profile from band-limited data using either the dense operator or linear operator.

# solve dense
minv_dense = pylops.avo.poststack.PoststackInversion(
    d, wav / 2, m0=mback, explicit=True, simultaneous=False
)[0]

# solve lop
minv = pylops.avo.poststack.PoststackInversion(
    d_dense,
    wav / 2,
    m0=mback,
    explicit=False,
    simultaneous=False,
    **dict(iter_lim=2000)
)[0]

# solve noisy
mn = pylops.avo.poststack.PoststackInversion(
    dn_dense, wav / 2, m0=mback, explicit=True, epsR=1e0, **dict(damp=1e-1)
)[0]

fig, axs = plt.subplots(1, 2, figsize=(6, 7), sharey=True)
axs[0].plot(d_dense, t0, "k", lw=4, label="Dense")
axs[0].plot(d, t0, "--r", lw=2, label="Lop")
axs[0].plot(dn_dense, t0, "-.g", lw=2, label="Noisy")
axs[0].set_title("Data")
axs[0].invert_yaxis()
axs[0].axis("tight")
axs[0].legend(loc=1)
axs[1].plot(m, t0, "k", lw=4, label="True")
axs[1].plot(mback, t0, "--b", lw=4, label="Back")
axs[1].plot(minv_dense, t0, "--m", lw=2, label="Inv Dense")
axs[1].plot(minv, t0, "--r", lw=2, label="Inv Lop")
axs[1].plot(mn, t0, "--g", lw=2, label="Inv Noisy")
axs[1].set_title("Model")
axs[1].axis("tight")
axs[1].legend(loc=1)
plt.tight_layout()
Data, Model

We see how inverting a dense matrix is in this case faster than solving for the linear operator (a good estimate of the model is in fact obtained only after 2000 iterations of lsqr). Nevertheless, having a linear operator is useful when we deal with larger dimensions (2d or 3d) and we want to couple our modelling operator with different types of spatial regularizations or preconditioning.

Before we move onto a 2d example, let’s consider the case of non-stationary wavelet and see how we can easily use the same routines in this case

# wavelet
ntwav = 41
f0s = np.flip(np.arange(nt0) * 0.05 + 3)
wavs = np.array([ricker(t0[:ntwav], f0)[0] for f0 in f0s])
wavc = np.argmax(wavs[0])

plt.figure(figsize=(5, 4))
plt.imshow(wavs.T, cmap="gray", extent=(t0[0], t0[-1], t0[ntwav], -t0[ntwav]))
plt.xlabel("t")
plt.title("Wavelets")
plt.axis("tight")

# operator
PPop = pylops.avo.poststack.PoststackLinearModelling(wavs / 2, nt0=nt0, explicit=True)

# data
d = PPop * m

# solve
minv = pylops.avo.poststack.PoststackInversion(
    d, wavs / 2, m0=mback, explicit=True, **dict(cond=1e-10)
)[0]

fig, axs = plt.subplots(1, 2, figsize=(6, 7), sharey=True)
axs[0].plot(d, t0, "k", lw=4)
axs[0].set_title("Data")
axs[0].invert_yaxis()
axs[0].axis("tight")
axs[1].plot(m, t0, "k", lw=4, label="True")
axs[1].plot(mback, t0, "--b", lw=4, label="Back")
axs[1].plot(minv, t0, "--r", lw=2, label="Inv")
axs[1].set_title("Model")
axs[1].axis("tight")
axs[1].legend(loc=1)
plt.tight_layout()
  • Wavelets
  • Data, Model

We move now to a 2d example. First of all the model is loaded and data generated.

# model
inputfile = "../testdata/avo/poststack_model.npz"

model = np.load(inputfile)
m = np.log(model["model"][:, ::3])
x, z = model["x"][::3] / 1000.0, model["z"] / 1000.0
nx, nz = len(x), len(z)

# smooth model
nsmoothz, nsmoothx = 60, 50
mback = filtfilt(np.ones(nsmoothz) / float(nsmoothz), 1, m, axis=0)
mback = filtfilt(np.ones(nsmoothx) / float(nsmoothx), 1, mback, axis=1)

# dense operator
PPop_dense = pylops.avo.poststack.PoststackLinearModelling(
    wav / 2, nt0=nz, spatdims=nx, explicit=True
)

# lop operator
PPop = pylops.avo.poststack.PoststackLinearModelling(wav / 2, nt0=nz, spatdims=nx)

# data
d = (PPop_dense * m.ravel()).reshape(nz, nx)
n = np.random.normal(0, 1e-1, d.shape)
dn = d + n

Finally we perform 4 different inversions:

  • trace-by-trace inversion with explicit solver and dense operator with noise-free data

  • trace-by-trace inversion with explicit solver and dense operator with noisy data

  • multi-trace regularized inversion with iterative solver and linear operator using the result of trace-by-trace inversion as starting guess

    \[J = ||\Delta \mathbf{d} - \mathbf{W} \Delta \mathbf{ai}||_2 + \epsilon_\nabla ^2 ||\nabla \mathbf{ai}||_2\]

    where \(\Delta \mathbf{d}=\mathbf{d}-\mathbf{W}\mathbf{AI_0}\) is the residual data

  • multi-trace blocky inversion with iterative solver and linear operator

# dense inversion with noise-free data
minv_dense = pylops.avo.poststack.PoststackInversion(
    d, wav / 2, m0=mback, explicit=True, simultaneous=False
)[0]

# dense inversion with noisy data
minv_dense_noisy = pylops.avo.poststack.PoststackInversion(
    dn, wav / 2, m0=mback, explicit=True, epsI=4e-2, simultaneous=False
)[0]

# spatially regularized lop inversion with noisy data
minv_lop_reg = pylops.avo.poststack.PoststackInversion(
    dn,
    wav / 2,
    m0=minv_dense_noisy,
    explicit=False,
    epsR=5e1,
    **dict(damp=np.sqrt(1e-4), iter_lim=80)
)[0]

# blockiness promoting inversion with noisy data
minv_lop_blocky = pylops.avo.poststack.PoststackInversion(
    dn,
    wav / 2,
    m0=mback,
    explicit=False,
    epsR=[0.4],
    epsRL1=[0.1],
    **dict(mu=0.1, niter_outer=5, niter_inner=10, iter_lim=5, damp=1e-3)
)[0]

fig, axs = plt.subplots(2, 4, figsize=(15, 9))
axs[0][0].imshow(d, cmap="gray", extent=(x[0], x[-1], z[-1], z[0]), vmin=-0.4, vmax=0.4)
axs[0][0].set_title("Data")
axs[0][0].axis("tight")
axs[0][1].imshow(
    dn, cmap="gray", extent=(x[0], x[-1], z[-1], z[0]), vmin=-0.4, vmax=0.4
)
axs[0][1].set_title("Noisy Data")
axs[0][1].axis("tight")
axs[0][2].imshow(
    m,
    cmap="gist_rainbow",
    extent=(x[0], x[-1], z[-1], z[0]),
    vmin=m.min(),
    vmax=m.max(),
)
axs[0][2].set_title("Model")
axs[0][2].axis("tight")
axs[0][3].imshow(
    mback,
    cmap="gist_rainbow",
    extent=(x[0], x[-1], z[-1], z[0]),
    vmin=m.min(),
    vmax=m.max(),
)
axs[0][3].set_title("Smooth Model")
axs[0][3].axis("tight")
axs[1][0].imshow(
    minv_dense,
    cmap="gist_rainbow",
    extent=(x[0], x[-1], z[-1], z[0]),
    vmin=m.min(),
    vmax=m.max(),
)
axs[1][0].set_title("Noise-free Inversion")
axs[1][0].axis("tight")
axs[1][1].imshow(
    minv_dense_noisy,
    cmap="gist_rainbow",
    extent=(x[0], x[-1], z[-1], z[0]),
    vmin=m.min(),
    vmax=m.max(),
)
axs[1][1].set_title("Trace-by-trace Noisy Inversion")
axs[1][1].axis("tight")
axs[1][2].imshow(
    minv_lop_reg,
    cmap="gist_rainbow",
    extent=(x[0], x[-1], z[-1], z[0]),
    vmin=m.min(),
    vmax=m.max(),
)
axs[1][2].set_title("Regularized Noisy Inversion - lop ")
axs[1][2].axis("tight")
axs[1][3].imshow(
    minv_lop_blocky,
    cmap="gist_rainbow",
    extent=(x[0], x[-1], z[-1], z[0]),
    vmin=m.min(),
    vmax=m.max(),
)
axs[1][3].set_title("Blocky Noisy Inversion - lop ")
axs[1][3].axis("tight")

fig, ax = plt.subplots(1, 1, figsize=(3, 7))
ax.plot(m[:, nx // 2], z, "k", lw=4, label="True")
ax.plot(mback[:, nx // 2], z, "--r", lw=4, label="Back")
ax.plot(minv_dense[:, nx // 2], z, "--b", lw=2, label="Inv Dense")
ax.plot(minv_dense_noisy[:, nx // 2], z, "--m", lw=2, label="Inv Dense noisy")
ax.plot(minv_lop_reg[:, nx // 2], z, "--g", lw=2, label="Inv Lop regularized")
ax.plot(minv_lop_blocky[:, nx // 2], z, "--y", lw=2, label="Inv Lop blocky")
ax.set_title("Model")
ax.invert_yaxis()
ax.axis("tight")
ax.legend()
plt.tight_layout()
  • Data, Noisy Data, Model, Smooth Model, Noise-free Inversion, Trace-by-trace Noisy Inversion, Regularized Noisy Inversion - lop , Blocky Noisy Inversion - lop
  • Model

That’s almost it. If you wonder how this can be applied to real data, head over to the following notebook where the open-source segyio library is used alongside pylops to create an end-to-end open-source seismic inversion workflow with SEG-Y input data.

Total running time of the script: ( 0 minutes 12.010 seconds)

Gallery generated by Sphinx-Gallery

08. Pre-stack (AVO) inversion

Pre-stack inversion represents one step beyond post-stack inversion in that not only the profile of acoustic impedance can be inferred from seismic data, rather a set of elastic parameters is estimated from pre-stack data (i.e., angle gathers) using the information contained in the so-called AVO (amplitude versus offset) response. Such elastic parameters represent vital information for more sophisticated geophysical subsurface characterization than it would be possible to achieve working with post-stack seismic data.

In this tutorial, the pylops.avo.prestack.PrestackLinearModelling operator is used for modelling of both 1d and 2d synthetic pre-stack seismic data using 1d profiles or 2d models of different subsurface elastic parameters (P-wave velocity, S-wave velocity, and density) as input.

\[d(t, \theta) = w(t) * \sum_{i=1}^N G_i(t, \theta) \frac{\mathrm{d}\ln m_i(t)}{\mathrm{d}t}\]

where \(\mathbf{m}(t)=[V_P(t), V_S(t), \rho(t)]\) is a vector containing three elastic parameters at time \(t\), \(G_i(t, \theta)\) are the coefficients of the AVO parametrization used to model pre-stack data and \(w(t)\) is the time domain seismic wavelet. In compact form:

\[\mathbf{d}= \mathbf{W} \mathbf{G} \mathbf{D} \mathbf{m}\]

where \(\mathbf{W}\) is a convolution operator, \(\mathbf{G}\) is the AVO modelling operator, \(\mathbf{D}\) is a block-diagonal derivative operator, and \(\mathbf{m}\) is the input model. Subsequently the elastic parameters are estimated via the pylops.avo.prestack.PrestackInversion module. Once again, a two-steps inversion strategy can also be used to deal with the case of noisy data.

import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import filtfilt

import pylops
from pylops.utils.wavelets import ricker

plt.close("all")
np.random.seed(0)

Let’s start with a 1d example. A synthetic profile of acoustic impedance is created and data is modelled using both the dense and linear operator version of pylops.avo.prestack.PrestackLinearModelling operator

# sphinx_gallery_thumbnail_number = 5

# model
nt0 = 301
dt0 = 0.004

t0 = np.arange(nt0) * dt0
vp = 1200 + np.arange(nt0) + filtfilt(np.ones(5) / 5.0, 1, np.random.normal(0, 80, nt0))
vs = 600 + vp / 2 + filtfilt(np.ones(5) / 5.0, 1, np.random.normal(0, 20, nt0))
rho = 1000 + vp + filtfilt(np.ones(5) / 5.0, 1, np.random.normal(0, 30, nt0))
vp[131:] += 500
vs[131:] += 200
rho[131:] += 100
vsvp = 0.5
m = np.stack((np.log(vp), np.log(vs), np.log(rho)), axis=1)

# background model
nsmooth = 50
mback = filtfilt(np.ones(nsmooth) / float(nsmooth), 1, m, axis=0)

# angles
ntheta = 21
thetamin, thetamax = 0, 40
theta = np.linspace(thetamin, thetamax, ntheta)

# wavelet
ntwav = 41
wav = ricker(t0[: ntwav // 2 + 1], 15)[0]

# lop
PPop = pylops.avo.prestack.PrestackLinearModelling(
    wav, theta, vsvp=vsvp, nt0=nt0, linearization="akirich"
)

# dense
PPop_dense = pylops.avo.prestack.PrestackLinearModelling(
    wav, theta, vsvp=vsvp, nt0=nt0, linearization="akirich", explicit=True
)

# data lop
dPP = PPop * m.ravel()
dPP = dPP.reshape(nt0, ntheta)

# data dense
dPP_dense = PPop_dense * m.T.ravel()
dPP_dense = dPP_dense.reshape(ntheta, nt0).T

# noisy data
dPPn_dense = dPP_dense + np.random.normal(0, 1e-2, dPP_dense.shape)

We can now invert our data and retrieve elastic profiles for both noise-free and noisy data using pylops.avo.prestack.PrestackInversion.

# dense
minv_dense, dPP_dense_res = pylops.avo.prestack.PrestackInversion(
    dPP_dense,
    theta,
    wav,
    m0=mback,
    linearization="akirich",
    explicit=True,
    returnres=True,
    **dict(cond=1e-10)
)

# lop
minv, dPP_res = pylops.avo.prestack.PrestackInversion(
    dPP,
    theta,
    wav,
    m0=mback,
    linearization="akirich",
    explicit=False,
    returnres=True,
    **dict(damp=1e-10, iter_lim=2000)
)

# dense noisy
minv_dense_noise, dPPn_dense_res = pylops.avo.prestack.PrestackInversion(
    dPPn_dense,
    theta,
    wav,
    m0=mback,
    linearization="akirich",
    explicit=True,
    returnres=True,
    **dict(cond=1e-1)
)

# lop noisy (with vertical smoothing)
minv_noise, dPPn_res = pylops.avo.prestack.PrestackInversion(
    dPPn_dense,
    theta,
    wav,
    m0=mback,
    linearization="akirich",
    explicit=False,
    epsR=5e-1,
    returnres=True,
    **dict(damp=1e-1, iter_lim=100)
)

The data, inverted models and residuals are now displayed.

# data and model
fig, (axd, axdn, axvp, axvs, axrho) = plt.subplots(1, 5, figsize=(8, 5), sharey=True)
axd.imshow(
    dPP_dense,
    cmap="gray",
    extent=(theta[0], theta[-1], t0[-1], t0[0]),
    vmin=-np.abs(dPP_dense).max(),
    vmax=np.abs(dPP_dense).max(),
)
axd.set_title("Data")
axd.axis("tight")
axdn.imshow(
    dPPn_dense,
    cmap="gray",
    extent=(theta[0], theta[-1], t0[-1], t0[0]),
    vmin=-np.abs(dPP_dense).max(),
    vmax=np.abs(dPP_dense).max(),
)
axdn.set_title("Noisy Data")
axdn.axis("tight")
axvp.plot(vp, t0, "k", lw=4, label="True")
axvp.plot(np.exp(mback[:, 0]), t0, "--r", lw=4, label="Back")
axvp.plot(np.exp(minv_dense[:, 0]), t0, "--b", lw=2, label="Inv Dense")
axvp.plot(np.exp(minv[:, 0]), t0, "--m", lw=2, label="Inv Lop")
axvp.plot(np.exp(minv_dense_noise[:, 0]), t0, "--c", lw=2, label="Noisy Dense")
axvp.plot(np.exp(minv_noise[:, 0]), t0, "--g", lw=2, label="Noisy Lop")
axvp.set_title(r"$V_P$")
axvs.plot(vs, t0, "k", lw=4, label="True")
axvs.plot(np.exp(mback[:, 1]), t0, "--r", lw=4, label="Back")
axvs.plot(np.exp(minv_dense[:, 1]), t0, "--b", lw=2, label="Inv Dense")
axvs.plot(np.exp(minv[:, 1]), t0, "--m", lw=2, label="Inv Lop")
axvs.plot(np.exp(minv_dense_noise[:, 1]), t0, "--c", lw=2, label="Noisy Dense")
axvs.plot(np.exp(minv_noise[:, 1]), t0, "--g", lw=2, label="Noisy Lop")
axvs.set_title(r"$V_S$")
axrho.plot(rho, t0, "k", lw=4, label="True")
axrho.plot(np.exp(mback[:, 2]), t0, "--r", lw=4, label="Back")
axrho.plot(np.exp(minv_dense[:, 2]), t0, "--b", lw=2, label="Inv Dense")
axrho.plot(np.exp(minv[:, 2]), t0, "--m", lw=2, label="Inv Lop")
axrho.plot(np.exp(minv_dense_noise[:, 2]), t0, "--c", lw=2, label="Noisy Dense")
axrho.plot(np.exp(minv_noise[:, 2]), t0, "--g", lw=2, label="Noisy Lop")
axrho.set_title(r"$\rho$")
axrho.legend(loc="center left", bbox_to_anchor=(1, 0.5))
axd.axis("tight")
plt.tight_layout()

# residuals
fig, axs = plt.subplots(1, 4, figsize=(8, 5), sharey=True)
fig.suptitle("Residuals", fontsize=14, fontweight="bold", y=0.95)
im = axs[0].imshow(
    dPP_dense_res,
    cmap="gray",
    extent=(theta[0], theta[-1], t0[-1], t0[0]),
    vmin=-0.1,
    vmax=0.1,
)
axs[0].set_title("Dense")
axs[0].set_xlabel(r"$\theta$")
axs[0].set_ylabel("t[s]")
axs[0].axis("tight")
axs[1].imshow(
    dPP_res,
    cmap="gray",
    extent=(theta[0], theta[-1], t0[-1], t0[0]),
    vmin=-0.1,
    vmax=0.1,
)
axs[1].set_title("Lop")
axs[1].set_xlabel(r"$\theta$")
axs[1].axis("tight")
axs[2].imshow(
    dPPn_dense_res,
    cmap="gray",
    extent=(theta[0], theta[-1], t0[-1], t0[0]),
    vmin=-0.1,
    vmax=0.1,
)
axs[2].set_title("Noisy Dense")
axs[2].set_xlabel(r"$\theta$")
axs[2].axis("tight")
axs[3].imshow(
    dPPn_res,
    cmap="gray",
    extent=(theta[0], theta[-1], t0[-1], t0[0]),
    vmin=-0.1,
    vmax=0.1,
)
axs[3].set_title("Noisy Lop")
axs[3].set_xlabel(r"$\theta$")
axs[3].axis("tight")
plt.tight_layout()
plt.subplots_adjust(top=0.85)
  • Data, Noisy Data, $V_P$, $V_S$, $\rho$
  • Residuals, Dense, Lop, Noisy Dense, Noisy Lop

Finally before moving to the 2d example, we consider the case when both PP and PS data are available. A joint PP-PS inversion can be easily solved as follows.

PSop = pylops.avo.prestack.PrestackLinearModelling(
    2 * wav, theta, vsvp=vsvp, nt0=nt0, linearization="ps"
)
PPPSop = pylops.VStack((PPop, PSop))

# data
dPPPS = PPPSop * m.ravel()
dPPPS = dPPPS.reshape(2, nt0, ntheta)

dPPPSn = dPPPS + np.random.normal(0, 1e-2, dPPPS.shape)

# Invert
minvPPSP, dPPPS_res = pylops.avo.prestack.PrestackInversion(
    dPPPS,
    theta,
    [wav, 2 * wav],
    m0=mback,
    linearization=["fatti", "ps"],
    epsR=5e-1,
    returnres=True,
    **dict(damp=1e-1, iter_lim=100)
)

# Data and model
fig, (axd, axdn, axvp, axvs, axrho) = plt.subplots(1, 5, figsize=(8, 5), sharey=True)
axd.imshow(
    dPPPSn[0],
    cmap="gray",
    extent=(theta[0], theta[-1], t0[-1], t0[0]),
    vmin=-np.abs(dPPPSn[0]).max(),
    vmax=np.abs(dPPPSn[0]).max(),
)
axd.set_title("PP Data")
axd.axis("tight")
axdn.imshow(
    dPPPSn[1],
    cmap="gray",
    extent=(theta[0], theta[-1], t0[-1], t0[0]),
    vmin=-np.abs(dPPPSn[1]).max(),
    vmax=np.abs(dPPPSn[1]).max(),
)
axdn.set_title("PS Data")
axdn.axis("tight")
axvp.plot(vp, t0, "k", lw=4, label="True")
axvp.plot(np.exp(mback[:, 0]), t0, "--r", lw=4, label="Back")
axvp.plot(np.exp(minv_noise[:, 0]), t0, "--g", lw=2, label="PP")
axvp.plot(np.exp(minvPPSP[:, 0]), t0, "--b", lw=2, label="PP+PS")
axvp.set_title(r"$V_P$")
axvs.plot(vs, t0, "k", lw=4, label="True")
axvs.plot(np.exp(mback[:, 1]), t0, "--r", lw=4, label="Back")
axvs.plot(np.exp(minv_noise[:, 1]), t0, "--g", lw=2, label="PP")
axvs.plot(np.exp(minvPPSP[:, 1]), t0, "--b", lw=2, label="PP+PS")
axvs.set_title(r"$V_S$")
axrho.plot(rho, t0, "k", lw=4, label="True")
axrho.plot(np.exp(mback[:, 2]), t0, "--r", lw=4, label="Back")
axrho.plot(np.exp(minv_noise[:, 2]), t0, "--g", lw=2, label="PP")
axrho.plot(np.exp(minvPPSP[:, 2]), t0, "--b", lw=2, label="PP+PS")
axrho.set_title(r"$\rho$")
axrho.legend(loc="center left", bbox_to_anchor=(1, 0.5))
axd.axis("tight")
plt.tight_layout()
PP Data, PS Data, $V_P$, $V_S$, $\rho$

We move now to a 2d example. First of all the model is loaded and data generated.

# model
inputfile = "../testdata/avo/poststack_model.npz"

model = np.load(inputfile)
x, z = model["x"][::6] / 1000.0, model["z"][:300] / 1000.0
nx, nz = len(x), len(z)
m = 1000 * model["model"][:300, ::6]

mvp = m.copy()
mvs = m / 2
mrho = m / 3 + 300
m = np.log(np.stack((mvp, mvs, mrho), axis=1))

# smooth model
nsmoothz, nsmoothx = 30, 25
mback = filtfilt(np.ones(nsmoothz) / float(nsmoothz), 1, m, axis=0)
mback = filtfilt(np.ones(nsmoothx) / float(nsmoothx), 1, mback, axis=2)

# dense operator
PPop_dense = pylops.avo.prestack.PrestackLinearModelling(
    wav,
    theta,
    vsvp=vsvp,
    nt0=nz,
    spatdims=(nx,),
    linearization="akirich",
    explicit=True,
)

# lop operator
PPop = pylops.avo.prestack.PrestackLinearModelling(
    wav, theta, vsvp=vsvp, nt0=nz, spatdims=(nx,), linearization="akirich"
)

# data
dPP = PPop_dense * m.swapaxes(0, 1).ravel()
dPP = dPP.reshape(ntheta, nz, nx).swapaxes(0, 1)
dPPn = dPP + np.random.normal(0, 5e-2, dPP.shape)

Finally we perform the same 4 different inversions as in the post-stack tutorial (see 07. Post-stack inversion for more details).

# dense inversion with noise-free data
minv_dense = pylops.avo.prestack.PrestackInversion(
    dPP, theta, wav, m0=mback, explicit=True, simultaneous=False
)

# dense inversion with noisy data
minv_dense_noisy = pylops.avo.prestack.PrestackInversion(
    dPPn, theta, wav, m0=mback, explicit=True, epsI=4e-2, simultaneous=False
)

# spatially regularized lop inversion with noisy data
minv_lop_reg = pylops.avo.prestack.PrestackInversion(
    dPPn,
    theta,
    wav,
    m0=minv_dense_noisy,
    explicit=False,
    epsR=1e1,
    **dict(damp=np.sqrt(1e-4), iter_lim=20)
)

# blockiness promoting inversion with noisy data
minv_blocky = pylops.avo.prestack.PrestackInversion(
    dPPn,
    theta,
    wav,
    m0=mback,
    explicit=False,
    epsR=0.4,
    epsRL1=0.1,
    **dict(mu=0.1, niter_outer=3, niter_inner=3, iter_lim=5, damp=1e-3)
)

Let’s now visualize the inverted elastic parameters for the different scenarios

def plotmodel(
    axs,
    m,
    x,
    z,
    vmin,
    vmax,
    params=("VP", "VS", "Rho"),
    cmap="gist_rainbow",
    title=None,
):
    """Quick visualization of model"""
    for ip, param in enumerate(params):
        axs[ip].imshow(
            m[:, ip], extent=(x[0], x[-1], z[-1], z[0]), vmin=vmin, vmax=vmax, cmap=cmap
        )
        axs[ip].set_title("%s - %s" % (param, title))
        axs[ip].axis("tight")
    plt.setp(axs[1].get_yticklabels(), visible=False)
    plt.setp(axs[2].get_yticklabels(), visible=False)


# data
fig = plt.figure(figsize=(8, 9))
ax1 = plt.subplot2grid((2, 3), (0, 0), colspan=3)
ax2 = plt.subplot2grid((2, 3), (1, 0))
ax3 = plt.subplot2grid((2, 3), (1, 1), sharey=ax2)
ax4 = plt.subplot2grid((2, 3), (1, 2), sharey=ax2)
ax1.imshow(
    dPP[:, 0], cmap="gray", extent=(x[0], x[-1], z[-1], z[0]), vmin=-0.4, vmax=0.4
)
ax1.vlines(
    [x[nx // 5], x[nx // 2], x[4 * nx // 5]],
    ymin=z[0],
    ymax=z[-1],
    colors="w",
    linestyles="--",
)
ax1.set_xlabel("x [km]")
ax1.set_ylabel("z [km]")
ax1.set_title(r"Stack ($\theta$=0)")
ax1.axis("tight")
ax2.imshow(
    dPP[:, :, nx // 5],
    cmap="gray",
    extent=(theta[0], theta[-1], z[-1], z[0]),
    vmin=-0.4,
    vmax=0.4,
)
ax2.set_xlabel(r"$\theta$")
ax2.set_ylabel("z [km]")
ax2.set_title(r"Gather (x=%.2f)" % x[nx // 5])
ax2.axis("tight")
ax3.imshow(
    dPP[:, :, nx // 2],
    cmap="gray",
    extent=(theta[0], theta[-1], z[-1], z[0]),
    vmin=-0.4,
    vmax=0.4,
)
ax3.set_xlabel(r"$\theta$")
ax3.set_title(r"Gather (x=%.2f)" % x[nx // 2])
ax3.axis("tight")
ax4.imshow(
    dPP[:, :, 4 * nx // 5],
    cmap="gray",
    extent=(theta[0], theta[-1], z[-1], z[0]),
    vmin=-0.4,
    vmax=0.4,
)
ax4.set_xlabel(r"$\theta$")
ax4.set_title(r"Gather (x=%.2f)" % x[4 * nx // 5])
ax4.axis("tight")
plt.setp(ax3.get_yticklabels(), visible=False)
plt.setp(ax4.get_yticklabels(), visible=False)

# noisy data
fig = plt.figure(figsize=(8, 9))
ax1 = plt.subplot2grid((2, 3), (0, 0), colspan=3)
ax2 = plt.subplot2grid((2, 3), (1, 0))
ax3 = plt.subplot2grid((2, 3), (1, 1), sharey=ax2)
ax4 = plt.subplot2grid((2, 3), (1, 2), sharey=ax2)
ax1.imshow(
    dPPn[:, 0], cmap="gray", extent=(x[0], x[-1], z[-1], z[0]), vmin=-0.4, vmax=0.4
)
ax1.vlines(
    [x[nx // 5], x[nx // 2], x[4 * nx // 5]],
    ymin=z[0],
    ymax=z[-1],
    colors="w",
    linestyles="--",
)
ax1.set_xlabel("x [km]")
ax1.set_ylabel("z [km]")
ax1.set_title(r"Noisy Stack ($\theta$=0)")
ax1.axis("tight")
ax2.imshow(
    dPPn[:, :, nx // 5],
    cmap="gray",
    extent=(theta[0], theta[-1], z[-1], z[0]),
    vmin=-0.4,
    vmax=0.4,
)
ax2.set_xlabel(r"$\theta$")
ax2.set_ylabel("z [km]")
ax2.set_title(r"Gather (x=%.2f)" % x[nx // 5])
ax2.axis("tight")
ax3.imshow(
    dPPn[:, :, nx // 2],
    cmap="gray",
    extent=(theta[0], theta[-1], z[-1], z[0]),
    vmin=-0.4,
    vmax=0.4,
)
ax3.set_title(r"Gather (x=%.2f)" % x[nx // 2])
ax3.set_xlabel(r"$\theta$")
ax3.axis("tight")
ax4.imshow(
    dPPn[:, :, 4 * nx // 5],
    cmap="gray",
    extent=(theta[0], theta[-1], z[-1], z[0]),
    vmin=-0.4,
    vmax=0.4,
)
ax4.set_xlabel(r"$\theta$")
ax4.set_title(r"Gather (x=%.2f)" % x[4 * nx // 5])
ax4.axis("tight")
plt.setp(ax3.get_yticklabels(), visible=False)
plt.setp(ax4.get_yticklabels(), visible=False)

# inverted models
fig, axs = plt.subplots(6, 3, figsize=(8, 19))
fig.suptitle("Model", fontsize=12, fontweight="bold", y=0.95)
plotmodel(axs[0], m, x, z, m.min(), m.max(), title="True")
plotmodel(axs[1], mback, x, z, m.min(), m.max(), title="Back")
plotmodel(axs[2], minv_dense, x, z, m.min(), m.max(), title="Dense")
plotmodel(axs[3], minv_dense_noisy, x, z, m.min(), m.max(), title="Dense noisy")
plotmodel(axs[4], minv_lop_reg, x, z, m.min(), m.max(), title="Lop regularized")
plotmodel(axs[5], minv_blocky, x, z, m.min(), m.max(), title="Lop blocky")
plt.tight_layout()
plt.subplots_adjust(top=0.92)

fig, axs = plt.subplots(1, 3, figsize=(8, 7))
for ip, param in enumerate(["VP", "VS", "Rho"]):
    axs[ip].plot(m[:, ip, nx // 2], z, "k", lw=4, label="True")
    axs[ip].plot(mback[:, ip, nx // 2], z, "--r", lw=4, label="Back")
    axs[ip].plot(minv_dense[:, ip, nx // 2], z, "--b", lw=2, label="Inv Dense")
    axs[ip].plot(
        minv_dense_noisy[:, ip, nx // 2], z, "--m", lw=2, label="Inv Dense noisy"
    )
    axs[ip].plot(
        minv_lop_reg[:, ip, nx // 2], z, "--g", lw=2, label="Inv Lop regularized"
    )
    axs[ip].plot(minv_blocky[:, ip, nx // 2], z, "--y", lw=2, label="Inv Lop blocky")
    axs[ip].set_title(param)
    axs[ip].invert_yaxis()
axs[2].legend(loc=8, fontsize="small")
plt.tight_layout()
  • Stack ($\theta$=0), Gather (x=1.42), Gather (x=2.41), Gather (x=3.37)
  • Noisy Stack ($\theta$=0), Gather (x=1.42), Gather (x=2.41), Gather (x=3.37)
  • Model, VP - True, VS - True, Rho - True, VP - Back, VS - Back, Rho - Back, VP - Dense, VS - Dense, Rho - Dense, VP - Dense noisy, VS - Dense noisy, Rho - Dense noisy, VP - Lop regularized, VS - Lop regularized, Rho - Lop regularized, VP - Lop blocky, VS - Lop blocky, Rho - Lop blocky
  • VP, VS, Rho

While the background model m0 has been provided in all the examples so far, it is worth showing that the module pylops.avo.prestack.PrestackInversion can also produce so-called relative elastic parameters (i.e., variations from an average medium property) when the background model m0 is not available.

dminv = pylops.avo.prestack.PrestackInversion(
    dPP, theta, wav, m0=None, explicit=True, simultaneous=False
)

fig, axs = plt.subplots(1, 3, figsize=(8, 3))
plotmodel(axs, dminv, x, z, -dminv.max(), dminv.max(), cmap="seismic", title="relative")

fig, axs = plt.subplots(1, 3, figsize=(8, 7))
for ip, param in enumerate(["VP", "VS", "Rho"]):
    axs[ip].plot(dminv[:, ip, nx // 2], z, "k", lw=2)
    axs[ip].set_title(param)
    axs[ip].invert_yaxis()
  • VP - relative, VS - relative, Rho - relative
  • VP, VS, Rho

Total running time of the script: ( 0 minutes 31.440 seconds)

Gallery generated by Sphinx-Gallery

09. Multi-Dimensional Deconvolution

This example shows how to set-up and run the pylops.waveeqprocessing.MDD inversion using synthetic data.

import warnings

import matplotlib.pyplot as plt
import numpy as np

import pylops
from pylops.utils.seismicevents import hyperbolic2d, makeaxis
from pylops.utils.tapers import taper3d
from pylops.utils.wavelets import ricker

warnings.filterwarnings("ignore")
plt.close("all")

# sphinx_gallery_thumbnail_number = 5

Let’s start by creating a set of hyperbolic events to be used as our MDC kernel

# Input parameters
par = {
    "ox": -150,
    "dx": 10,
    "nx": 31,
    "oy": -250,
    "dy": 10,
    "ny": 51,
    "ot": 0,
    "dt": 0.004,
    "nt": 300,
    "f0": 20,
    "nfmax": 200,
}

t0_m = [0.2]
vrms_m = [700.0]
amp_m = [1.0]

t0_G = [0.2, 0.5, 0.7]
vrms_G = [800.0, 1200.0, 1500.0]
amp_G = [1.0, 0.6, 0.5]

# Taper
tap = taper3d(par["nt"], [par["ny"], par["nx"]], (5, 5), tapertype="hanning")

# Create axis
t, t2, x, y = makeaxis(par)

# Create wavelet
wav = ricker(t[:41], f0=par["f0"])[0]

# Generate model
m, mwav = hyperbolic2d(x, t, t0_m, vrms_m, amp_m, wav)

# Generate operator
G, Gwav = np.zeros((par["ny"], par["nx"], par["nt"])), np.zeros(
    (par["ny"], par["nx"], par["nt"])
)
for iy, y0 in enumerate(y):
    G[iy], Gwav[iy] = hyperbolic2d(x - y0, t, t0_G, vrms_G, amp_G, wav)
G, Gwav = G * tap, Gwav * tap

# Add negative part to data and model
m = np.concatenate((np.zeros((par["nx"], par["nt"] - 1)), m), axis=-1)
mwav = np.concatenate((np.zeros((par["nx"], par["nt"] - 1)), mwav), axis=-1)
Gwav2 = np.concatenate((np.zeros((par["ny"], par["nx"], par["nt"] - 1)), Gwav), axis=-1)

# Define MDC linear operator
Gwav_fft = np.fft.rfft(Gwav2, 2 * par["nt"] - 1, axis=-1)
Gwav_fft = (Gwav_fft[..., : par["nfmax"]]).transpose(2, 0, 1)

MDCop = pylops.waveeqprocessing.MDC(
    Gwav_fft,
    nt=2 * par["nt"] - 1,
    nv=1,
    dt=0.004,
    dr=1.0,
)

# Create data
d = MDCop * m.T.ravel()
d = d.reshape(2 * par["nt"] - 1, par["ny"]).T

Let’s display what we have so far: operator, input model, and data

fig, axs = plt.subplots(1, 2, figsize=(8, 6))
axs[0].imshow(
    Gwav2[int(par["ny"] / 2)].T,
    aspect="auto",
    interpolation="nearest",
    cmap="gray",
    vmin=-np.abs(Gwav2.max()),
    vmax=np.abs(Gwav2.max()),
    extent=(x.min(), x.max(), t2.max(), t2.min()),
)
axs[0].set_title("G - inline view", fontsize=15)
axs[0].set_xlabel(r"$x_R$")
axs[1].set_ylabel(r"$t$")
axs[1].imshow(
    Gwav2[:, int(par["nx"] / 2)].T,
    aspect="auto",
    interpolation="nearest",
    cmap="gray",
    vmin=-np.abs(Gwav2.max()),
    vmax=np.abs(Gwav2.max()),
    extent=(y.min(), y.max(), t2.max(), t2.min()),
)
axs[1].set_title("G - inline view", fontsize=15)
axs[1].set_xlabel(r"$x_S$")
axs[1].set_ylabel(r"$t$")
fig.tight_layout()

fig, axs = plt.subplots(1, 2, figsize=(8, 6))
axs[0].imshow(
    mwav.T,
    aspect="auto",
    interpolation="nearest",
    cmap="gray",
    vmin=-np.abs(mwav.max()),
    vmax=np.abs(mwav.max()),
    extent=(x.min(), x.max(), t2.max(), t2.min()),
)
axs[0].set_title(r"$m$", fontsize=15)
axs[0].set_xlabel(r"$x_R$")
axs[1].set_ylabel(r"$t$")
axs[1].imshow(
    d.T,
    aspect="auto",
    interpolation="nearest",
    cmap="gray",
    vmin=-np.abs(d.max()),
    vmax=np.abs(d.max()),
    extent=(x.min(), x.max(), t2.max(), t2.min()),
)
axs[1].set_title(r"$d$", fontsize=15)
axs[1].set_xlabel(r"$x_S$")
axs[1].set_ylabel(r"$t$")
fig.tight_layout()
  • G - inline view, G - inline view
  • $m$, $d$

We are now ready to feed our operator to pylops.waveeqprocessing.MDD and invert back for our input model

minv, madj, psfinv, psfadj = pylops.waveeqprocessing.MDD(
    Gwav,
    d[:, par["nt"] - 1 :],
    dt=par["dt"],
    dr=par["dx"],
    nfmax=par["nfmax"],
    wav=wav,
    twosided=True,
    add_negative=True,
    adjoint=True,
    psf=True,
    dottest=False,
    **dict(damp=1e-4, iter_lim=20, show=0)
)

fig = plt.figure(figsize=(8, 6))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))
ax1.imshow(
    madj.T,
    aspect="auto",
    interpolation="nearest",
    cmap="gray",
    vmin=-np.abs(madj.max()),
    vmax=np.abs(madj.max()),
    extent=(x.min(), x.max(), t2.max(), t2.min()),
)
ax1.set_title("Adjoint m", fontsize=15)
ax1.set_xlabel(r"$x_V$")
axs[1].set_ylabel(r"$t$")
ax2.imshow(
    minv.T,
    aspect="auto",
    interpolation="nearest",
    cmap="gray",
    vmin=-np.abs(minv.max()),
    vmax=np.abs(minv.max()),
    extent=(x.min(), x.max(), t2.max(), t2.min()),
)
ax2.set_title("Inverted m", fontsize=15)
ax2.set_xlabel(r"$x_V$")
axs[1].set_ylabel(r"$t$")
ax3.plot(
    madj[int(par["nx"] / 2)] / np.abs(madj[int(par["nx"] / 2)]).max(), t2, "r", lw=5
)
ax3.plot(
    minv[int(par["nx"] / 2)] / np.abs(minv[int(par["nx"] / 2)]).max(), t2, "k", lw=3
)
ax3.set_ylim([t2[-1], t2[0]])
fig.tight_layout()

fig, axs = plt.subplots(1, 2, figsize=(8, 6))
axs[0].imshow(
    psfinv[int(par["nx"] / 2)].T,
    aspect="auto",
    interpolation="nearest",
    vmin=-np.abs(psfinv.max()),
    vmax=np.abs(psfinv.max()),
    cmap="gray",
    extent=(x.min(), x.max(), t2.max(), t2.min()),
)
axs[0].set_title("Inverted psf - inline view", fontsize=15)
axs[0].set_xlabel(r"$x_V$")
axs[1].set_ylabel(r"$t$")
axs[1].imshow(
    psfinv[:, int(par["nx"] / 2)].T,
    aspect="auto",
    interpolation="nearest",
    vmin=-np.abs(psfinv.max()),
    vmax=np.abs(psfinv.max()),
    cmap="gray",
    extent=(y.min(), y.max(), t2.max(), t2.min()),
)
axs[1].set_title("Inverted psf - xline view", fontsize=15)
axs[1].set_xlabel(r"$x_V$")
axs[1].set_ylabel(r"$t$")
fig.tight_layout()
  • Adjoint m, Inverted m
  • Inverted psf - inline view, Inverted psf - xline view

We repeat the same procedure but this time we will add a preconditioning by means of causality_precond parameter, which enforces the inverted model to be zero in the negative part of the time axis (as expected by theory). This preconditioning will have the effect of speeding up the convergence of the iterative solver and thus reduce the computation time of the deconvolution

minvprec = pylops.waveeqprocessing.MDD(
    Gwav,
    d[:, par["nt"] - 1 :],
    dt=par["dt"],
    dr=par["dx"],
    nfmax=par["nfmax"],
    wav=wav,
    twosided=True,
    add_negative=True,
    adjoint=False,
    psf=False,
    causality_precond=True,
    dottest=False,
    **dict(damp=1e-4, iter_lim=50, show=0)
)

fig = plt.figure(figsize=(8, 6))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))
ax1.imshow(
    madj.T,
    aspect="auto",
    interpolation="nearest",
    cmap="gray",
    vmin=-np.abs(madj.max()),
    vmax=np.abs(madj.max()),
    extent=(x.min(), x.max(), t2.max(), t2.min()),
)
ax1.set_title("Adjoint m", fontsize=15)
ax1.set_xlabel(r"$x_V$")
axs[1].set_ylabel(r"$t$")
ax2.imshow(
    minvprec.T,
    aspect="auto",
    interpolation="nearest",
    cmap="gray",
    vmin=-np.abs(minvprec.max()),
    vmax=np.abs(minvprec.max()),
    extent=(x.min(), x.max(), t2.max(), t2.min()),
)
ax2.set_title("Inverted m", fontsize=15)
ax2.set_xlabel(r"$x_V$")
axs[1].set_ylabel(r"$t$")
ax3.plot(
    madj[int(par["nx"] / 2)] / np.abs(madj[int(par["nx"] / 2)]).max(), t2, "r", lw=5
)
ax3.plot(
    minvprec[int(par["nx"] / 2)] / np.abs(minv[int(par["nx"] / 2)]).max(), t2, "k", lw=3
)
ax3.set_ylim([t2[-1], t2[0]])
fig.tight_layout()
Adjoint m, Inverted m

Total running time of the script: ( 0 minutes 11.100 seconds)

Gallery generated by Sphinx-Gallery

10. Marchenko redatuming by inversion

This example shows how to set-up and run the pylops.waveeqprocessing.Marchenko inversion using synthetic data.

# sphinx_gallery_thumbnail_number = 5
# pylint: disable=C0103
import warnings

import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import convolve

from pylops.waveeqprocessing import Marchenko

warnings.filterwarnings("ignore")
plt.close("all")

Let’s start by defining some input parameters and loading the test data

# Input parameters
inputfile = "../testdata/marchenko/input.npz"

vel = 2400.0  # velocity
toff = 0.045  # direct arrival time shift
nsmooth = 10  # time window smoothing
nfmax = 1000  # max frequency for MDC (#samples)
niter = 10  # iterations

inputdata = np.load(inputfile)

# Receivers
r = inputdata["r"]
nr = r.shape[1]
dr = r[0, 1] - r[0, 0]

# Sources
s = inputdata["s"]
ns = s.shape[1]
ds = s[0, 1] - s[0, 0]

# Virtual points
vs = inputdata["vs"]

# Density model
rho = inputdata["rho"]
z, x = inputdata["z"], inputdata["x"]

# Reflection data (R[s, r, t]) and subsurface fields
R = inputdata["R"][:, :, :-100]
R = np.swapaxes(R, 0, 1)  # just because of how the data was saved

Gsub = inputdata["Gsub"][:-100]
G0sub = inputdata["G0sub"][:-100]
wav = inputdata["wav"]
wav_c = np.argmax(wav)

t = inputdata["t"][:-100]
ot, dt, nt = t[0], t[1] - t[0], len(t)

Gsub = np.apply_along_axis(convolve, 0, Gsub, wav, mode="full")
Gsub = Gsub[wav_c:][:nt]
G0sub = np.apply_along_axis(convolve, 0, G0sub, wav, mode="full")
G0sub = G0sub[wav_c:][:nt]

plt.figure(figsize=(10, 5))
plt.imshow(rho, cmap="gray", extent=(x[0], x[-1], z[-1], z[0]))
plt.scatter(s[0, 5::10], s[1, 5::10], marker="*", s=150, c="r", edgecolors="k")
plt.scatter(r[0, ::10], r[1, ::10], marker="v", s=150, c="b", edgecolors="k")
plt.scatter(vs[0], vs[1], marker=".", s=250, c="m", edgecolors="k")
plt.axis("tight")
plt.xlabel("x [m]")
plt.ylabel("y [m]")
plt.title("Model and Geometry")
plt.xlim(x[0], x[-1])

fig, axs = plt.subplots(1, 3, sharey=True, figsize=(12, 7))
axs[0].imshow(
    R[0].T, cmap="gray", vmin=-1e-2, vmax=1e-2, extent=(r[0, 0], r[0, -1], t[-1], t[0])
)
axs[0].set_title("R shot=0")
axs[0].set_xlabel(r"$x_R$")
axs[0].set_ylabel(r"$t$")
axs[0].axis("tight")
axs[0].set_ylim(1.5, 0)
axs[1].imshow(
    R[ns // 2].T,
    cmap="gray",
    vmin=-1e-2,
    vmax=1e-2,
    extent=(r[0, 0], r[0, -1], t[-1], t[0]),
)
axs[1].set_title("R shot=%d" % (ns // 2))
axs[1].set_xlabel(r"$x_R$")
axs[1].set_ylabel(r"$t$")
axs[1].axis("tight")
axs[1].set_ylim(1.5, 0)
axs[2].imshow(
    R[-1].T, cmap="gray", vmin=-1e-2, vmax=1e-2, extent=(r[0, 0], r[0, -1], t[-1], t[0])
)
axs[2].set_title("R shot=%d" % ns)
axs[2].set_xlabel(r"$x_R$")
axs[2].axis("tight")
axs[2].set_ylim(1.5, 0)
fig.tight_layout()

fig, axs = plt.subplots(1, 2, sharey=True, figsize=(8, 6))
axs[0].imshow(
    Gsub, cmap="gray", vmin=-1e6, vmax=1e6, extent=(r[0, 0], r[0, -1], t[-1], t[0])
)
axs[0].set_title("G")
axs[0].set_xlabel(r"$x_R$")
axs[0].set_ylabel(r"$t$")
axs[0].axis("tight")
axs[0].set_ylim(1.5, 0)
axs[1].imshow(
    G0sub, cmap="gray", vmin=-1e6, vmax=1e6, extent=(r[0, 0], r[0, -1], t[-1], t[0])
)
axs[1].set_title("G0")
axs[1].set_xlabel(r"$x_R$")
axs[1].set_ylabel(r"$t$")
axs[1].axis("tight")
axs[1].set_ylim(1.5, 0)
fig.tight_layout()
  • Model and Geometry
  • R shot=0, R shot=50, R shot=101
  • G, G0

Let’s now create an object of the pylops.waveeqprocessing.Marchenko class and apply redatuming for a single subsurface point vs.

# direct arrival window
trav = np.sqrt((vs[0] - r[0]) ** 2 + (vs[1] - r[1]) ** 2) / vel

MarchenkoWM = Marchenko(
    R, dt=dt, dr=dr, nfmax=nfmax, wav=wav, toff=toff, nsmooth=nsmooth
)

(
    f1_inv_minus,
    f1_inv_plus,
    p0_minus,
    g_inv_minus,
    g_inv_plus,
) = MarchenkoWM.apply_onepoint(
    trav,
    G0=G0sub.T,
    rtm=True,
    greens=True,
    dottest=True,
    **dict(iter_lim=niter, show=True)
)
g_inv_tot = g_inv_minus + g_inv_plus
Dot test passed, v^H(Opu)=405.16509639614344 - u^H(Op^Hv)=405.1650963961446
Dot test passed, v^H(Opu)=172.06507561051328 - u^H(Op^Hv)=172.06507561051316

LSQR            Least-squares solution of  Ax = b
The matrix A has 282598 rows and 282598 columns
damp = 0.00000000000000e+00   calc_var =        0
atol = 1.00e-06                 conlim = 1.00e+08
btol = 1.00e-06               iter_lim =       10

   Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
     0  0.00000e+00   3.134e+07  3.134e+07    1.0e+00  3.3e-08
     1  0.00000e+00   1.374e+07  1.374e+07    4.4e-01  9.3e-01   1.1e+00  1.0e+00
     2  0.00000e+00   7.770e+06  7.770e+06    2.5e-01  3.9e-01   1.8e+00  2.2e+00
     3  0.00000e+00   5.750e+06  5.750e+06    1.8e-01  3.3e-01   2.1e+00  3.4e+00
     4  0.00000e+00   3.930e+06  3.930e+06    1.3e-01  3.4e-01   2.5e+00  5.1e+00
     5  0.00000e+00   3.042e+06  3.042e+06    9.7e-02  2.6e-01   2.9e+00  6.8e+00
     6  0.00000e+00   2.423e+06  2.423e+06    7.7e-02  2.2e-01   3.3e+00  8.6e+00
     7  0.00000e+00   1.675e+06  1.675e+06    5.3e-02  2.5e-01   3.6e+00  1.1e+01
     8  0.00000e+00   1.248e+06  1.248e+06    4.0e-02  2.0e-01   3.9e+00  1.3e+01
     9  0.00000e+00   1.004e+06  1.004e+06    3.2e-02  1.5e-01   4.2e+00  1.4e+01
    10  0.00000e+00   7.762e+05  7.762e+05    2.5e-02  1.8e-01   4.4e+00  1.6e+01

LSQR finished
The iteration limit has been reached

istop =       7   r1norm = 7.8e+05   anorm = 4.4e+00   arnorm = 6.1e+05
itn   =      10   r2norm = 7.8e+05   acond = 1.6e+01   xnorm  = 3.6e+07

We can now compare the result of Marchenko redatuming via LSQR with standard redatuming

fig, axs = plt.subplots(1, 3, sharey=True, figsize=(12, 7))
axs[0].imshow(
    p0_minus.T,
    cmap="gray",
    vmin=-5e5,
    vmax=5e5,
    extent=(r[0, 0], r[0, -1], t[-1], -t[-1]),
)
axs[0].set_title(r"$p_0^-$")
axs[0].set_xlabel(r"$x_R$")
axs[0].set_ylabel(r"$t$")
axs[0].axis("tight")
axs[0].set_ylim(1.2, 0)
axs[1].imshow(
    g_inv_minus.T,
    cmap="gray",
    vmin=-5e5,
    vmax=5e5,
    extent=(r[0, 0], r[0, -1], t[-1], -t[-1]),
)
axs[1].set_title(r"$g^-$")
axs[1].set_xlabel(r"$x_R$")
axs[1].set_ylabel(r"$t$")
axs[1].axis("tight")
axs[1].set_ylim(1.2, 0)
axs[2].imshow(
    g_inv_plus.T,
    cmap="gray",
    vmin=-5e5,
    vmax=5e5,
    extent=(r[0, 0], r[0, -1], t[-1], -t[-1]),
)
axs[2].set_title(r"$g^+$")
axs[2].set_xlabel(r"$x_R$")
axs[2].set_ylabel(r"$t$")
axs[2].axis("tight")
axs[2].set_ylim(1.2, 0)
fig.tight_layout()

fig = plt.figure(figsize=(12, 7))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))
ax1.imshow(
    Gsub, cmap="gray", vmin=-5e5, vmax=5e5, extent=(r[0, 0], r[0, -1], t[-1], t[0])
)
ax1.set_title(r"$G_{true}$")
axs[0].set_xlabel(r"$x_R$")
axs[0].set_ylabel(r"$t$")
ax1.axis("tight")
ax1.set_ylim(1.2, 0)
ax2.imshow(
    g_inv_tot.T,
    cmap="gray",
    vmin=-5e5,
    vmax=5e5,
    extent=(r[0, 0], r[0, -1], t[-1], -t[-1]),
)
ax2.set_title(r"$G_{est}$")
axs[1].set_xlabel(r"$x_R$")
axs[1].set_ylabel(r"$t$")
ax2.axis("tight")
ax2.set_ylim(1.2, 0)
ax3.plot(Gsub[:, nr // 2] / Gsub.max(), t, "r", lw=5)
ax3.plot(g_inv_tot[nr // 2, nt - 1 :] / g_inv_tot.max(), t, "k", lw=3)
ax3.set_ylim(1.2, 0)
fig.tight_layout()
  • $p_0^-$, $g^-$, $g^+$
  • $G_{true}$, $G_{est}$

Note that Marchenko redatuming can also be applied simultaneously to multiple subsurface points. Use pylops.waveeqprocessing.Marchenko.apply_multiplepoints instead of pylops.waveeqprocessing.Marchenko.apply_onepoint.

Total running time of the script: ( 0 minutes 7.702 seconds)

Gallery generated by Sphinx-Gallery

11. Radon filtering

In this example we will be taking advantage of the pylops.signalprocessing.Radon2D operator to perform filtering of unwanted events from a seismic data. For those of you not familiar with seismic data, let’s imagine that we have a data composed of a certain number of flat events and a parabolic event , we are after a transform that allows us to separate such an event from the others and filter it out. Those of you with a geophysics background may immediately realize this is the case of seismic angle (or offset) gathers after migration and those events with parabolic moveout are generally residual multiples that we would like to suppress prior to performing further analysis of our data.

The Radon transform is actually a very good transform to perform such a separation. We can thus devise a simple workflow that takes our data as input, applies a Radon transform, filters some of the events out and goes back to the original domain.

import matplotlib.pyplot as plt
import numpy as np

import pylops
from pylops.utils.wavelets import ricker

plt.close("all")
np.random.seed(0)

Let’s first create a data composed on 3 linear events and a parabolic event.

par = {"ox": 0, "dx": 2, "nx": 121, "ot": 0, "dt": 0.004, "nt": 100, "f0": 30}

# linear events
v = 1500  # m/s
t0 = [0.1, 0.2, 0.3]  # s
theta = [0, 0, 0]
amp = [1.0, -2, 0.5]

# parabolic event
tp0 = [0.13]  # s
px = [0]  # s/m
pxx = [5e-7]  # s²/m²
ampp = [0.7]

# create axis
taxis, taxis2, xaxis, yaxis = pylops.utils.seismicevents.makeaxis(par)

# create wavelet
wav = ricker(taxis[:41], f0=par["f0"])[0]

# generate model
y = (
    pylops.utils.seismicevents.linear2d(xaxis, taxis, v, t0, theta, amp, wav)[1]
    + pylops.utils.seismicevents.parabolic2d(xaxis, taxis, tp0, px, pxx, ampp, wav)[1]
)

We can now create the pylops.signalprocessing.Radon2D operator. We also apply its adjoint to the data to obtain a representation of those 3 linear events overlapping to a parabolic event in the Radon domain. Similarly, we feed the operator to a sparse solver like pylops.optimization.sparsity.FISTA to obtain a sparse represention of the data in the Radon domain. At this point we try to filter out the unwanted event. We can see how this is much easier for the sparse transform as each event has a much more compact representation in the Radon domain than for the adjoint transform.

# radon operator
npx = 61
pxmax = 5e-4  # s/m
px = np.linspace(-pxmax, pxmax, npx)

Rop = pylops.signalprocessing.Radon2D(
    taxis, xaxis, px, kind="linear", interp="nearest", centeredh=False, dtype="float64"
)

# adjoint Radon transform
xadj = Rop.H * y

# sparse Radon transform
xinv, niter, cost = pylops.optimization.sparsity.fista(
    Rop, y.ravel(), niter=15, eps=1e1
)
xinv = xinv.reshape(Rop.dims)

# filtering
xfilt = np.zeros_like(xadj)
xfilt[npx // 2 - 3 : npx // 2 + 4] = xadj[npx // 2 - 3 : npx // 2 + 4]

yfilt = Rop * xfilt

# filtering on sparse transform
xinvfilt = np.zeros_like(xinv)
xinvfilt[npx // 2 - 3 : npx // 2 + 4] = xinv[npx // 2 - 3 : npx // 2 + 4]

yinvfilt = Rop * xinvfilt

Finally we visualize our results.

pclip = 0.7
fig, axs = plt.subplots(1, 5, sharey=True, figsize=(12, 5))
axs[0].imshow(
    y.T,
    cmap="gray",
    vmin=-pclip * np.abs(y).max(),
    vmax=pclip * np.abs(y).max(),
    extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]),
)
axs[0].set(xlabel="$x$ [m]", ylabel="$t$ [s]", title="Data")
axs[0].axis("tight")
axs[1].imshow(
    xadj.T,
    cmap="gray",
    vmin=-pclip * np.abs(xadj).max(),
    vmax=pclip * np.abs(xadj).max(),
    extent=(px[0], px[-1], taxis[-1], taxis[0]),
)
axs[1].axvline(px[npx // 2 - 3], color="r", linestyle="--")
axs[1].axvline(px[npx // 2 + 3], color="r", linestyle="--")
axs[1].set(xlabel="$p$ [s/m]", title="Radon")
axs[1].axis("tight")
axs[2].imshow(
    yfilt.T,
    cmap="gray",
    vmin=-pclip * np.abs(yfilt).max(),
    vmax=pclip * np.abs(yfilt).max(),
    extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]),
)
axs[2].set(xlabel="$x$ [m]", title="Filtered data")
axs[2].axis("tight")
axs[3].imshow(
    xinv.T,
    cmap="gray",
    vmin=-pclip * np.abs(xinv).max(),
    vmax=pclip * np.abs(xinv).max(),
    extent=(px[0], px[-1], taxis[-1], taxis[0]),
)
axs[3].axvline(px[npx // 2 - 3], color="r", linestyle="--")
axs[3].axvline(px[npx // 2 + 3], color="r", linestyle="--")
axs[3].set(xlabel="$p$ [s/m]", title="Sparse Radon")
axs[3].axis("tight")
axs[4].imshow(
    yinvfilt.T,
    cmap="gray",
    vmin=-pclip * np.abs(y).max(),
    vmax=pclip * np.abs(y).max(),
    extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]),
)

axs[4].set(xlabel="$x$ [m]", title="Sparse filtered data")
axs[4].axis("tight")
plt.tight_layout()
Data, Radon, Filtered data, Sparse Radon, Sparse filtered data

As expected, the Radon domain is a suitable domain for this type of filtering and the sparse transform improves the ability to filter out parabolic events with small curvature.

On the other hand, it is important to note that we have not been able to correctly preserve the amplitudes of each event. This is because the sparse Radon transform can only identify a sparsest response that explain the data within a certain threshold. For this reason a more suitable approach for preserving amplitudes could be to apply a parabolic Raodn transform with the aim of reconstructing only the unwanted event and apply an adaptive subtraction between the input data and the reconstructed unwanted event.

Total running time of the script: ( 0 minutes 26.034 seconds)

Gallery generated by Sphinx-Gallery

12. Seismic regularization

The problem of seismic data regularization (or interpolation) is a very simple one to write, yet ill-posed and very hard to solve.

The forward modelling operator is a simple pylops.Restriction operator which is applied along the spatial direction(s).

\[\mathbf{y} = \mathbf{R} \mathbf{x}\]

Here \(\mathbf{y} = [\mathbf{y}_{R_1}^T, \mathbf{y}_{R_2}^T,\ldots, \mathbf{y}_{R_N^T}]^T\) where each vector \(\mathbf{y}_{R_i}\) contains all time samples recorded in the seismic data at the specific receiver \(R_i\). Similarly, \(\mathbf{x} = [\mathbf{x}_{r_1}^T, \mathbf{x}_{r_2}^T,\ldots, \mathbf{x}_{r_M}^T]\), contains all traces at the regularly and finely sampled receiver locations \(r_i\).

By inverting such an equation we can create a regularized data with densely and regularly spatial direction(s).

import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import convolve

import pylops
from pylops.utils.seismicevents import linear2d, makeaxis
from pylops.utils.wavelets import ricker

np.random.seed(0)
plt.close("all")

Let’s start by creating a very simple 2d data composed of 3 linear events input parameters

par = {"ox": 0, "dx": 2, "nx": 70, "ot": 0, "dt": 0.004, "nt": 80, "f0": 20}

v = 1500
t0_m = [0.1, 0.2, 0.28]
theta_m = [0, 30, -80]
phi_m = [0]
amp_m = [1.0, -2, 0.5]

# axis
taxis, t2, xaxis, y = makeaxis(par)

# wavelet
wav = ricker(taxis[:41], f0=par["f0"])[0]

# model
_, x = linear2d(xaxis, taxis, v, t0_m, theta_m, amp_m, wav)

We can now define the spatial locations along which the data has been sampled. In this specific example we will assume that we have access only to 40% of the ‘original’ locations.

perc_subsampling = 0.6
nxsub = int(np.round(par["nx"] * perc_subsampling))

iava = np.sort(np.random.permutation(np.arange(par["nx"]))[:nxsub])

# restriction operator
Rop = pylops.Restriction((par["nx"], par["nt"]), iava, axis=0, dtype="float64")

# data
y = Rop * x.ravel()
y = y.reshape(nxsub, par["nt"])

# mask
ymask = Rop.mask(x.ravel())

# inverse
xinv = Rop / y.ravel()
xinv = xinv.reshape(par["nx"], par["nt"])

fig, axs = plt.subplots(1, 2, sharey=True, figsize=(5, 4))
axs[0].imshow(
    x.T, cmap="gray", vmin=-2, vmax=2, extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0])
)
axs[0].set_title("Model")
axs[0].axis("tight")
axs[1].imshow(
    ymask.T,
    cmap="gray",
    vmin=-2,
    vmax=2,
    extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]),
)
axs[1].set_title("Masked model")
axs[1].axis("tight")
plt.tight_layout()
Model, Masked model

As we can see, inverting the restriction operator is not possible without adding any prior information into the inverse problem. In the following we will consider two possible routes:

  • regularized inversion with second derivative along the spatial axis

    \[J = \|\mathbf{y} - \mathbf{R} \mathbf{x}\|_2 + \epsilon_\nabla ^2 \|\nabla \mathbf{x}\|_2\]
  • sparsity-promoting inversion with pylops.FFT2 operator used as sparsyfing transform

    \[J = \|\mathbf{y} - \mathbf{R} \mathbf{F}^H \mathbf{x}\|_2 + \epsilon \|\mathbf{F}^H \mathbf{x}\|_1\]
# smooth inversion
D2op = pylops.SecondDerivative((par["nx"], par["nt"]), axis=0, dtype="float64")

xsmooth, _, _ = pylops.waveeqprocessing.SeismicInterpolation(
    y,
    par["nx"],
    iava,
    kind="spatial",
    **dict(epsRs=[np.sqrt(0.1)], damp=np.sqrt(1e-4), iter_lim=50, show=0)
)

# sparse inversion with FFT2
nfft = 2**8
FFTop = pylops.signalprocessing.FFT2D(
    dims=[par["nx"], par["nt"]], nffts=[nfft, nfft], sampling=[par["dx"], par["dt"]]
)
X = FFTop * x.ravel()
X = np.reshape(X, (nfft, nfft))

xl1, Xl1, cost = pylops.waveeqprocessing.SeismicInterpolation(
    y,
    par["nx"],
    iava,
    kind="fk",
    nffts=(nfft, nfft),
    sampling=(par["dx"], par["dt"]),
    **dict(niter=50, eps=1e-1)
)

fig, axs = plt.subplots(1, 4, sharey=True, figsize=(13, 4))
axs[0].imshow(
    x.T, cmap="gray", vmin=-2, vmax=2, extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0])
)
axs[0].set_title("Model")
axs[0].axis("tight")
axs[1].imshow(
    ymask.T,
    cmap="gray",
    vmin=-2,
    vmax=2,
    extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]),
)
axs[1].set_title("Masked model")
axs[1].axis("tight")
axs[2].imshow(
    xsmooth.T,
    cmap="gray",
    vmin=-2,
    vmax=2,
    extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]),
)
axs[2].set_title("Smoothed model")
axs[2].axis("tight")
axs[3].imshow(
    xl1.T,
    cmap="gray",
    vmin=-2,
    vmax=2,
    extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]),
)
axs[3].set_title("L1 model")
axs[3].axis("tight")

fig, axs = plt.subplots(1, 3, figsize=(10, 2))
axs[0].imshow(
    np.fft.fftshift(np.abs(X[:, : nfft // 2 - 1]), axes=0).T,
    extent=(
        np.fft.fftshift(FFTop.f1)[0],
        np.fft.fftshift(FFTop.f1)[-1],
        FFTop.f2[nfft // 2 - 1],
        FFTop.f2[0],
    ),
)
axs[0].set_title("Model in f-k domain")
axs[0].axis("tight")
axs[0].set_xlim(-0.1, 0.1)
axs[0].set_ylim(50, 0)
axs[1].imshow(
    np.fft.fftshift(np.abs(Xl1[:, : nfft // 2 - 1]), axes=0).T,
    extent=(
        np.fft.fftshift(FFTop.f1)[0],
        np.fft.fftshift(FFTop.f1)[-1],
        FFTop.f2[nfft // 2 - 1],
        FFTop.f2[0],
    ),
)
axs[1].set_title("Reconstructed model in f-k domain")
axs[1].axis("tight")
axs[1].set_xlim(-0.1, 0.1)
axs[1].set_ylim(50, 0)
axs[2].plot(cost, "k", lw=3)
axs[2].set_title("FISTA convergence")
plt.tight_layout()
  • Model, Masked model, Smoothed model, L1 model
  • Model in f-k domain, Reconstructed model in f-k domain, FISTA convergence

We see how adding prior information to the inversion can help improving the estimate of the regularized seismic data. Nevertheless, in both cases the reconstructed data is not perfect. A better sparsyfing transform could in fact be chosen here to be the linear pylops.signalprocessing.Radon2D transform in spite of the pylops.FFT2 transform.

npx = 40
pxmax = 1e-3
px = np.linspace(-pxmax, pxmax, npx)
Radop = pylops.signalprocessing.Radon2D(taxis, xaxis, px, engine="numba")

RRop = Rop * Radop

# adjoint
Xadj_fromx = Radop.H * x.ravel()
Xadj_fromx = Xadj_fromx.reshape(npx, par["nt"])

Xadj = RRop.H * y.ravel()
Xadj = Xadj.reshape(npx, par["nt"])

# L1 inverse
xl1, Xl1, cost = pylops.waveeqprocessing.SeismicInterpolation(
    y,
    par["nx"],
    iava,
    kind="radon-linear",
    spataxis=xaxis,
    taxis=taxis,
    paxis=px,
    centeredh=True,
    **dict(niter=50, eps=1e-1)
)

fig, axs = plt.subplots(2, 3, sharey=True, figsize=(12, 7))
axs[0][0].imshow(
    x.T, cmap="gray", vmin=-2, vmax=2, extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0])
)
axs[0][0].set_title("Data", fontsize=12)
axs[0][0].axis("tight")
axs[0][1].imshow(
    ymask.T,
    cmap="gray",
    vmin=-2,
    vmax=2,
    extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]),
)
axs[0][1].set_title("Masked data", fontsize=12)
axs[0][1].axis("tight")
axs[0][2].imshow(
    xl1.T,
    cmap="gray",
    vmin=-2,
    vmax=2,
    extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]),
)
axs[0][2].set_title("Reconstructed data", fontsize=12)
axs[0][2].axis("tight")
axs[1][0].imshow(
    Xadj_fromx.T,
    cmap="gray",
    vmin=-70,
    vmax=70,
    extent=(px[0], px[-1], taxis[-1], taxis[0]),
)
axs[1][0].set_title("Adj. Radon on data", fontsize=12)
axs[1][0].axis("tight")
axs[1][1].imshow(
    Xadj.T, cmap="gray", vmin=-50, vmax=50, extent=(px[0], px[-1], taxis[-1], taxis[0])
)
axs[1][1].set_title("Adj. Radon on subsampled data", fontsize=12)
axs[1][1].axis("tight")
axs[1][2].imshow(
    Xl1.T, cmap="gray", vmin=-0.2, vmax=0.2, extent=(px[0], px[-1], taxis[-1], taxis[0])
)
axs[1][2].set_title("Inverse Radon on subsampled data", fontsize=12)
axs[1][2].axis("tight")
plt.tight_layout()
Data, Masked data, Reconstructed data, Adj. Radon on data, Adj. Radon on subsampled data, Inverse Radon on subsampled data

Finally, let’s take now a more realistic dataset. We will use once again the linear pylops.signalprocessing.Radon2D transform but we will take advantnge of the pylops.signalprocessing.Sliding2D operator to perform such a transform locally instead of globally to the entire dataset.

inputfile = "../testdata/marchenko/input.npz"
inputdata = np.load(inputfile)

x = inputdata["R"][50, :, ::2]
x = x / np.abs(x).max()
taxis, xaxis = inputdata["t"][::2], inputdata["r"][0]

par = {}
par["nx"], par["nt"] = x.shape
par["dx"] = inputdata["r"][0, 1] - inputdata["r"][0, 0]
par["dt"] = inputdata["t"][1] - inputdata["t"][0]

# add wavelet
wav = inputdata["wav"][::2]
wav_c = np.argmax(wav)
x = np.apply_along_axis(convolve, 1, x, wav, mode="full")
x = x[:, wav_c:][:, : par["nt"]]

# gain
gain = np.tile((taxis**2)[:, np.newaxis], (1, par["nx"])).T
x = x * gain

# subsampling locations
perc_subsampling = 0.5
Nsub = int(np.round(par["nx"] * perc_subsampling))
iava = np.sort(np.random.permutation(np.arange(par["nx"]))[:Nsub])

# restriction operator
Rop = pylops.Restriction((par["nx"], par["nt"]), iava, axis=0, dtype="float64")

y = Rop * x.ravel()
xadj = Rop.H * y.ravel()

y = y.reshape(Nsub, par["nt"])
xadj = xadj.reshape(par["nx"], par["nt"])

# apply mask
ymask = Rop.mask(x.ravel())

# sliding windows with radon transform
dx = par["dx"]
nwins = 4
nwin = 27
nover = 3
npx = 31
pxmax = 5e-4
px = np.linspace(-pxmax, pxmax, npx)
dimsd = x.shape
dims = (nwins * npx, dimsd[1])

Op = pylops.signalprocessing.Radon2D(
    taxis,
    np.linspace(-par["dx"] * nwin // 2, par["dx"] * nwin // 2, nwin),
    px,
    centeredh=True,
    kind="linear",
    engine="numba",
)
Slidop = pylops.signalprocessing.Sliding2D(
    Op, dims, dimsd, nwin, nover, tapertype="cosine"
)

# adjoint
RSop = Rop * Slidop

Xadj_fromx = Slidop.H * x.ravel()
Xadj_fromx = Xadj_fromx.reshape(npx * nwins, par["nt"])

Xadj = RSop.H * y.ravel()
Xadj = Xadj.reshape(npx * nwins, par["nt"])

# inverse
xl1, Xl1, _ = pylops.waveeqprocessing.SeismicInterpolation(
    y,
    par["nx"],
    iava,
    kind="sliding",
    spataxis=xaxis,
    taxis=taxis,
    paxis=px,
    nwins=nwins,
    nwin=nwin,
    nover=nover,
    **dict(niter=50, eps=1e-2)
)

fig, axs = plt.subplots(2, 3, sharey=True, figsize=(12, 14))
axs[0][0].imshow(
    x.T,
    cmap="gray",
    vmin=-0.1,
    vmax=0.1,
    extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]),
)
axs[0][0].set_title("Data")
axs[0][0].axis("tight")
axs[0][1].imshow(
    ymask.T,
    cmap="gray",
    vmin=-0.1,
    vmax=0.1,
    extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]),
)
axs[0][1].set_title("Masked data")
axs[0][1].axis("tight")
axs[0][2].imshow(
    xl1.T,
    cmap="gray",
    vmin=-0.1,
    vmax=0.1,
    extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]),
)
axs[0][2].set_title("Reconstructed data")
axs[0][2].axis("tight")
axs[1][0].imshow(
    Xadj_fromx.T,
    cmap="gray",
    vmin=-1,
    vmax=1,
    extent=(px[0], px[-1], taxis[-1], taxis[0]),
)
axs[1][0].set_title("Adjoint Radon on data")
axs[1][0].axis("tight")
axs[1][1].imshow(
    Xadj.T,
    cmap="gray",
    vmin=-0.6,
    vmax=0.6,
    extent=(px[0], px[-1], taxis[-1], taxis[0]),
)
axs[1][1].set_title("Adjoint Radon on subsampled data")
axs[1][1].axis("tight")
axs[1][2].imshow(
    Xl1.T,
    cmap="gray",
    vmin=-0.03,
    vmax=0.03,
    extent=(px[0], px[-1], taxis[-1], taxis[0]),
)
axs[1][2].set_title("Inverse Radon on subsampled data")
axs[1][2].axis("tight")
plt.tight_layout()
Data, Masked data, Reconstructed data, Adjoint Radon on data, Adjoint Radon on subsampled data, Inverse Radon on subsampled data

As expected the linear pylops.signalprocessing.Radon2D is able to locally explain events in the input data and leads to a satisfactory recovery. Note that increasing the number of iterations and sliding windows can further refine the result, especially the accuracy of weak events, as shown in this companion notebook.

Total running time of the script: ( 0 minutes 8.310 seconds)

Gallery generated by Sphinx-Gallery

13. Deghosting

Single-component seismic data can be decomposed in their up- and down-going constituents in a model driven fashion. This task can be achieved by defining an f-k propagator (or ghost model) and solving an inverse problem as described in pylops.waveeqprocessing.Deghosting.

import matplotlib.pyplot as plt

# sphinx_gallery_thumbnail_number = 3
import numpy as np
from scipy.sparse.linalg import lsqr

import pylops

np.random.seed(0)
plt.close("all")

Let’s start by loading the input dataset and geometry

inputfile = "../testdata/updown/input.npz"
inputdata = np.load(inputfile)

vel_sep = 2400.0  # velocity at separation level
clip = 1e-1  # plotting clip


# Receivers
r = inputdata["r"]
nr = r.shape[1]
dr = r[0, 1] - r[0, 0]

# Sources
s = inputdata["s"]

# Model
rho = inputdata["rho"]

# Axes
t = inputdata["t"]
nt, dt = len(t), t[1] - t[0]
x, z = inputdata["x"], inputdata["z"]
dx, dz = x[1] - x[0], z[1] - z[0]

# Data
p = inputdata["p"].T
p /= p.max()

fig = plt.figure(figsize=(9, 4))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=4)
ax2 = plt.subplot2grid((1, 5), (0, 4))
ax1.imshow(rho, cmap="gray", extent=(x[0], x[-1], z[-1], z[0]))
ax1.scatter(r[0, ::5], r[1, ::5], marker="v", s=150, c="b", edgecolors="k")
ax1.scatter(s[0], s[1], marker="*", s=250, c="r", edgecolors="k")
ax1.axis("tight")
ax1.set_xlabel("x [m]")
ax1.set_ylabel("y [m]")
ax1.set_title("Model and Geometry")
ax1.set_xlim(x[0], x[-1])
ax1.set_ylim(z[-1], z[0])
ax2.plot(rho[:, len(x) // 2], z, "k", lw=2)
ax2.set_ylim(z[-1], z[0])
ax2.set_yticks([])
plt.tight_layout()
Model and Geometry

To be able to deghost the input dataset, we need to remove its direct arrival. In this example we will create a mask based on the analytical traveltime of the direct arrival.

direct = np.sqrt(np.sum((s[:, np.newaxis] - r) ** 2, axis=0)) / vel_sep

# Window
off = 0.035
direct_off = direct + off
win = np.zeros((nt, nr))
iwin = np.round(direct_off / dt).astype(int)
for i in range(nr):
    win[iwin[i] :, i] = 1

fig, axs = plt.subplots(1, 2, sharey=True, figsize=(8, 7))
axs[0].imshow(
    p.T,
    cmap="gray",
    vmin=-clip * np.abs(p).max(),
    vmax=clip * np.abs(p).max(),
    extent=(r[0, 0], r[0, -1], t[-1], t[0]),
)
axs[0].plot(r[0], direct_off, "r", lw=2)
axs[0].set_title(r"$P$")
axs[0].axis("tight")
axs[1].imshow(
    win * p.T,
    cmap="gray",
    vmin=-clip * np.abs(p).max(),
    vmax=clip * np.abs(p).max(),
    extent=(r[0, 0], r[0, -1], t[-1], t[0]),
)
axs[1].set_title(r"Windowed $P$")
axs[1].axis("tight")
axs[1].set_ylim(1, 0)
plt.tight_layout()
$P$, Windowed $P$

We can now perform deghosting

pup, pdown = pylops.waveeqprocessing.Deghosting(
    p.T,
    nt,
    nr,
    dt,
    dr,
    vel_sep,
    r[1, 0] + dz,
    win=win,
    npad=5,
    ntaper=11,
    solver=lsqr,
    dottest=False,
    dtype="complex128",
    **dict(damp=1e-10, iter_lim=60)
)

fig, axs = plt.subplots(1, 3, sharey=True, figsize=(12, 7))
axs[0].imshow(
    p.T,
    cmap="gray",
    vmin=-clip * np.abs(p).max(),
    vmax=clip * np.abs(p).max(),
    extent=(r[0, 0], r[0, -1], t[-1], t[0]),
)
axs[0].set_title(r"$P$")
axs[0].axis("tight")
axs[1].imshow(
    pup,
    cmap="gray",
    vmin=-clip * np.abs(p).max(),
    vmax=clip * np.abs(p).max(),
    extent=(r[0, 0], r[0, -1], t[-1], t[0]),
)
axs[1].set_title(r"$P^-$")
axs[1].axis("tight")
axs[2].imshow(
    pdown,
    cmap="gray",
    vmin=-clip * np.abs(p).max(),
    vmax=clip * np.abs(p).max(),
    extent=(r[0, 0], r[0, -1], t[-1], t[0]),
)
axs[2].set_title(r"$P^+$")
axs[2].axis("tight")
axs[2].set_ylim(1, 0)

plt.figure(figsize=(14, 3))
plt.plot(t, p[nr // 2], "k", lw=2, label=r"$p$")
plt.plot(t, pup[:, nr // 2], "r", lw=2, label=r"$p^-$")
plt.xlim(0, t[200])
plt.ylim(-0.2, 0.2)
plt.legend()
plt.tight_layout()

plt.figure(figsize=(14, 3))
plt.plot(t, pdown[:, nr // 2], "b", lw=2, label=r"$p^+$")
plt.plot(t, pup[:, nr // 2], "r", lw=2, label=r"$p^-$")
plt.xlim(0, t[200])
plt.ylim(-0.2, 0.2)
plt.legend()
plt.tight_layout()
  • $P$, $P^-$, $P^+$
  • deghosting
  • deghosting

To see more examples head over to the following notebook: notebook1.

Total running time of the script: ( 0 minutes 6.874 seconds)

Gallery generated by Sphinx-Gallery

14. Seismic wavefield decomposition

Multi-component seismic data can be decomposed in their up- and down-going constituents in a purely data driven fashion. This task can be accurately achieved by linearly combining the input pressure and particle velocity data in the frequency-wavenumber described in details in pylops.waveeqprocessing.UpDownComposition2D and pylops.waveeqprocessing.WavefieldDecomposition.

In this tutorial we will consider a simple synthetic data composed of six events (three up-going and three down-going). We will first combine them to create pressure and particle velocity data and then show how we can retrieve their directional constituents both by directly combining the input data as well as by setting an inverse problem. The latter approach results vital in case of spatial aliasing, as applying simple scaled summation in the frequency-wavenumber would result in sub-optimal decomposition due to the superposition of different frequency-wavenumber pairs at some (aliased) locations.

import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import filtfilt

import pylops
from pylops.utils.seismicevents import hyperbolic2d, makeaxis
from pylops.utils.wavelets import ricker

np.random.seed(0)
plt.close("all")

Let’s first the input up- and down-going wavefields

par = {"ox": -220, "dx": 5, "nx": 89, "ot": 0, "dt": 0.004, "nt": 200, "f0": 40}

t0_plus = np.array([0.2, 0.5, 0.7])
t0_minus = t0_plus + 0.04
vrms = np.array([1400.0, 1500.0, 2000.0])
amp = np.array([1.0, -0.6, 0.5])
vel_sep = 1000.0  # velocity at separation level
rho_sep = 1000.0  # density at separation level

# Create axis
t, t2, x, y = makeaxis(par)

# Create wavelet
wav = ricker(t[:41], f0=par["f0"])[0]

# Create data
_, p_minus = hyperbolic2d(x, t, t0_minus, vrms, amp, wav)
_, p_plus = hyperbolic2d(x, t, t0_plus, vrms, amp, wav)

We can now combine them to create pressure and particle velocity data

critical = 1.1
ntaper = 51
nfft = 2**10

# 2d fft operator
FFTop = pylops.signalprocessing.FFT2D(
    dims=[par["nx"], par["nt"]], nffts=[nfft, nfft], sampling=[par["dx"], par["dt"]]
)

# obliquity factor
[Kx, F] = np.meshgrid(FFTop.f1, FFTop.f2, indexing="ij")
k = F / vel_sep
Kz = np.sqrt((k**2 - Kx**2).astype(np.complex128))
Kz[np.isnan(Kz)] = 0
OBL = rho_sep * (np.abs(F) / Kz)
OBL[Kz == 0] = 0

mask = np.abs(Kx) < critical * np.abs(F) / vel_sep
OBL *= mask
OBL = filtfilt(np.ones(ntaper) / float(ntaper), 1, OBL, axis=0)
OBL = filtfilt(np.ones(ntaper) / float(ntaper), 1, OBL, axis=1)

# composition operator
UPop = pylops.waveeqprocessing.UpDownComposition2D(
    par["nt"],
    par["nx"],
    par["dt"],
    par["dx"],
    rho_sep,
    vel_sep,
    nffts=(nfft, nfft),
    critical=critical * 100.0,
    ntaper=ntaper,
    dtype="complex128",
)

# wavefield modelling
d = UPop * np.concatenate((p_plus.ravel(), p_minus.ravel())).ravel()
d = np.real(d.reshape(2 * par["nx"], par["nt"]))
p, vz = d[: par["nx"]], d[par["nx"] :]

# obliquity scaled vz
VZ = FFTop * vz.ravel()
VZ = VZ.reshape(nfft, nfft)

VZ_obl = OBL * VZ
vz_obl = FFTop.H * VZ_obl.ravel()
vz_obl = np.real(vz_obl.reshape(par["nx"], par["nt"]))

fig, axs = plt.subplots(1, 4, figsize=(10, 5))
axs[0].imshow(
    p.T,
    aspect="auto",
    vmin=-1,
    vmax=1,
    interpolation="nearest",
    cmap="gray",
    extent=(x.min(), x.max(), t.max(), t.min()),
)
axs[0].set_title(r"$p$", fontsize=15)
axs[0].set_xlabel("x")
axs[0].set_ylabel("t")
axs[1].imshow(
    vz_obl.T,
    aspect="auto",
    vmin=-1,
    vmax=1,
    interpolation="nearest",
    cmap="gray",
    extent=(x.min(), x.max(), t.max(), t.min()),
)
axs[1].set_title(r"$v_z^{obl}$", fontsize=15)
axs[1].set_xlabel("x")
axs[1].set_ylabel("t")
axs[2].imshow(
    p_plus.T,
    aspect="auto",
    vmin=-1,
    vmax=1,
    interpolation="nearest",
    cmap="gray",
    extent=(x.min(), x.max(), t.max(), t.min()),
)
axs[2].set_title(r"$p^+$", fontsize=15)
axs[2].set_xlabel("x")
axs[2].set_ylabel("t")
axs[3].imshow(
    p_minus.T,
    aspect="auto",
    interpolation="nearest",
    cmap="gray",
    extent=(x.min(), x.max(), t.max(), t.min()),
    vmin=-1,
    vmax=1,
)
axs[3].set_title(r"$p^-$", fontsize=15)
axs[3].set_xlabel("x")
axs[3].set_ylabel("t")
plt.tight_layout()
$p$, $v_z^{obl}$, $p^+$, $p^-$

Wavefield separation is first performed using the analytical expression for combining pressure and particle velocity data in the wavenumber-frequency domain

pup_sep, pdown_sep = pylops.waveeqprocessing.WavefieldDecomposition(
    p,
    vz,
    par["nt"],
    par["nx"],
    par["dt"],
    par["dx"],
    rho_sep,
    vel_sep,
    nffts=(nfft, nfft),
    kind="analytical",
    critical=critical * 100,
    ntaper=ntaper,
    dtype="complex128",
)
fig = plt.figure(figsize=(12, 5))
axs0 = plt.subplot2grid((2, 5), (0, 0), rowspan=2)
axs1 = plt.subplot2grid((2, 5), (0, 1), rowspan=2)
axs2 = plt.subplot2grid((2, 5), (0, 2), colspan=3)
axs3 = plt.subplot2grid((2, 5), (1, 2), colspan=3)
axs0.imshow(
    pup_sep.T, cmap="gray", vmin=-1, vmax=1, extent=(x.min(), x.max(), t.max(), t.min())
)
axs0.set_title(r"$p^-$ analytical")
axs0.axis("tight")
axs1.imshow(
    pdown_sep.T,
    cmap="gray",
    vmin=-1,
    vmax=1,
    extent=(x.min(), x.max(), t.max(), t.min()),
)
axs1.set_title(r"$p^+$ analytical")
axs1.axis("tight")
axs2.plot(t, p[par["nx"] // 2], "r", lw=2, label=r"$p$")
axs2.plot(t, vz_obl[par["nx"] // 2], "--b", lw=2, label=r"$v_z^{obl}$")
axs2.set_ylim(-1, 1)
axs2.set_title("Data at x=%.2f" % x[par["nx"] // 2])
axs2.set_xlabel("t [s]")
axs2.legend()
axs3.plot(t, pup_sep[par["nx"] // 2], "r", lw=2, label=r"$p^-$ ana")
axs3.plot(t, pdown_sep[par["nx"] // 2], "--b", lw=2, label=r"$p^+$ ana")
axs3.set_title("Separated wavefields at x=%.2f" % x[par["nx"] // 2])
axs3.set_xlabel("t [s]")
axs3.set_ylim(-1, 1)
axs3.legend()
plt.tight_layout()
$p^-$ analytical, $p^+$ analytical, Data at x=0.00, Separated wavefields at x=0.00

We repeat the same exercise but this time we invert the composition operator pylops.waveeqprocessing.UpDownComposition2D

pup_inv, pdown_inv = pylops.waveeqprocessing.WavefieldDecomposition(
    p,
    vz,
    par["nt"],
    par["nx"],
    par["dt"],
    par["dx"],
    rho_sep,
    vel_sep,
    nffts=(nfft, nfft),
    kind="inverse",
    critical=critical * 100,
    ntaper=ntaper,
    scaling=1.0 / vz.max(),
    dtype="complex128",
    **dict(damp=1e-10, iter_lim=20)
)

fig = plt.figure(figsize=(12, 5))
axs0 = plt.subplot2grid((2, 5), (0, 0), rowspan=2)
axs1 = plt.subplot2grid((2, 5), (0, 1), rowspan=2)
axs2 = plt.subplot2grid((2, 5), (0, 2), colspan=3)
axs3 = plt.subplot2grid((2, 5), (1, 2), colspan=3)
axs0.imshow(
    pup_inv.T, cmap="gray", vmin=-1, vmax=1, extent=(x.min(), x.max(), t.max(), t.min())
)
axs0.set_title(r"$p^-$ inverse")
axs0.axis("tight")
axs1.imshow(
    pdown_inv.T,
    cmap="gray",
    vmin=-1,
    vmax=1,
    extent=(x.min(), x.max(), t.max(), t.min()),
)
axs1.set_title(r"$p^+$ inverse")
axs1.axis("tight")
axs2.plot(t, p[par["nx"] // 2], "r", lw=2, label=r"$p$")
axs2.plot(t, vz_obl[par["nx"] // 2], "--b", lw=2, label=r"$v_z^{obl}$")
axs2.set_ylim(-1, 1)
axs2.set_title("Data at x=%.2f" % x[par["nx"] // 2])
axs2.set_xlabel("t [s]")
axs2.legend()
axs3.plot(t, pup_inv[par["nx"] // 2], "r", lw=2, label=r"$p^-$ inv")
axs3.plot(t, pdown_inv[par["nx"] // 2], "--b", lw=2, label=r"$p^+$ inv")
axs3.set_title("Separated wavefields at x=%.2f" % x[par["nx"] // 2])
axs3.set_xlabel("t [s]")
axs3.set_ylim(-1, 1)
axs3.legend()
plt.tight_layout()
$p^-$ inverse, $p^+$ inverse, Data at x=0.00, Separated wavefields at x=0.00

The up- and down-going constituents have been succesfully separated in both cases. Finally, we use the pylops.waveeqprocessing.UpDownComposition2D operator to reconstruct the particle velocity wavefield from its up- and down-going pressure constituents

PtoVop = pylops.waveeqprocessing.PressureToVelocity(
    par["nt"],
    par["nx"],
    par["dt"],
    par["dx"],
    rho_sep,
    vel_sep,
    nffts=(nfft, nfft),
    critical=critical * 100.0,
    ntaper=ntaper,
    topressure=False,
)

vdown_rec = (PtoVop * pdown_inv.ravel()).reshape(par["nx"], par["nt"])
vup_rec = (PtoVop * pup_inv.ravel()).reshape(par["nx"], par["nt"])
vz_rec = np.real(vdown_rec - vup_rec)

fig, axs = plt.subplots(1, 3, figsize=(13, 6))
axs[0].imshow(
    vz.T,
    cmap="gray",
    vmin=-1e-6,
    vmax=1e-6,
    extent=(x.min(), x.max(), t.max(), t.min()),
)
axs[0].set_title(r"$vz$")
axs[0].axis("tight")
axs[1].imshow(
    vz_rec.T, cmap="gray", vmin=-1e-6, vmax=1e-6, extent=(x.min(), x.max(), t[-1], t[0])
)
axs[1].set_title(r"$vz rec$")
axs[1].axis("tight")
axs[2].imshow(
    vz.T - vz_rec.T,
    cmap="gray",
    vmin=-1e-6,
    vmax=1e-6,
    extent=(x.min(), x.max(), t[-1], t[0]),
)
axs[2].set_title(r"$error$")
axs[2].axis("tight")
plt.tight_layout()
$vz$, $vz rec$, $error$

To see more examples, including applying wavefield separation and regularization simultaneously, as well as 3D examples, head over to the following notebooks: notebook1 and notebook2

Total running time of the script: ( 0 minutes 15.497 seconds)

Gallery generated by Sphinx-Gallery

15. Least-squares migration

Seismic migration is the process by which seismic data are manipulated to create an image of the subsurface reflectivity.

While traditionally solved as the adjont of the demigration operator, it is becoming more and more common to solve the underlying inverse problem in the quest for more accurate and detailed subsurface images.

Indipendently of the choice of the modelling operator (i.e., ray-based or full wavefield-based), the demigration/migration process can be expressed as a linear operator of such a kind:

\[d(\mathbf{x_r}, \mathbf{x_s}, t) = w(t) * \int\limits_V G(\mathbf{x}, \mathbf{x_s}, t) G(\mathbf{x_r}, \mathbf{x}, t) m(\mathbf{x})\,\mathrm{d}\mathbf{x}\]

where \(m(\mathbf{x})\) is the reflectivity at every location in the subsurface, \(G(\mathbf{x}, \mathbf{x_s}, t)\) and \(G(\mathbf{x_r}, \mathbf{x}, t)\) are the Green’s functions from source-to-subsurface-to-receiver and finally \(w(t)\) is the wavelet. Ultimately, while the Green’s functions can be computed in many different ways, solving this system of equations for the reflectivity model is what we generally refer to as Least-squares migration (LSM).

In this tutorial we will consider the most simple scenario where we use an eikonal solver to compute the Green’s functions and show how we can use the pylops.waveeqprocessing.LSM operator to perform LSM.

import matplotlib.pyplot as plt
import numpy as np
from scipy.sparse.linalg import lsqr

import pylops

plt.close("all")
np.random.seed(0)

To start we create a simple model with 2 interfaces

# Velocity Model
nx, nz = 81, 60
dx, dz = 4, 4
x, z = np.arange(nx) * dx, np.arange(nz) * dz
v0 = 1000  # initial velocity
kv = 0.0  # gradient
vel = np.outer(np.ones(nx), v0 + kv * z)

# Reflectivity Model
refl = np.zeros((nx, nz))
refl[:, 30] = -1
refl[:, 50] = 0.5

# Receivers
nr = 11
rx = np.linspace(10 * dx, (nx - 10) * dx, nr)
rz = 20 * np.ones(nr)
recs = np.vstack((rx, rz))
dr = recs[0, 1] - recs[0, 0]

# Sources
ns = 10
sx = np.linspace(dx * 10, (nx - 10) * dx, ns)
sz = 10 * np.ones(ns)
sources = np.vstack((sx, sz))
ds = sources[0, 1] - sources[0, 0]
plt.figure(figsize=(10, 5))
im = plt.imshow(vel.T, cmap="summer", extent=(x[0], x[-1], z[-1], z[0]))
plt.scatter(recs[0], recs[1], marker="v", s=150, c="b", edgecolors="k")
plt.scatter(sources[0], sources[1], marker="*", s=150, c="r", edgecolors="k")
cb = plt.colorbar(im)
cb.set_label("[m/s]")
plt.axis("tight")
plt.xlabel("x [m]"), plt.ylabel("y [m]")
plt.title("Velocity")
plt.xlim(x[0], x[-1])
plt.tight_layout()

plt.figure(figsize=(10, 5))
im = plt.imshow(refl.T, cmap="gray", extent=(x[0], x[-1], z[-1], z[0]))
plt.scatter(recs[0], recs[1], marker="v", s=150, c="b", edgecolors="k")
plt.scatter(sources[0], sources[1], marker="*", s=150, c="r", edgecolors="k")
plt.colorbar(im)
plt.axis("tight")
plt.xlabel("x [m]"), plt.ylabel("y [m]")
plt.title("Reflectivity")
plt.xlim(x[0], x[-1])
plt.tight_layout()
  • Velocity
  • Reflectivity

We can now create our LSM object and invert for the reflectivity using two different solvers: scipy.sparse.linalg.lsqr (LS solution) and pylops.optimization.sparsity.fista (LS solution with sparse model).

nt = 651
dt = 0.004
t = np.arange(nt) * dt
wav, wavt, wavc = pylops.utils.wavelets.ricker(t[:41], f0=20)


lsm = pylops.waveeqprocessing.LSM(
    z,
    x,
    t,
    sources,
    recs,
    v0,
    wav,
    wavc,
    mode="analytic",
    engine="numba",
)

d = lsm.Demop * refl

madj = lsm.Demop.H * d

minv = lsm.solve(d.ravel(), solver=lsqr, **dict(iter_lim=100))
minv = minv.reshape(nx, nz)

minv_sparse = lsm.solve(
    d.ravel(), solver=pylops.optimization.sparsity.fista, **dict(eps=1e2, niter=100)
)
minv_sparse = minv_sparse.reshape(nx, nz)

# demigration
d = d.reshape(ns, nr, nt)

dadj = lsm.Demop * madj  # (ns * nr, nt)
dadj = dadj.reshape(ns, nr, nt)

dinv = lsm.Demop * minv
dinv = dinv.reshape(ns, nr, nt)

dinv_sparse = lsm.Demop * minv_sparse
dinv_sparse = dinv_sparse.reshape(ns, nr, nt)

# sphinx_gallery_thumbnail_number = 2
fig, axs = plt.subplots(2, 2, figsize=(10, 8))
axs[0][0].imshow(refl.T, cmap="gray", vmin=-1, vmax=1)
axs[0][0].axis("tight")
axs[0][0].set_title(r"$m$")
axs[0][1].imshow(madj.T, cmap="gray", vmin=-madj.max(), vmax=madj.max())
axs[0][1].set_title(r"$m_{adj}$")
axs[0][1].axis("tight")
axs[1][0].imshow(minv.T, cmap="gray", vmin=-1, vmax=1)
axs[1][0].axis("tight")
axs[1][0].set_title(r"$m_{inv}$")
axs[1][1].imshow(minv_sparse.T, cmap="gray", vmin=-1, vmax=1)
axs[1][1].axis("tight")
axs[1][1].set_title(r"$m_{FISTA}$")
plt.tight_layout()

fig, axs = plt.subplots(1, 4, figsize=(10, 4))
axs[0].imshow(d[0, :, :300].T, cmap="gray", vmin=-d.max(), vmax=d.max())
axs[0].set_title(r"$d$")
axs[0].axis("tight")
axs[1].imshow(dadj[0, :, :300].T, cmap="gray", vmin=-dadj.max(), vmax=dadj.max())
axs[1].set_title(r"$d_{adj}$")
axs[1].axis("tight")
axs[2].imshow(dinv[0, :, :300].T, cmap="gray", vmin=-d.max(), vmax=d.max())
axs[2].set_title(r"$d_{inv}$")
axs[2].axis("tight")
axs[3].imshow(dinv_sparse[0, :, :300].T, cmap="gray", vmin=-d.max(), vmax=d.max())
axs[3].set_title(r"$d_{fista}$")
axs[3].axis("tight")
plt.tight_layout()

fig, axs = plt.subplots(1, 4, figsize=(10, 4))
axs[0].imshow(d[ns // 2, :, :300].T, cmap="gray", vmin=-d.max(), vmax=d.max())
axs[0].set_title(r"$d$")
axs[0].axis("tight")
axs[1].imshow(dadj[ns // 2, :, :300].T, cmap="gray", vmin=-dadj.max(), vmax=dadj.max())
axs[1].set_title(r"$d_{adj}$")
axs[1].axis("tight")
axs[2].imshow(dinv[ns // 2, :, :300].T, cmap="gray", vmin=-d.max(), vmax=d.max())
axs[2].set_title(r"$d_{inv}$")
axs[2].axis("tight")
axs[3].imshow(dinv_sparse[ns // 2, :, :300].T, cmap="gray", vmin=-d.max(), vmax=d.max())
axs[3].set_title(r"$d_{fista}$")
axs[3].axis("tight")
plt.tight_layout()
  • $m$, $m_{adj}$, $m_{inv}$, $m_{FISTA}$
  • $d$, $d_{adj}$, $d_{inv}$, $d_{fista}$
  • $d$, $d_{adj}$, $d_{inv}$, $d_{fista}$

This was just a short teaser, for a more advanced set of examples of 2D and 3D traveltime-based LSM head over to this notebook.

Total running time of the script: ( 0 minutes 7.017 seconds)

Gallery generated by Sphinx-Gallery

16. CT Scan Imaging

This tutorial considers a very well-known inverse problem from the field of medical imaging.

We will be using the pylops.signalprocessing.Radon2D operator to model a sinogram, which is a graphic representation of the raw data obtained from a CT scan. The sinogram is further inverted using both a L2 solver and a TV-regularized solver like Split-Bregman.

import matplotlib.pyplot as plt

# sphinx_gallery_thumbnail_number = 2
import numpy as np
from numba import jit

import pylops

plt.close("all")
np.random.seed(10)

Let’s start by loading the Shepp-Logan phantom model. We can then construct the sinogram by providing a custom-made function to the pylops.signalprocessing.Radon2D that samples parametric curves of such a type:

\[t(r,\theta; x) = \tan(90°-\theta)x + \frac{r}{\sin(\theta)}\]

where \(\theta\) is the angle between the x-axis (\(x\)) and the perpendicular to the summation line and \(r\) is the distance from the origin of the summation line.

@jit(nopython=True)
def radoncurve(x, r, theta):
    return (
        (r - ny // 2) / (np.sin(theta) + 1e-15)
        + np.tan(np.pi / 2.0 - theta) * x
        + ny // 2
    )


x = np.load("../testdata/optimization/shepp_logan_phantom.npy").T
x = x / x.max()
nx, ny = x.shape

ntheta = 151
theta = np.linspace(0.0, np.pi, ntheta, endpoint=False)

RLop = pylops.signalprocessing.Radon2D(
    np.arange(ny),
    np.arange(nx),
    theta,
    kind=radoncurve,
    centeredh=True,
    interp=False,
    engine="numba",
    dtype="float64",
)

y = RLop.H * x

We can now first perform the adjoint, which in the medical imaging literature is also referred to as back-projection.

This is the first step of a common reconstruction technique, named filtered back-projection, which simply applies a correction filter in the frequency domain to the adjoint model.

xrec = RLop * y

fig, axs = plt.subplots(1, 3, figsize=(10, 4))
axs[0].imshow(x.T, vmin=0, vmax=1, cmap="gray")
axs[0].set_title("Model")
axs[0].axis("tight")
axs[1].imshow(y.T, cmap="gray")
axs[1].set_title("Data")
axs[1].axis("tight")
axs[2].imshow(xrec.T, cmap="gray")
axs[2].set_title("Adjoint model")
axs[2].axis("tight")
fig.tight_layout()
Model, Data, Adjoint model

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.

Dop = [
    pylops.FirstDerivative(
        (nx, ny), axis=0, edge=True, kind="backward", dtype=np.float64
    ),
    pylops.FirstDerivative(
        (nx, ny), axis=1, edge=True, kind="backward", dtype=np.float64
    ),
]
D2op = pylops.Laplacian(dims=(nx, ny), edge=True, dtype=np.float64)

# L2
xinv_sm = pylops.optimization.leastsquares.regularized_inversion(
    RLop.H, y.ravel(), [D2op], epsRs=[1e1], **dict(iter_lim=20)
)[0]
xinv_sm = np.real(xinv_sm.reshape(nx, ny))

# TV
mu = 1.5
lamda = [1.0, 1.0]
niter = 3
niterinner = 4

xinv = pylops.optimization.sparsity.splitbregman(
    RLop.H,
    y.ravel(),
    Dop,
    niter_outer=niter,
    niter_inner=niterinner,
    mu=mu,
    epsRL1s=lamda,
    tol=1e-4,
    tau=1.0,
    show=False,
    **dict(iter_lim=20, damp=1e-2)
)[0]
xinv = np.real(xinv.reshape(nx, ny))

fig, axs = plt.subplots(1, 3, figsize=(10, 4))
axs[0].imshow(x.T, vmin=0, vmax=1, cmap="gray")
axs[0].set_title("Model")
axs[0].axis("tight")
axs[1].imshow(xinv_sm.T, vmin=0, vmax=1, cmap="gray")
axs[1].set_title("L2 Inversion")
axs[1].axis("tight")
axs[2].imshow(xinv.T, vmin=0, vmax=1, cmap="gray")
axs[2].set_title("TV-Reg Inversion")
axs[2].axis("tight")
fig.tight_layout()
Model, L2 Inversion, TV-Reg Inversion

Total running time of the script: ( 0 minutes 12.445 seconds)

Gallery generated by Sphinx-Gallery

17. Real/Complex Inversion

In this tutorial we will discuss two equivalent approaches to the solution of inverse problems with real-valued model vector and complex-valued data vector. In other words, we consider a modelling operator \(\mathbf{A}:\mathbb{F}^m \to \mathbb{C}^n\) (which could be the case for example for the real FFT).

Mathematically speaking, this problem can be solved equivalently by inverting the complex-valued problem:

\[\mathbf{y} = \mathbf{A} \mathbf{x}\]

or the real-valued augmented system

\[\begin{split}\DeclareMathOperator{\Real}{Re} \DeclareMathOperator{\Imag}{Im} \begin{bmatrix} \Real(\mathbf{y}) \\ \Imag(\mathbf{y}) \end{bmatrix} = \begin{bmatrix} \Real(\mathbf{A}) \\ \Imag(\mathbf{A}) \end{bmatrix} \mathbf{x}\end{split}\]

Whilst we already know how to solve the first problem, let’s see how we can solve the second one by taking advantage of the real method of the pylops.LinearOperator object. We will also wrap our linear operator into a pylops.MemoizeOperator which remembers the last N model and data vectors and by-passes the computation of the forward and/or adjoint pass whenever the same pair reappears. This is very useful in our case when we want to compute the real and the imag components of

import matplotlib.pyplot as plt
import numpy as np

import pylops

plt.close("all")
np.random.seed(0)

To start we create the forward problem

n = 5
x = np.arange(n) + 1.0

# make A
Ar = np.random.normal(0, 1, (n, n))
Ai = np.random.normal(0, 1, (n, n))
A = Ar + 1j * Ai
Aop = pylops.MatrixMult(A, dtype=np.complex128)
y = Aop @ x

Let’s check we can solve this problem using the first formulation

A1op = Aop.toreal(forw=False, adj=True)
xinv = A1op.div(y)

print(f"xinv={xinv}\n")
xinv=[1.+1.83186799e-15j 2.-8.88178420e-16j 3.-8.88178420e-16j
 4.-2.22044605e-16j 5.+6.66133815e-16j]

Let’s now see how we formulate the second problem

Amop = pylops.MemoizeOperator(Aop, max_neval=10)
Arop = Amop.toreal()
Aiop = Amop.toimag()

A1op = pylops.VStack([Arop, Aiop])
y1 = np.concatenate([np.real(y), np.imag(y)])
xinv1 = np.real(A1op.div(y1))

print(f"xinv1={xinv1}\n")
xinv1=[1. 2. 3. 4. 5.]

Total running time of the script: ( 0 minutes 0.010 seconds)

Gallery generated by Sphinx-Gallery

18. Deblending

The cocktail party problem arises when sounds from different sources mix before reaching our ears (or any recording device), requiring the brain (or any hardware in the recording device) to estimate individual sources from the received mixture. In seismic acquisition, an analog problem is present when multiple sources are fired simultaneously. This family of acquisition methods is usually referred to as simultaneous shooting and the problem of separating the blendend shot gathers into their individual components is called deblending. Whilst various firing strategies can be adopted, in this example we consider the continuos blending problem where a single source is fired sequentially at an interval shorter than the amount of time required for waves to travel into the Earth and come back.

Simply stated the forward problem can be written as:

\[\mathbf{d}^b = \boldsymbol\Phi \mathbf{d}\]

Here \(\mathbf{d} = [\mathbf{d}_1^T, \mathbf{d}_2^T,\ldots, \mathbf{d}_N^T]^T\) is a stack of \(N\) individual shot gathers, \(\boldsymbol\Phi=[\boldsymbol\Phi_1, \boldsymbol\Phi_2,\ldots, \boldsymbol\Phi_N]\) is the blending operator, \(\mathbf{d}^b\) is the so-called supergather than contains all shots superimposed to each other.

In order to successfully invert this severely underdetermined problem, two key ingredients must be introduced:

  • the firing time of each source (i.e., shifts of the blending operator) must be chosen to be dithered around a nominal regular, periodic firing interval. In our case, we consider shots of duration \(T=4s\), regular firing time of \(T_s=2s\) and a dithering code as follows \(\Delta t = U(-1,1)\);

  • prior information about the data to reconstruct, either in the form of regularization or preconditioning must be introduced. In our case we will use a patch-FK transform as preconditioner and solve the problem imposing sparsity in the transformed domain.

In other words, we aim to solve the following problem:

\[J = \|\mathbf{d}^b - \boldsymbol\Phi \mathbf{S}^H \mathbf{x}\|_2 + \epsilon \|\mathbf{x}\|_1\]

for which we will use the pylops.optimization.sparsity.fista solver.

import matplotlib.pyplot as plt
import numpy as np
from scipy.sparse.linalg import lobpcg as sp_lobpcg

import pylops

np.random.seed(10)
plt.close("all")

Let’s start by defining a blending operator

def Blending(nt, ns, dt, overlap, times, dtype="float64"):
    """Blending operator"""
    pad = int(overlap * nt)
    OpShiftPad = []
    for i in range(ns):
        PadOp = pylops.Pad(nt, (pad * i, pad * (ns - 1 - i)), dtype=dtype)
        ShiftOp = pylops.signalprocessing.Shift(
            pad * (ns - 1) + nt, times[i], axis=0, sampling=dt, real=False, dtype=dtype
        )
        OpShiftPad.append(ShiftOp * PadOp)
    return pylops.HStack(OpShiftPad)

We can now load and display a small portion of the MobilAVO dataset composed of 60 shots and a single receiver. This data is unblended.

data = np.load("../testdata/deblending/mobil.npy")
ns, nt = data.shape

dt = 0.004
t = np.arange(nt) * dt

fig, ax = plt.subplots(1, 1, figsize=(12, 8))
ax.imshow(
    data.T,
    cmap="gray",
    vmin=-50,
    vmax=50,
    extent=(0, ns, t[-1], 0),
    interpolation="none",
)
ax.set_title("CRG")
ax.set_xlabel("#Src")
ax.set_ylabel("t [s]")
ax.axis("tight")
plt.tight_layout()
CRG

We are now ready to define the blending operator, blend our data, and apply the adjoint of the blending operator to it. This is usually referred as pseudo-deblending: as we will see brings back each source to its own nominal firing time, but since sources partially overlap in time, it will also generate some burst like noise in the data. Deblending can hopefully fix this.

overlap = 0.5
pad = int(overlap * nt)
ignition_times = 2.0 * np.random.rand(ns) - 1.0
Bop = Blending(nt, ns, dt, overlap, ignition_times, dtype="complex128")
data_blended = Bop * data.ravel()
data_pseudo = Bop.H * data_blended.ravel()
data_pseudo = data_pseudo.reshape(ns, nt)

fig, ax = plt.subplots(1, 1, figsize=(12, 8))
ax.imshow(
    data_pseudo.T.real,
    cmap="gray",
    vmin=-50,
    vmax=50,
    extent=(0, ns, t[-1], 0),
    interpolation="none",
)
ax.set_title("Pseudo-deblended CRG")
ax.set_xlabel("#Src")
ax.set_ylabel("t [s]")
ax.axis("tight")
plt.tight_layout()
Pseudo-deblended CRG

We are finally ready to solve our deblending inverse problem

# Patched FK
dimsd = data.shape
nwin = (20, 80)
nover = (10, 40)
nop = (128, 128)
nop1 = (128, 65)
nwins = (5, 24)
dims = (nwins[0] * nop1[0], nwins[1] * nop1[1])

Fop = pylops.signalprocessing.FFT2D(nwin, nffts=nop, real=True)
Sop = pylops.signalprocessing.Patch2D(
    Fop.H, dims, dimsd, nwin, nover, nop1, tapertype="hanning"
)
# Overall operator
Op = Bop * Sop

# Compute max eigenvalue (we do this explicitly to be able to run this fast)
Op1 = pylops.LinearOperator(Op.H * Op, explicit=False)
X = np.random.rand(Op1.shape[0], 1).astype(Op1.dtype)
maxeig = sp_lobpcg(Op1, X=X, maxiter=5, tol=1e-10)[0][0]
alpha = 1.0 / maxeig

# Deblend
niter = 60
decay = (np.exp(-0.05 * np.arange(niter)) + 0.2) / 1.2

with pylops.disabled_ndarray_multiplication():
    p_inv = pylops.optimization.sparsity.fista(
        Op,
        data_blended.ravel(),
        niter=niter,
        eps=5e0,
        alpha=alpha,
        decay=decay,
        show=True,
    )[0]
data_inv = Sop * p_inv
data_inv = data_inv.reshape(ns, nt)

fig, axs = plt.subplots(1, 4, sharey=False, figsize=(12, 8))
axs[0].imshow(
    data.T.real,
    cmap="gray",
    extent=(0, ns, t[-1], 0),
    vmin=-50,
    vmax=50,
    interpolation="none",
)
axs[0].set_title("CRG")
axs[0].set_xlabel("#Src")
axs[0].set_ylabel("t [s]")
axs[0].axis("tight")
axs[1].imshow(
    data_pseudo.T.real,
    cmap="gray",
    extent=(0, ns, t[-1], 0),
    vmin=-50,
    vmax=50,
    interpolation="none",
)
axs[1].set_title("Pseudo-deblended CRG")
axs[1].set_xlabel("#Src")
axs[1].axis("tight")
axs[2].imshow(
    data_inv.T.real,
    cmap="gray",
    extent=(0, ns, t[-1], 0),
    vmin=-50,
    vmax=50,
    interpolation="none",
)
axs[2].set_xlabel("#Src")
axs[2].set_title("Deblended CRG")
axs[2].axis("tight")
axs[3].imshow(
    data.T.real - data_inv.T.real,
    cmap="gray",
    extent=(0, ns, t[-1], 0),
    vmin=-50,
    vmax=50,
    interpolation="none",
)
axs[3].set_xlabel("#Src")
axs[3].set_title("Blending error")
axs[3].axis("tight")
plt.tight_layout()
CRG, Pseudo-deblended CRG, Deblended CRG, Blending error
FISTA (soft thresholding)
--------------------------------------------------------------------------------
The Operator Op has 30500 rows and 998400 cols
eps = 5.000000e+00      tol = 1.000000e-10      niter = 60
alpha = 3.261681e-01    thresh = 8.154204e-01
--------------------------------------------------------------------------------
   Itn          x[0]              r2norm     r12norm     xupdate
     1   0.00e+00+0.00e+00j    3.463e+06   5.081e+06   1.251e+03
     2   0.00e+00+0.00e+00j    1.955e+06   4.282e+06   6.844e+02
     3   0.00e+00+0.00e+00j    1.134e+06   3.898e+06   5.635e+02
     4   0.00e+00+0.00e+00j    7.080e+05   3.681e+06   4.576e+02
     5   0.00e+00+0.00e+00j    4.873e+05   3.510e+06   3.831e+02
     6   0.00e+00+0.00e+00j    3.688e+05   3.348e+06   3.357e+02
     7   0.00e+00+0.00e+00j    3.014e+05   3.192e+06   3.048e+02
     8   0.00e+00+0.00e+00j    2.600e+05   3.048e+06   2.826e+02
     9   0.00e+00+0.00e+00j    2.316e+05   2.923e+06   2.642e+02
    10   0.00e+00+0.00e+00j    2.105e+05   2.816e+06   2.479e+02
    11   0.00e+00+0.00e+00j    1.934e+05   2.727e+06   2.327e+02
    21   0.00e+00+0.00e+00j    1.028e+05   2.435e+06   1.202e+02
    31   -0.00e+00+0.00e+00j    6.359e+04   2.450e+06   8.024e+01
    41   -0.00e+00+0.00e+00j    4.281e+04   2.490e+06   6.262e+01
    51   -0.00e+00+0.00e+00j    3.179e+04   2.520e+06   5.211e+01
    52   -0.00e+00+0.00e+00j    3.101e+04   2.523e+06   5.126e+01
    53   -0.00e+00+0.00e+00j    3.027e+04   2.525e+06   5.040e+01
    54   -0.00e+00+0.00e+00j    2.957e+04   2.527e+06   4.960e+01
    55   -0.00e+00+0.00e+00j    2.890e+04   2.529e+06   4.879e+01
    56   -0.00e+00+0.00e+00j    2.827e+04   2.531e+06   4.804e+01
    57   -0.00e+00+0.00e+00j    2.768e+04   2.533e+06   4.733e+01
    58   -0.00e+00+0.00e+00j    2.712e+04   2.535e+06   4.663e+01
    59   -0.00e+00+0.00e+00j    2.660e+04   2.536e+06   4.598e+01
    60   -0.00e+00+0.00e+00j    2.610e+04   2.538e+06   4.531e+01

Iterations = 60        Total time (s) = 44.24
--------------------------------------------------------------------------------

Finally, let’s look a bit more at what really happened under the hood. We display a number of patches and their associated FK spectrum

Sop1 = pylops.signalprocessing.Patch2D(
    Fop.H, dims, dimsd, nwin, nover, nop1, tapertype=None
)

# Original
p = Sop1.H * data.ravel()
preshape = p.reshape(nwins[0], nwins[1], nop1[0], nop1[1])

ix = 16
fig, axs = plt.subplots(2, 4, figsize=(12, 5))
fig.suptitle("Data patches")
for i in range(4):
    axs[0][i].imshow(np.fft.fftshift(np.abs(preshape[i, ix]).T, axes=1))
    axs[0][i].axis("tight")
    axs[1][i].imshow(
        np.real((Fop.H * preshape[i, ix].ravel()).reshape(nwin)).T,
        cmap="gray",
        vmin=-30,
        vmax=30,
        interpolation="none",
    )
    axs[1][i].axis("tight")
plt.tight_layout()

# Pseudo-deblended
p_pseudo = Sop1.H * data_pseudo.ravel()
p_pseudoreshape = p_pseudo.reshape(nwins[0], nwins[1], nop1[0], nop1[1])

ix = 16
fig, axs = plt.subplots(2, 4, figsize=(12, 5))
fig.suptitle("Pseudo-deblended patches")
for i in range(4):
    axs[0][i].imshow(np.fft.fftshift(np.abs(p_pseudoreshape[i, ix]).T, axes=1))
    axs[0][i].axis("tight")
    axs[1][i].imshow(
        np.real((Fop.H * p_pseudoreshape[i, ix].ravel()).reshape(nwin)).T,
        cmap="gray",
        vmin=-30,
        vmax=30,
        interpolation="none",
    )
    axs[1][i].axis("tight")
plt.tight_layout()

# Deblended
p_inv = Sop1.H * data_inv.ravel()
p_invreshape = p_inv.reshape(nwins[0], nwins[1], nop1[0], nop1[1])

ix = 16
fig, axs = plt.subplots(2, 4, figsize=(12, 5))
fig.suptitle("Deblended patches")
for i in range(4):
    axs[0][i].imshow(np.fft.fftshift(np.abs(p_invreshape[i, ix]).T, axes=1))
    axs[0][i].axis("tight")
    axs[1][i].imshow(
        np.real((Fop.H * p_invreshape[i, ix].ravel()).reshape(nwin)).T,
        cmap="gray",
        vmin=-30,
        vmax=30,
        interpolation="none",
    )
    axs[1][i].axis("tight")
plt.tight_layout()
  • Data patches
  • Pseudo-deblended patches
  • Deblended patches

Total running time of the script: ( 0 minutes 51.398 seconds)

Gallery generated by Sphinx-Gallery

19. Automatic Differentiation

This tutorial focuses on the use of pylops.TorchOperator to allow performing Automatic Differentiation (AD) on chains of operators which can be:

  • native PyTorch mathematical operations (e.g., torch.log, torch.sin, torch.tan, torch.pow, …)

  • neural network operators in torch.nn

  • PyLops linear operators

This opens up many opportunities, such as easily including linear regularization terms to nonlinear cost functions or using linear preconditioners with nonlinear modelling operators.

import matplotlib.pyplot as plt
import numpy as np
import torch
import torch.nn as nn
from torch.autograd import gradcheck

import pylops

plt.close("all")
np.random.seed(10)
torch.manual_seed(10)
<torch._C.Generator object at 0x7f686f92f450>

In this example we consider a simple multidimensional functional:

\[\mathbf{y} = \mathbf{A} sin(\mathbf{x})\]

and we use AD to compute the gradient with respect to the input vector evaluated at \(\mathbf{x}=\mathbf{x}_0\) : \(\mathbf{g} = d\mathbf{y} / d\mathbf{x} |_{\mathbf{x}=\mathbf{x}_0}\).

Let’s start by defining the Jacobian:

\[\begin{split}\textbf{J} = \begin{bmatrix} dy_1 / dx_1 & ... & dy_1 / dx_M \\ ... & ... & ... \\ dy_N / dx_1 & ... & dy_N / dx_M \end{bmatrix} = \begin{bmatrix} a_{11} cos(x_1) & ... & a_{1M} cos(x_M) \\ ... & ... & ... \\ a_{N1} cos(x_1) & ... & a_{NM} cos(x_M) \end{bmatrix} = \textbf{A} cos(\mathbf{x})\end{split}\]

Since both input and output are multidimensional, PyTorch backward actually computes the product between the transposed Jacobian and a vector \(\mathbf{v}\): \(\mathbf{g}=\mathbf{J^T} \mathbf{v}\).

To validate the correctness of the AD result, we can in this simple case also compute the Jacobian analytically and apply it to the same vector \(\mathbf{v}\) that we have provided to PyTorch backward.

nx, ny = 10, 6
x0 = torch.arange(nx, dtype=torch.double, requires_grad=True)

# Forward
A = np.random.normal(0.0, 1.0, (ny, nx))
At = torch.from_numpy(A)
Aop = pylops.TorchOperator(pylops.MatrixMult(A))
y = Aop.apply(torch.sin(x0))

# AD
v = torch.ones(ny, dtype=torch.double)
y.backward(v, retain_graph=True)
adgrad = x0.grad

# Analytical
J = At * torch.cos(x0)
anagrad = torch.matmul(J.T, v)

print("Input: ", x0)
print("AD gradient: ", adgrad)
print("Analytical gradient: ", anagrad)
Input:  tensor([0., 1., 2., 3., 4., 5., 6., 7., 8., 9.], dtype=torch.float64,
       requires_grad=True)
AD gradient:  tensor([ 0.1539, -0.2356,  1.4944, -3.1160, -2.1699,  0.4492,  1.6168,  1.9101,
        -0.4259,  1.7154], dtype=torch.float64)
Analytical gradient:  tensor([ 0.1539, -0.2356,  1.4944, -3.1160, -2.1699,  0.4492,  1.6168,  1.9101,
        -0.4259,  1.7154], dtype=torch.float64, grad_fn=<MvBackward>)

Similarly we can use the torch.autograd.gradcheck directly from PyTorch. Note that doubles must be used for this to succeed with very small eps and atol

input = (
    torch.arange(nx, dtype=torch.double, requires_grad=True),
    Aop.matvec,
    Aop.rmatvec,
    Aop.device,
    "cpu",
)
test = gradcheck(Aop.Top, input, eps=1e-6, atol=1e-4)
print(test)
True

Note that while matrix-vector multiplication could have been performed using the native PyTorch operator torch.matmul, in this case we have shown that we are also able to use a PyLops operator wrapped in pylops.TorchOperator. As already mentioned, this gives us the ability to use much more complex linear operators provided by PyLops within a chain of mixed linear and nonlinear AD-enabled operators. To conclude, let’s see how we can chain a torch convolutional network with PyLops pylops.Smoothing2D operator. First of all, we consider a single training sample.

class Network(nn.Module):
    def __init__(self, input_channels):
        super(Network, self).__init__()
        self.conv1 = nn.Conv2d(
            input_channels, input_channels // 2, kernel_size=3, padding=1
        )
        self.conv2 = nn.Conv2d(
            input_channels // 2, input_channels // 4, kernel_size=3, padding=1
        )
        self.activation = nn.LeakyReLU(0.2)
        self.maxpool = nn.MaxPool2d(kernel_size=2, stride=2)

    def forward(self, x):
        x = self.conv1(x)
        x = self.activation(x)
        x = self.conv2(x)
        x = self.activation(x)
        return x


net = Network(4)
Cop = pylops.TorchOperator(pylops.Smoothing2D((5, 5), dims=(32, 32)))

# Forward
x = torch.randn(1, 4, 32, 32).requires_grad_()
y = Cop.apply(net(x).view(-1)).reshape(32, 32)

# Backward
loss = y.sum()
loss.backward()

fig, axs = plt.subplots(1, 2, figsize=(12, 3))
axs[0].imshow(y.detach().numpy())
axs[0].set_title("Forward")
axs[0].axis("tight")
axs[1].imshow(x.grad.reshape(4 * 32, 32).T)
axs[1].set_title("Gradient")
axs[1].axis("tight")
plt.tight_layout()
Forward, Gradient

And finally we do the same with a batch of 3 training samples.

net = Network(4)
Cop = pylops.TorchOperator(pylops.Smoothing2D((5, 5), dims=(32, 32)), batch=True)

# Forward
x = torch.randn(3, 4, 32, 32).requires_grad_()
y = Cop.apply(net(x).reshape(3, 32 * 32)).reshape(3, 32, 32)

# Backward
loss = y.sum()
loss.backward()

fig, axs = plt.subplots(1, 2, figsize=(12, 3))
axs[0].imshow(y[0].detach().numpy())
axs[0].set_title("Forward")
axs[0].axis("tight")
axs[1].imshow(x.grad[0].reshape(4 * 32, 32).T)
axs[1].set_title("Gradient")
axs[1].axis("tight")
plt.tight_layout()
Forward, Gradient

Total running time of the script: ( 0 minutes 0.531 seconds)

Gallery generated by Sphinx-Gallery

Gallery generated by Sphinx-Gallery

Frequenty Asked Questions

1. Can I visualize my operator?

Yes, you can. Every operator has a method called todense that will return the dense matrix equivalent of the operator. Note, however, that in order to do so we need to allocate a numpy array of the size of your operator and apply the operator N times, where N is the number of columns of the operator. The allocation can be very heavy on your memory and the computation may take long time, so use it with care only for small toy examples to understand what your operator looks like. This method should however not be abused, as the reason of working with linear operators is indeed that you don’t really need to access the explicit matrix representation of an operator.

2. Can I have an older version of cupy or cusignal installed in my system ( cupy-cudaXX<8.1.0 or cusignal>=0.16.0 )?

Yes. Nevertheless you need to tell PyLops that you don’t want to use its cupy backend by setting the environment variable CUPY_PYLOPS=0 or CUPY_SIGNAL=0. Failing to do so will lead to an error when you import pylops because some of the cupyx routines that we use are not available in earlier version of cupy.

PyLops API

The Application Programming Interface (API) of PyLops can be loosely seen as composed of a stack of three main layers:

  • Linear operators: building blocks for the setting up of inverse problems

  • Solvers: interfaces to a variety of solvers, providing an easy way to augment an inverse problem with additional regularization and/or preconditioning term

  • Applications: high-level interfaces allowing users to easily setup and solve specific problems (while hiding the non-needed details - i.e., creation and setup of linear operators and solvers).

Linear operators

Templates

LinearOperator(*args, **kwargs)

Common interface for performing matrix-vector products.

FunctionOperator(*args, **kwargs)

Function Operator.

MemoizeOperator(*args, **kwargs)

Memoize Operator.

TorchOperator(*args, **kwargs)

Wrap a PyLops operator into a Torch function.

Basic operators

MatrixMult(*args, **kwargs)

Matrix multiplication.

Identity(*args, **kwargs)

Identity operator.

Zero(*args, **kwargs)

Zero operator.

Diagonal(*args, **kwargs)

Diagonal operator.

Transpose(*args, **kwargs)

Transpose operator.

Flip(*args, **kwargs)

Flip along an axis.

Roll(*args, **kwargs)

Roll along an axis.

Pad(*args, **kwargs)

Pad operator.

Sum(*args, **kwargs)

Sum operator.

Symmetrize(*args, **kwargs)

Symmetrize along an axis.

Restriction(*args, **kwargs)

Restriction (or sampling) operator.

Regression(*args, **kwargs)

Polynomial regression.

LinearRegression(taxis[, dtype])

Linear regression.

CausalIntegration(*args, **kwargs)

Causal integration.

Spread(*args, **kwargs)

Spread operator.

VStack(*args, **kwargs)

Vertical stacking.

HStack(*args, **kwargs)

Horizontal stacking.

Block(ops[, nproc, dtype])

Block operator.

BlockDiag(*args, **kwargs)

Block-diagonal operator.

Kronecker(*args, **kwargs)

Kronecker operator.

Real(*args, **kwargs)

Real operator.

Imag(*args, **kwargs)

Imag operator.

Conj(*args, **kwargs)

Complex conjugate operator.

Smoothing and derivatives

Smoothing1D(nsmooth, dims[, axis, dtype])

1D Smoothing.

Smoothing2D(nsmooth, dims[, axes, dtype])

2D Smoothing.

FirstDerivative(*args, **kwargs)

First derivative.

SecondDerivative(*args, **kwargs)

Second derivative.

Laplacian(dims[, axes, weights, sampling, ...])

Laplacian.

Gradient(dims[, sampling, edge, kind, dtype])

Gradient.

FirstDirectionalDerivative(dims, v[, ...])

First Directional derivative.

SecondDirectionalDerivative(dims, v[, ...])

Second Directional derivative.

Signal processing

Convolve1D(*args, **kwargs)

1D convolution operator.

Convolve2D(dims, h[, offset, axes, method, ...])

2D convolution operator.

ConvolveND(*args, **kwargs)

ND convolution operator.

Interp(dims, iava[, axis, kind, dtype, name])

Interpolation operator.

Bilinear(*args, **kwargs)

Bilinear interpolation operator.

FFT(dims[, axis, nfft, sampling, norm, ...])

One dimensional Fast-Fourier Transform.

FFT2D(dims[, axes, nffts, sampling, norm, ...])

Two dimensional Fast-Fourier Transform.

FFTND(dims[, axes, nffts, sampling, norm, ...])

N-dimensional Fast-Fourier Transform.

Shift(dims, shift[, axis, nfft, sampling, ...])

Shift operator

DWT(*args, **kwargs)

One dimensional Wavelet operator.

DWT2D(*args, **kwargs)

Two dimensional Wavelet operator.

Seislet(*args, **kwargs)

Two dimensional Seislet operator.

Radon2D(taxis, haxis, pxaxis[, kind, ...])

Two dimensional Radon transform.

Radon3D(taxis, hyaxis, hxaxis, pyaxis, pxaxis)

Three dimensional Radon transform.

ChirpRadon2D(*args, **kwargs)

2D Chirp Radon transform

ChirpRadon3D(*args, **kwargs)

3D Chirp Radon transform

Sliding1D(Op, dim, dimd, nwin, nover[, ...])

1D Sliding transform operator.

Sliding2D(Op, dims, dimsd, nwin, nover[, ...])

2D Sliding transform operator.

Sliding3D(Op, dims, dimsd, nwin, nover, nop)

3D Sliding transform operator.w

Patch2D(Op, dims, dimsd, nwin, nover, nop[, ...])

2D Patch transform operator.

Patch3D(Op, dims, dimsd, nwin, nover, nop[, ...])

3D Patch transform operator.

Fredholm1(*args, **kwargs)

Fredholm integral of first kind.

Wave-Equation processing

PressureToVelocity(nt, nr, dt, dr, rho, vel)

Pressure to Vertical velocity conversion.

UpDownComposition2D(nt, nr, dt, dr, rho, vel)

2D Up-down wavefield composition.

UpDownComposition3D(nt, nr, dt, dr, rho, vel)

3D Up-down wavefield composition.

MDC(G, nt, nv[, dt, dr, twosided, ...])

Multi-dimensional convolution.

PhaseShift(vel, dz, nt, freq, kx[, ky, ...])

Phase shift operator

Kirchhoff(*args, **kwargs)

Kirchhoff Demigration operator.

AcousticWave2D(*args, **kwargs)

Devito Acoustic propagator.

Geophysical subsurface characterization

avo.AVOLinearModelling(*args, **kwargs)

AVO Linearized modelling.

poststack.PoststackLinearModelling(wav, nt0)

Post-stack linearized seismic modelling operator.

prestack.PrestackLinearModelling(wav, theta)

Pre-stack linearized seismic modelling operator.

prestack.PrestackWaveletModelling(m, theta, nwav)

Pre-stack linearized seismic modelling operator for wavelet.

Solvers

Template

Solver(Op[, callbacks])

This is a template class which a user must subclass when implementing a new solver.

Basic

CG(Op[, callbacks])

Conjugate gradient

CGLS(Op[, callbacks])

Conjugate gradient least squares

LSQR(Op)

Solve an overdetermined system of equations given an operator Op and data y using LSQR iterations.

cg(Op, y[, x0, niter, tol, show, itershow, ...])

Conjugate gradient

cgls(Op, y[, x0, niter, damp, tol, show, ...])

Conjugate gradient least squares

lsqr(Op, y[, x0, damp, atol, btol, conlim, ...])

LSQR

Least-squares

NormalEquationsInversion(Op[, callbacks])

Inversion of normal equations.

RegularizedInversion(Op[, callbacks])

Regularized inversion.

PreconditionedInversion(Op[, callbacks])

Preconditioned inversion.

normal_equations_inversion(Op, y, Regs[, ...])

Inversion of normal equations.

regularized_inversion(Op, y, Regs[, x0, ...])

Regularized inversion.

preconditioned_inversion(Op, y, P[, x0, ...])

Preconditioned inversion.

Sparsity

IRLS(Op[, callbacks])

Iteratively reweighted least squares.

OMP(Op[, callbacks])

Orthogonal Matching Pursuit (OMP).

ISTA(Op[, callbacks])

Iterative Shrinkage-Thresholding Algorithm (ISTA).

FISTA(Op[, callbacks])

Fast Iterative Shrinkage-Thresholding Algorithm (FISTA).

SPGL1(Op[, callbacks])

Spectral Projected-Gradient for L1 norm.

SplitBregman(Op[, callbacks])

Split Bregman for mixed L2-L1 norms.

irls(Op, y[, x0, nouter, threshR, epsR, ...])

Iteratively reweighted least squares.

omp(Op, y[, niter_outer, niter_inner, ...])

Orthogonal Matching Pursuit (OMP).

ista(Op, y[, x0, niter, SOp, eps, alpha, ...])

Iterative Shrinkage-Thresholding Algorithm (ISTA).

fista(Op, y[, x0, niter, SOp, eps, alpha, ...])

Fast Iterative Shrinkage-Thresholding Algorithm (FISTA).

spgl1(Op, y[, x0, SOp, tau, sigma, show])

Spectral Projected-Gradient for L1 norm.

splitbregman(Op, y, RegsL1[, x0, ...])

Split Bregman for mixed L2-L1 norms.

Callbacks

Callbacks()

This is a template class which a user must subclass when implementing callbacks for a solver.

MetricsCallback(xtrue[, Op, which])

Metrics callback

Applications

Wave-Equation processing

SeismicInterpolation(data, nrec, iava[, ...])

Seismic interpolation (or regularization).

Deghosting(p, nt, nr, dt, dr, vel, zrec[, ...])

Wavefield deghosting.

WavefieldDecomposition(p, vz, nt, nr, dt, ...)

Up-down wavefield decomposition.

MDD(G, d[, dt, dr, nfmax, wav, twosided, ...])

Multi-dimensional deconvolution.

Marchenko(R[, dt, nt, dr, nfmax, wav, toff, ...])

Marchenko redatuming

LSM(z, x, t, srcs, recs, vel, wav, wavcenter)

Least-squares Migration (LSM).

Geophysical subsurface characterization

poststack.PoststackInversion(data, wav[, ...])

Post-stack linearized seismic inversion.

prestack.PrestackInversion(data, theta, wav)

Pre-stack linearized seismic inversion.

PyLops Utilities

Alongside with its Linear Operators and Solvers, PyLops contains also a number of auxiliary routines performing universal tasks that are used by several operators or simply within one or more Tutorials for the preparation of input data and subsequent visualization of results.

Dot-test

dottest(Op[, nr, nc, rtol, atol, ...])

Dot test.

Decorators

add_ndarray_support_to_solver(func)

Decorator which converts a solver-type function that only supports a 1d-array into one that supports one (dimsd-shaped) ndarray.

disable_ndarray_multiplication(func)

Decorator which disables ndarray multiplication.

reshaped([func, forward, swapaxis])

Decorator for the common reshape/flatten pattern used in many operators.

Describe

describe(Op)

Describe a PyLops operator

Estimators

trace_hutchinson(Op[, neval, batch_size, ...])

Trace of linear operator using the Hutchinson method.

trace_hutchpp(Op[, neval, sampler, backend])

Trace of linear operator using the Hutch++ method.

trace_nahutchpp(Op[, neval, sampler, c1, ...])

Trace of linear operator using the NA-Hutch++ method.

Metrics

mae(xref, xcmp)

Mean Absolute Error (MAE)

mse(xref, xcmp)

Mean Square Error (MSE)

snr(xref, xcmp)

Signal to Noise Ratio (SNR)

psnr(xref, xcmp[, xmax])

Peak Signal to Noise Ratio (PSNR)

Geophysical Reservoir characterization

avo.zoeppritz_scattering(vp1, vs1, rho1, ...)

Zoeppritz solution.

avo.zoeppritz_element(vp1, vs1, rho1, vp0, ...)

Single element of Zoeppritz solution.

avo.zoeppritz_pp(vp1, vs1, rho1, vp0, vs0, ...)

PP reflection coefficient from the Zoeppritz scattering matrix.

avo.approx_zoeppritz_pp(vp1, vs1, rho1, vp0, ...)

PP reflection coefficient from the approximate Zoeppritz equation.

avo.akirichards(theta, vsvp[, n])

Three terms Aki-Richards approximation.

avo.fatti(theta, vsvp[, n])

Three terms Fatti approximation.

avo.ps(theta, vsvp[, n])

PS reflection coefficient

Scalability test

scalability_test(Op, x[, workers, forward])

Scalability test.

Sliding and Patching

sliding1d.sliding1d_design(dimd, nwin, ...)

Design Sliding1D operator

sliding2d.sliding2d_design(dimsd, nwin, ...)

Design Sliding2D operator

sliding3d.sliding3d_design(dimsd, nwin, ...)

Design Sliding3D operator

patch2d.patch2d_design(dimsd, nwin, nover, nop)

Design Patch2D operator

patch3d.patch3d_design(dimsd, nwin, nover, nop)

Design Patch3D operator

Synthetics

seismicevents.makeaxis(par)

Create axes t, x, and y axes

seismicevents.linear2d(x, t, v, t0, theta, ...)

Linear 2D events

seismicevents.parabolic2d(x, t, t0, px, pxx, ...)

Parabolic 2D events

seismicevents.hyperbolic2d(x, t, t0, vrms, ...)

Hyperbolic 2D events

seismicevents.linear3d(x, y, t, v, t0, ...)

Linear 3D events

seismicevents.hyperbolic3d(x, y, t, t0, ...)

Hyperbolic 3D events

marchenko.directwave(wav, trav, nt, dt[, ...])

Analytical direct wave in acoustic media

Signal-processing

signalprocessing.convmtx(h, n)

Convolution matrix

signalprocessing.nonstationary_convmtx(H, n)

Convolution matrix from a bank of filters

signalprocessing.dip_estimate(d[, dz, dx, ...])

Local dip estimation

signalprocessing.slope_estimate(d[, dz, dx, ...])

Local slope estimation

Tapers

tapers.taper2d(nt, nmask, ntap[, tapertype])

2D taper

tapers.taper3d(nt, nmask, ntap[, tapertype])

3D taper

tapers.tapernd(nmask, ntap[, tapertype])

ND taper

Wavelets

wavelets.gaussian(t[, std])

Gaussian wavelet

wavelets.klauder(t[, f, taper])

Klauder wavelet

wavelets.ormsby(t[, f, taper])

Ormsby wavelet

wavelets.ricker(t[, f0, taper])

Ricker wavelet

Implementing new operators

Users are welcome to create new operators and add them to the PyLops library.

In this tutorial, we will go through the key steps in the definition of an operator, using the pylops.Diagonal as an example. This is a very simple operator that applies a diagonal matrix to the model in forward mode and to the data in adjoint mode.

Creating the operator

The first thing we need to do is to create a new file with the name of the operator we would like to implement. Note that as the operator will be a class, we need to follow the UpperCaseCamelCase convention both for the class itself and for the filename.

Once we have created the file, we will start by importing the modules that will be needed by the operator. While this varies from operator to operator, you will always need to import the pylops.LinearOperator class, which will be used as parent class for any of our operators:

from pylops import LinearOperator

This class is a child of the scipy.sparse.linalg.LinearOperator class itself which implements the same methods of its parent class as well as an additional method for quick inversion: such method can be easily accessed by using \ between the operator and the data (e.g., A \ y).

After that we define our new object:

class Diagonal(LinearOperator):

followed by a numpydoc docstring (starting with r""" and ending with """) containing the documentation of the operator. Such docstring should contain at least a short description of the operator, a Parameters section with a detailed description of the input parameters and a Notes section providing a mathematical explanation of the operator. Take a look at some of the core operators of PyLops to get a feeling of the level of details of the mathematical explanation.

Initialization (__init__)

We then need to create the __init__ where the input parameters are passed and saved as members of our class. While the input parameters change from operator to operator, it is always required to create three members, the first called shape with a tuple containing the dimensions of the operator in the data and model space, the second called dtype with the data type object (np.dtype) of the model and data, and the third called explicit with a boolean (True or False) identifying if the operator can be inverted by a direct solver or requires an iterative solver. This member is True if the operator has also a member A that contains the matrix to be inverted like for example in the pylops.MatrixMult operator, and it will be False otherwise. In this case we have another member called d which is equal to the input vector containing the diagonal elements of the matrix we want to multiply to the model and data.

def __init__(self, d, dtype=None):
    self.d = d.ravel()
    self.shape = (len(self.d), len(self.d))
    self.dtype = np.dtype(dtype)
    self.explicit = False

Alternatively, since version 2.0.0, the recommended way of initializing operators derived from the base pylops.LinearOperator class is to invoke super to assign the required attributes:

def __init__(self, d, dtype=None):
    self.d = d.ravel()
    super().__init__(dtype=np.dtype(dtype), shape=(len(self.d), len(self.d)))

In this case, there is no need to declare explicit as it already defaults to False. Since version 2.0.0, every pylops.LinearOperator class is imbued with dims, dimsd, clinear and explicit, in addition to the required dtype and shape. See the docs of pylops.LinearOperator for more information about what these attributes mean.

Forward mode (_matvec)

We can then move onto writing the forward mode in the method _matvec. In other words, we will need to write the piece of code that will implement the following operation \(\mathbf{y} = \mathbf{A}\mathbf{x}\). Such method is always composed of two inputs (the object itself self and the input model x). In our case the code to be added to the forward is very simple, we will just need to apply element-wise multiplication between the model \(\mathbf{x}\) and the elements along the diagonal contained in the array \(\mathbf{d}\). We will finally need to return the result of this operation:

def _matvec(self, x):
    return self.d * x
Adjoint mode (_rmatvec)

Finally we need to implement the adjoint mode in the method _rmatvec. In other words, we will need to write the piece of code that will implement the following operation \(\mathbf{x} = \mathbf{A}^H\mathbf{y}\). Such method is also composed of two inputs (the object itself self and the input data y). In our case the code to be added to the forward is the same as the one from the forward (but this will obviously be different from operator to operator):

def _rmatvec(self, x):
    return self.d * x

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 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.

We can generate random vectors \(\mathbf{u}\) and \(\mathbf{v}\) and verify the the following equality within a numerical tolerance:

\[(\mathbf{A}*\mathbf{u})^H*\mathbf{v} = \mathbf{u}^H*(\mathbf{A}^H*\mathbf{v})\]

The method pylops.utils.dottest implements such a test for you, all you need to do is create a new test within an existing test_*.py file in the pytests folder (or in a new file).

Generally a test file will start with a number of dictionaries containing different parameters we would like to use in the testing of one or more operators. The test itself starts with a decorator that contains a list of all (or some) of dictionaries that will would like to use for our specific operator, followed by the definition of the test

@pytest.mark.parametrize("par", [(par1),(par2)])
def test_Diagonal(par):

At this point we can first of all create the operator and run the pylops.utils.dottest preceded by the assert command. Moreover, the forward and adjoint methods should tested towards expected outputs or even better, when the operator allows it (i.e., operator is invertible), a small inversion should be run and the inverted model tested towards the input model.

"""Dot-test and inversion for diagonal operator
"""
d = np.arange(par['nx']) + 1.

Dop = Diagonal(d)
assert dottest(Dop, par['nx'], par['nx'],
               complexflag=0 if par['imag'] == 1 else 3)

x = np.ones(par['nx'])
xlsqr = lsqr(Dop, Dop * x, damp=1e-20, iter_lim=300, show=0)[0]
assert_array_almost_equal(x, xlsqr, decimal=4)

Documenting the operator

Once the operator has been created, we can add it to the documentation of PyLops. To do so, simply add the name of the operator within the index.rst file in docs/source/api directory.

Moreover, in order to facilitate the user of your operator by other users, a simple example should be provided as part of the Sphinx-gallery of the documentation of the PyLops library. The directory examples contains several scripts that can be used as template.

Final checklist

Before submitting your new operator for review, use the following checklist to ensure that your code adheres to the guidelines of PyLops:

  • you have created a new file containing a single class (or a function when the new operator is a simple combination of existing operators - see pylops.Laplacian for an example of such operator) and added to a new or existing directory within the pylops package.

  • the new class contains at least __init__, _matvec and _matvec methods.

  • the new class (or function) has a numpydoc docstring documenting at least the input Parameters and with a Notes section providing a mathematical explanation of the operator

  • a new test has been added to an existing test_*.py file within the pytests folder. The test should verify that the new operator passes the pylops.utils.dottest. Moreover it is advisable to create a small toy example where the operator is applied in forward mode and the resulting data is inverted using \ from pylops.LinearOperator.

  • the new operator is used within at least one example (in examples directory) or one tutorial (in tutorials directory).

Implementing new solvers

Users are welcome to create new solvers and add them to the PyLops library.

In this tutorial, we will go through the key steps in the definition of a solver, using the pylops.optimization.basic.CG as an example.

Note

In case the solver that you are planning to create falls within the category of proximal solvers, we encourage to consider adding it to the PyProximal project.

Creating the solver

The first thing we need to do is to locate a file containing solvers in the same family of the solver we plan to include, or create a new file with the name of the solver we would like to implement (or preferably its family). Note that as the solver will be a class, we need to follow the UpperCaseCamelCase convention for the class itself but not for the filename.

At this point we can start by importing the modules that will be needed by the solver. This varies from solver to solver, however you will always need to import the pylops.optimization.basesolver.Solver which will be used as parent class for any of our solvers. Moreover, we always recommend to import pylops.utils.backend.get_array_module as solvers should be written in such a way that it can work both with numpy and cupy arrays. See later for details.

import time

import numpy as np

from pylops.optimization.basesolver import Solver
from pylops.utils.backend import get_array_module

After that we define our new object:

class CG(Solver):

followed by a numpydoc docstring (starting with r""" and ending with """) containing the documentation of the solver. Such docstring should contain at least a short description of the solver, a Parameters section with a description of the input parameters of the associated _init__ method and a Notes section providing a reference to the original solver and possibly a concise mathematical explanation of the solver. Take a look at some of the core solver of PyLops to get a feeling of the level of details of the mathematical explanation.

As for any Python class, our solver will need an __init__ method. In this case, however, we will just rely on that of the base class. Two input parameters are passed to the __init__ method and saved as members of our class, namely the operator \(\mathbf{Op}\) associated with the system of equations we wish to solve, \(\mathbf{y}=\mathbf{Op}\,\mathbf{x}\), and optionally a pylops.optimization.callback.Callbacks object. Moreover, an additional parameters is created that contains the current time (this is used later to report the execution time of the solver). Here is the __init__ method of the base class:

def __init__(self, Op, callbacks=None):
    self.Op = Op
    self.callbacks = callbacks
    self._registercallbacks()
    self.tstart = time.time()

We can now move onto writing the setup of the solver in the method setup. We will need to write a piece of code that prepares the solver prior to being able to apply a step. In general, this requires defining the data vector y and the initial guess of the solver x0 (if not provided, this will be automatically set to be a zero vector), alongside various hyperparameters of the solver — e.g., those involved in the stopping criterion. For example in this case we only have two parameters: niter refers to the maximum allowed number of iterations, and tol to tolerance on the residual norm (the solver will be stopped if this is smaller than the chosen tolerance). Moreover, we always have the possibility to decide whether we want to operate the solver (in this case its setup part) in verbose or silent mode. This is driven by the show parameter. We will soon discuss how to choose what to print on screen in case of verbose mode (show=True).

The setup method can be loosely seen as composed of three parts. First, the data vector and hyperparameters are stored as members of the class. Moreover the type of the y vector is checked to evaluate whether to use numpy or cupy for algebraic operations (this is done by self.ncp = get_array_module(y)). Second, the starting guess is initialized using either the provided vector x0 or a zero vector. Finally, a number of variables are initialized to be used inside the step method to keep track of the optimization process. Moreover, note that the setup method returns the created starting guess x (does not store it as member of the class).

def setup(self, y, x0=None, niter=None, tol=1e-4, show=False):

    self.y = y
    self.tol = tol
    self.niter = niter
    self.ncp = get_array_module(y)

    # initialize solver
    if x0 is None:
        x = self.ncp.zeros(self.Op.shape[1], dtype=self.y.dtype)
        self.r = self.y.copy()
    else:
        x = x0.copy()
        self.r = self.y - self.Op.matvec(x)
    self.c = self.r.copy()
    self.kold = self.ncp.abs(self.r.dot(self.r.conj()))

    # create variables to track the residual norm and iterations
    self.cost = []
    self.cost.append(np.sqrt(self.kold))
    self.iiter = 0

    # print setup
    if show:
        self._print_setup(np.iscomplexobj(x))
    return x

At this point, we need to implement the core of the solver, the step method. Here, we take the input at the previous iterate, update it following the rule of the solver of choice, and return it. The other input parameter required by this method is show to choose whether we want to print a report of the step on screen or not. However, if appropriate, a user can add additional input parameters. For CG, the step is:

def step(self, x, show=False):
    Opc = self.Op.matvec(self.c)
    cOpc = self.ncp.abs(self.c.dot(Opc.conj()))
    a = self.kold / cOpc
    x += a * self.c
    self.r -= a * Opc
    k = self.ncp.abs(self.r.dot(self.r.conj()))
    b = k / self.kold
    self.c = self.r + b * self.c
    self.kold = k
    self.iiter += 1
    self.cost.append(np.sqrt(self.kold))
    if show:
        self._print_step(x)
    return x

Similarly, we also implement a run method that is in charge of running a number of iterations by repeatedly calling the step method. It is also usually convenient to implement a finalize method; this method can do any required post-processing that should not be applied at the end of each step, rather at the end of the entire optimization process. For CG, this is as simple as converting the cost variable from a list to a numpy array. For more details, see our implementations for CG.

Last but not least, we can wrap it all up in the solve method. This method takes as input the data, the initial model and the same hyperparameters of the setup method and runs the entire optimization process. For CG:

def solve(self, y, x0=None, niter=10, tol=1e-4, show=False, itershow=[10, 10, 10]):
    x = self.setup(y=y, x0=x0, niter=niter, tol=tol, show=show)
    x = self.run(x, niter, show=show, itershow=itershow)
    self.finalize(show)
    return x, self.iiter, self.cost

And that’s it, we have implemented our first solver operator!

Although the methods that we just described are enough to implement any solver of choice, we find important to provide users with feedback during the inversion process. Imagine that the modelling operator is very expensive and can take minutes (or even hours to run), we don’t want to leave a user waiting for hours before they can tell if the solver has done something meaningful. To avoid such scenario, we can implement so called _print_* methods where *=solver, setup, step, finalize that print on screen some useful information (e.g., first value of the current estimate, norm of residual, etc.). The solver and finalize print are already implemented in the base class, the other two must be implemented when creating a new solver. When these methods are implemented and a user passes show=True to the associated method, our solver will provide such information on screen throughout the inverse process. To better understand how to write such methods, we suggest to look into the source code of the CG method.

Finally, to be backward compatible with versions of PyLops <v2.0.0, we also want to create a function with the same name of the class-based solver (but in small letters) which simply instantiates the solver and runs it. This function is usually placed in the same file of the class-based solver and snake_case should be used for its name. This function generally takes all the mandatory and optional parameters of the solver as input and returns some of the most valuable properties of the class-based solver object. An example for CG is:

def cg(Op, y, x0, niter=10, tol=1e-4, show=False, itershow=[10, 10, 10], callback=None):
    cgsolve = CG(Op)
    if callback is not None:
        cgsolve.callback = callback
    x, iiter, cost = cgsolve.solve(
        y=y, x0=x0, tol=tol, niter=niter, show=show, itershow=itershow
    )
    return x, iiter, cost

Testing the solver

Being able to write a solver is not yet a guarantee of the fact that the solver is correct, or in other words that the solver can converge to a correct solution (at least in the case of full rank operator).

We encourage to create a new test within an existing test_*.py file in the pytests folder (or in a new file). We also encourage to test the function-bases solver, as this will implicitly test the underlying class-based solver.

Generally a test file will start with a number of dictionaries containing different parameters we would like to use in the testing of one or more solvers. The test itself starts with a decorator that contains a list of all (or some) of dictionaries that will would like to use for our specific operator, followed by the definition of the test:

@pytest.mark.parametrize("par", [(par1),(par2)])
def test_CG(par):

At this point we can first create a full-rank operator, an input vector and compute the associated data. We can then run the solver for a certain number of iterations, checking that the solution agrees with the true x within a certain tolerance:

"""CG with linear operator
"""
np.random.seed(10)

A = np.random.normal(0, 10, (par["ny"], par["nx"]))
A = np.conj(A).T @ A  # to ensure definite positive matrix
Aop = MatrixMult(A, dtype=par["dtype"])

x = np.ones(par["nx"])
x0 = np.random.normal(0, 10, par["nx"])

y = Aop * x
xinv = cg(Aop, y, x0=x0, niter=par["nx"], tol=1e-5, show=True)[0]
assert_array_almost_equal(x, xinv, decimal=4)

Documenting the solver

Once the solver has been created, we can add it to the documentation of PyLops. To do so, simply add the name of the operator within the index.rst file in docs/source/api directory.

Moreover, in order to facilitate the user of your operator by other users, a simple example should be provided as part of the Sphinx-gallery of the documentation of the PyLops library. The directory examples contains several scripts that can be used as template.

Final checklist

Before submitting your new solver for review, use the following checklist to ensure that your code adheres to the guidelines of PyLops:

  • you have added the new solver to a new or existing file in the optimization directory within the pylops package.

  • the new class contains at least __init__, setup, step, run, finalize, and solve methods.

  • each of the above methods have a numpydoc docstring documenting at least the input Parameters and the __init__ method contains also a Notes section providing a mathematical explanation of the solver.

  • a new test has been added to an existing test_*.py file within the pytests folder. The test should verify that the new solver converges to the true solution for a well-designed inverse problem (i.e., full rank operator).

  • the new solver is used within at least one example (in examples directory) or one tutorial (in tutorials directory).

Contributing

Contributions are welcome and greatly appreciated!

The best way to get in touch with the core developers and maintainers is to join the PyLops slack channel as well as open new Issues directly from the GitHub repo.

Moreover, take a look at the Roadmap page for a list of current ideas for improvements and additions to the PyLops library.

Welcomed contributions

Bug reports

Report bugs at https://github.com/PyLops/pylops/issues.

If you are playing with the PyLops library and find a bug, please report it including:

  • Your operating system name and version.

  • Any details about your Python environment.

  • Detailed steps to reproduce the bug.

New operators and features

Open an issue at https://github.com/PyLops/pylops/issues with tag enhancement.

If you are proposing a new operator or a new feature:

  • Explain in detail how it should work.

  • Keep the scope as narrow as possible, to make it easier to implement.

Fix issues

There is always a backlog of issues that need to be dealt with. Look through the GitHub Issues for operator/feature requests or bugfixes.

Add examples or improve documentation

Writing new operators is not the only way to get involved and contribute. Create examples with existing operators as well as improving the documentation of existing operators is as important as making new operators and very much encouraged.

Step-by-step instructions for contributing

Ready to contribute?

  1. Follow all instructions in Step-by-step installation for developers.

  2. Create a branch for local development, usually starting from the dev branch:

>> git checkout -b name-of-your-branch dev

Now you can make your changes locally.

3. When you’re done making changes, check that your code follows the guidelines for Implementing new operators and that the both old and new tests pass successfully:

>> make tests
  1. Run flake8 to check the quality of your code:

>> make lint

Note that PyLops does not enforce full compliance with flake8, rather this is used as a guideline and will also be run as part of our CI. Make sure to limit to a minimum flake8 warnings before making a PR.

  1. Update the docs

>> make docupdate
  1. Commit your changes and push your branch to GitHub:

>> git add .
>> git commit -m "Your detailed description of your changes."
>> git push origin name-of-your-branch

Remember to add -u when pushing the branch for the first time. We recommend using Conventional Commits to format your commit messages, but this is not enforced.

  1. Submit a pull request through the GitHub website.

Pull Request Guidelines

Before you submit a pull request, check that it meets these guidelines:

  1. The pull request should include new tests for all the core routines that have been developed.

  2. If the pull request adds functionality, the docs should be updated accordingly.

  3. Ensure that the updated code passes all tests.

Project structure

This repository is organized as follows: * pylops: Python library containing various linear operators and auxiliary routines * pytests: set of pytests * testdata: sample datasets used in pytests and documentation * docs: Sphinx documentation * examples: set of python script examples for each linear operator to be embedded in documentation using sphinx-gallery * tutorials: set of python script tutorials to be embedded in documentation using sphinx-gallery

Changelog

Version 2.0.0

Released on: 12/08/2022

PyLops has undergone significant changes in this release, including new LinearOperator s, more features, new examples and bugfixes. To aid users in navigating the breaking changes, we provide the following document MIGRATION_V1_V2.md.

New Features

Documentation

Version 1.18.3

Released on: 30/07/2022

  • Refractored pylops.utils.dottest, and added two new optional input parameters (atol and rtol)

  • Added optional parameter densesolver to pylops.LinearOperator.div

  • Fixed pylops.optimization.basic.lsqr, pylops.optimization.sparsity.ISTA, and pylops.optimization.sparsity.FISTA to work with cupy arrays. This change was required by how recent cupy versions handle scalars, which are not converted directly into float types, rather kept as cupy arrays.

  • Fix bug in pylops.waveeqprocessing.Deghosting introduced in commit 7e596d4.

Version 1.18.2

Released on: 29/04/2022

  • Refractored pylops.utils.dottest, and added two new optional input parameters (atol and rtol)

  • Added optional parameter densesolver to pylops.LinearOperator.div

Version 1.18.1

Released on: 29/04/2022

  • !DELETED! due to a mistake in the release process

Version 1.18.0

Released on: 19/02/2022

  • Added NMO example to gallery

  • Extended pylops.Laplacian to N-dimensional arrays

  • Added forward kind to pylops.SecondDerivative and pylops.Laplacian

  • Added chirp-sliding kind to pylops.waveeqprocessing.seismicinterpolation.SeismicInterpolation

  • Fixed bug due to the new internal structure of LinearOperator submodule introduced in scipy1.8.0

Version 1.17.0

Released on: 29/01/2022

Version 1.16.0

Released on: 11/12/2021

  • Added pylops.utils.estimators submodule for trace estimation

  • Added x0 in pylops.optimization.sparsity.ISTA and pylops.optimization.sparsity.FISTA to handle non-zero initial guess

  • Modified pylops.optimization.sparsity.ISTA and pylops.optimization.sparsity.FISTA to handle multiple right hand sides

  • Modified creation of haxis in pylops.signalprocessing.Radon2D and pylops.signalprocessing.Radon3D to allow for uncentered spatial axes

  • Fixed _rmatvec for explicit in pylops.LinearOperator._ColumnLinearOperator

Version 1.15.0

Released on: 23/10/2021

Version 1.14.0

Released on: 09/07/2021

  • Added pylops.optimization.solver.lsqr solver

  • Added utility routine pylops.utils.scalability_test for scalability tests when using multiprocessing

  • Added pylops.avo.avo.ps AVO modelling option and restructured pylops.avo.prestack.PrestackLinearModelling to allow passing any function handle that can perform AVO modelling apart from those directly available

  • Added R-linear operators (when setting the property clinear=False of a linear operator). pylops.basicoperators.Real, pylops.basicoperators.Imag, and pylops.basicoperators.Conj

  • Added possibility to run operators pylops.basicoperators.HStack, pylops.basicoperators.VStack, pylops.basicoperators.Block pylops.basicoperators.BlockDiag, and pylops.signalprocessing.Sliding3D using multiprocessing

  • Added dtype to vector X when using scipy.sparse.linalg.lobpcg in eigs method of pylops.LinearOperator

  • Use kind=forward fot FirstDerivative in pylops.avo.poststack.PoststackInversion inversion when dealing with L1 regularized inversion as it makes the inverse problem more stable (no ringing in solution)

  • Changed cost in pylops.optimization.solver.cg and pylops.optimization.solver.cgls to be L2 norms of residuals

  • Fixed pylops.utils.dottest.dottest for imaginary vectors and to ensure u and v vectors are of same dtype of the operator

Version 1.13.0

Released on: 26/03/2021

Version 1.12.0

Released on: 22/11/2020

Version 1.11.1

Released on: 24/10/2020

  • Fixed import of pyfttw when not available in :py:class:``pylops.signalprocessing.ChirpRadon3D`

Version 1.11.0

Released on: 24/10/2020

  • Added pylops.signalprocessing.ChirpRadon2D and pylops.signalprocessing.ChirpRadon3D operators.

  • Fixed bug in the inferred dimensions for regularization data creation in pylops.optimization.leastsquares.NormalEquationsInversion, pylops.optimization.leastsquares.RegularizedInversion, and pylops.optimization.sparsity.SplitBregman.

  • Changed dtype of pylops.HStack to allow automatic inference from dtypes of input operator.

  • Modified dtype of pylops.waveeqprocessing.Marchenko operator to ensure that outputs of forward and adjoint are real arrays.

  • Reverted to previous complex-friendly implementation of pylops.optimization.sparsity._softthreshold to avoid division by 0.

Version 1.10.0

Released on: 13/08/2020

Version 1.9.1

Released on: 25/05/2020

  • Changed internal behaviour of pylops.sparsity.OMP when niter_inner=0. Automatically reverts to Matching Pursuit algorithm.

  • Changed handling of dtype in pylops.signalprocessing.FFT and pylops.signalprocessing.FFT2D to ensure that the type of the input vector is retained when applying forward and adjoint.

  • Added dtype parameter to the FFT calls in the definition of the pylops.waveeqprocessing.MDD operation. This ensure that the type of the real part of G input is enforced to the output vectors of the forward and adjoint operations.

Version 1.9.0

Released on: 13/04/2020

Version 1.8.0

Released on: 12/01/2020

Version 1.7.0

Released on: 10/11/2019

Version 1.6.0

Released on: 10/08/2019

Version 1.5.0

Released on: 30/06/2019

Version 1.4.0

Released on: 01/05/2019

Version 1.3.0

Released on: 24/02/2019

Version 1.2.0

Released on: 13/01/2019

  • Added pylops.LinearOperator.eigs and pylops.LinearOperator.cond methods to estimate estimate eigenvalues and conditioning number using scipy wrapping of ARPACK

  • Modified default dtype for all operators to be float64 (or complex128) to be consistent with default dtypes used by numpy (and scipy) for real and complex floating point numbers.

  • Added pylops.Flip operator

  • Added pylops.Symmetrize operator

  • Added pylops.Block operator

  • Added pylops.Regression operator performing polynomial regression and modified pylops.LinearRegression to be a simple wrapper of pylops.Regression when order=1

  • Modified pylops.MatrixMult operator to work with both numpy ndarrays and scipy sparse matrices

  • Added pylops.avo.prestack.PrestackInversion routine

  • Added possibility to have a data weight via Weight input parameter to pylops.optimization.leastsquares.NormalEquationsInversion and pylops.optimization.leastsquares.RegularizedInversion solvers

  • Added pylops.optimization.sparsity.IRLS solver

Version 1.1.0

Released on: 13/12/2018

Version 1.0.1

Released on: 09/12/2018

  • Changed module from lops to pylops for consistency with library name (and pip install).

  • Removed quickplots from utilities and matplotlib from requirements of PyLops.

Version 1.0.0

Released on: 04/12/2018

  • First official release.

Roadmap

This roadmap is aimed at providing an high-level overview on the bug fixes, improvements and new functionality that are planned for the PyLops library.

Any of the fixes/additions mentioned in the roadmap are directly linked to a Github Issue that provides more details onto the reason and initial thoughts for the implementation of such a fix/addition.

Striked tasks have been completed and related github issue closed with more details regarding how this task has been carried out.

Library structure

  • Create a child repository and python library called geolops (just a suggestion) where geoscience-related operators and examples are moved across, keeping the core pylops library very generic and multi-purpose - Issue #22.

Code cleaning

  • Change all ``np.flatten()`` into ``np.ravel()`` - Issue #24.

  • Fix all if: return ... else: ... statements to enforce a single return with the same number of outputs - Issue #26.

  • Protected attributes and @property attributes in linear operator classes? - Issue #27.

Code optimization

  • Investigate speed-up given by decorating _matvec and _rmatvec methods with numba @jit and @stencil decorators - Issue #23.

  • Replace np.fft.* routines used in several submodules with double engine, numpy and pyFFTW - Issue #20.

Modules

avo
basicoperators
  • Create Kronecker operator - Issue #28.

  • Deal with edges in FirstDerivative and SecondDerivative operators - Issue #34.

optimization
signalprocessing
  • Compare performance in FTT operator of performing np.swap+np.fft.fft(…, axis=-1) versus np.fft.fft(…, axis=chosen) - Issue #33.

  • Add Wavelet operator performing the wavelet transform - Issue #21.

  • Fredholm1 operator applying Fredholm integrals of first kind - Issue #31.

  • Fredholm2 operators applying Fredholm integrals of second kind - Issue #31.

utils

Nothing so far

waveeqprocessing
  • numpy.matmul as a way to speed up integral computation (i.e., inner for loop) in MDC operator - Issue #32.

  • NMO operator performing NMO modelling - Issue #29.

  • WavefieldDecomposition operator performing acoustic wavefield separation by inversion - Issue #30.

Papers using PyLops

This section lists various conference abstracts and papers using the PyLops framework. If you publish a paper using PyLops, we’d love to hear about it!

2021

  • Vakalis, S., D. Chen, M. Yan, and J. Nanzer, 2021, Image enhancement in active incoherent millimeter-wave imaging. Passive and Active Millimeter-Wave Imaging XXIV. doi: 10.1117/12.2585650  
  • Li, X., T. Becker, M. Ravasi, J. Robertsson, and D.-J. van Manen, 2021, Closed-aperture unbounded acoustics experimentation using multidimensional deconvolution. The Journal of the Acoustical Society of America, 149, 1813–1828. doi: 10.1121/10.0003706  
  • Kuijpers, D., I. Vasconcelos, and P. Putzky, 2021, Reconstructing missing seismic data using Deep Learning. arXiv: 2101.09554  
  • Ravasi, M., and C. Birnie, 2021, A Joint Inversion-Segmentation approach to Assisted Seismic Interpretation. arXiv: 2102.03860  
  • Haindl, C., M. Ravasi, and F. Broggini, 2021, Handling gaps in acquisition geometries — Improving Marchenko-based imaging using sparsity-promoting inversion and joint inversion of time-lapse data. Geophysics, 86 (2), S143-S154. doi: 10.1190/geo2020-0036.1  
  • Ulrich, I. E., A. Zunino, C. Boehm, and A. Fichtner, 2021, Sparsifying regularizations for stochastic sample average minimization in ultrasound computed tomography. Medical Imaging 2021: Ultrasonic Imaging and Tomography. doi: 10.1117/12.2580926  
  • Nightingale, James., R. Hayes, A. Kelly, A. Amvrosiadis, A. Etherington, Q. He, N. Li, X. Cao, J. Frawley, S. Cole, A. Enia, C. Frenk, D. Harvey, R. Li, R. Massey, M. Negrello, and A. Robertson, 2020, PyAutoLens: Open-Source Strong Gravitational Lensing. Journal of Open Source Software, 6, 2825. doi: 10.21105/joss.02825  

2020

  • Feng, R., T. Mejer Hansen, D. Grana, and N. Balling, 2020, An unsupervised deep-learning method for porosity estimation based on poststack seismic data. Geophysics, 85 (6), M97-M105. doi: 10.1190/geo2020-0121.1  
  • Zhang, M., 2020, Marchenko Green’s functions from compressive sensing acquisition. SEG Technical Program Expanded Abstracts 2020. doi: 10.1190/segam2020-3424845.1  
  • Vargas, D., and I. Vasconcelos, 2020, Rayleigh-Marchenko Redatuming Using Scattered Fields in Highly Complex Media. EAGE 2020 Annual Conference & Exhibition Online. doi: 10.3997/2214-4609.202011347  
  • Ravasi, M., and I. Vasconcelos, 2020, Implementation of Large-Scale Integral Operators with Modern HPC Solutions. EAGE 2020 Annual Conference & Exhibition Online. doi: 10.3997/2214-4609.202010529  
  • Ravasi, M., and I. Vasconcelos, 2020, PyLops—A linear-operator Python library for scalable algebra and optimization. SoftwareX, 11, 100361. doi: 10.1016/j.softx.2019.100361  

2019

  • Ruan, J., and I. Vasconcelos, 2019, Data- and prior-driven sampling and wavefield reconstruction for sparse, irregularly-sampled, higher-order gradient data. SEG Technical Program Expanded Abstracts 2019. doi: 10.1190/segam2019-3216425.1  

How to cite

When using PyLops in scientific publications, please cite the following paper:

  • Ravasi, M., and I. Vasconcelos, 2020, PyLops—A linear-operator Python library for scalable algebra and optimization. SoftwareX, 11, 100361. doi: 10.1016/j.softx.2019.100361  

Contributors