Merge for v2.2.0

.github/workflows/build.yaml

@@ -7,22 +7,27 @@ jobs:
     strategy:
       matrix:
         platform: [ ubuntu-latest, macos-latest ]
-        python-version: ["3.7", "3.8", "3.9", "3.10"]
+        python-version: ["3.8", "3.9", "3.10"]
 
     runs-on: ${{ matrix.platform }}
     steps:
     - uses: actions/checkout@v3
+    - name: Get history and tags for SCM versioning to work
+      run: |
+        git fetch --prune --unshallow
+        git fetch --depth=1 origin +refs/tags/*:refs/tags/*
     - name: Set up Python ${{ matrix.python-version }}
       uses: actions/setup-python@v4
       with:
         python-version: ${{ matrix.python-version }}
     - name: Install dependencies
       run: |
-        python -m pip install --upgrade pip
+        python -m pip install --upgrade pip setuptools
         pip install flake8 pytest
         if [ -f requirements.txt ]; then pip install -r requirements-dev.txt; fi
     - name: Install pylops
       run: |
+        python -m setuptools_scm
         pip install .
     - name: Test with pytest
       run: |

.github/workflows/codacy-coverage-reporter.yaml

@@ -6,22 +6,27 @@ on: [push, pull_request_target]
 
 jobs:
   build:
+    strategy:
+      matrix:
+        platform: [ ubuntu-latest, ]
+        python-version: ["3.8", ]
 
-    runs-on: ubuntu-latest
+    runs-on: ${{ matrix.platform }}
     steps:
     - uses: actions/checkout@v3
-    - name: Set up Python
+    - name: Set up Python ${{ matrix.python-version }}
       uses: actions/setup-python@v4
       with:
-        python-version: 3.9
+        python-version: ${{ matrix.python-version }}
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
-        pip install flake8 pytest coverage
+        pip install flake8 pytest
         if [ -f requirements.txt ]; then pip install -r requirements-dev.txt; fi
     - name: Install pylops
       run: |
         pip install .
+        pip install coverage
     - name: Code coverage with coverage
       run: |
         coverage run -m pytest

.pre-commit-config.yaml

@@ -15,7 +15,7 @@ repos:
           - --line-length=88
 
   - repo: https://github.com/pycqa/isort
-    rev: 5.10.1
+    rev: 5.12.0
     hooks:
       - id: isort
         name: isort (python)

.readthedocs.yaml

@@ -0,0 +1,23 @@
+# .readthedocs.yaml
+# Read the Docs configuration file
+# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details
+
+# Required
+version: 2
+
+# Set the version of Python and other tools you might need
+build:
+  os: ubuntu-20.04
+  tools:
+    python: "3.9"
+
+# Build documentation in the docs/ directory with Sphinx
+sphinx:
+  configuration: docs/source/conf.py
+
+# Declare the Python requirements required to build your docs
+python:
+  install:
+    - requirements: requirements-dev.txt
+    - method: pip
+      path: .

CHANGELOG.md

@@ -1,3 +1,18 @@
+# 2.2.0
+
+* Added `pylops.signalprocessing.NonStationaryConvolve3D` operator
+* Added nd-array capabilities to `pylops.basicoperators.Identity` and `pylops.basicoperators.Zero`
+* Added second implementation in `pylops.waveeqprocessing.BlendingContinuous` which is more
+  performant when dealing with small number of receivers
+* Added `forceflat` property to operators with ambiguous `rmatvec` (`pylops.basicoperators.Block`,
+  `pylops.basicoperators.Bilinear`, `pylops.basicoperators.BlockDiag`, `pylops.basicoperators.HStack`,
+  `pylops.basicoperators.MatrixMult`, `pylops.basicoperators.VStack`, and `pylops.basicoperators.Zero`)
+* Improved `dynamic` mode of `pylops.waveeqprocessing.Kirchhoff` operator
+* Modified `pylops.signalprocessing.Convolve1D` to allow both filters that are both shorter and longer of the
+  input vector
+* Modified all solvers to use `matvec/rmatvec` instead of `@/.H @` to improve performance
+
+
 # 2.1.0
 * Added `pylops.signalprocessing.DCT`, `pylops.signalprocessing.NonStationaryConvolve1D`,
   `pylops.signalprocessing.NonStationaryConvolve2D`, `pylops.signalprocessing.NonStationaryFilters1D`, and

MIGRATION_V1_V2.md

@@ -10,7 +10,7 @@ should be used as a checklist when converting a piece of code using PyLops from
   for every operator by default. While the change is mostly backwards compatible, there are some operators (e.g. the ``Bilinear``
   transpose/conjugate) which can output reshaped arrays instead of 1d-arrays. To ensure no breakage, you can entirely disable this
   feature either globally by setting ``pylops.set_ndarray_multiplication(False)``, or locally with the context manager
-  ``pylops.disabled_ndarray_multiplication()``. Both will revert to v1.x behavior. At this time, PyLops sparse solvers do
+  ``pylops.disabled_ndarray_multiplication()``. Both will revert to v1.x behavior. At this time, PyLops solvers do
   *not* support N-D array multiplication.
 
   See the table at the end of this document for support ndarray operations.

Makefile

@@ -31,7 +31,7 @@ dev-install_conda:
 
 tests:
 	make pythoncheck
-	$(PYTHON) setup.py test
+	pytest
 
 doc:
 	cd docs  && rm -rf source/api/generated && rm -rf source/gallery &&\

README.md

@@ -147,3 +147,5 @@ A list of video tutorials to learn more about PyLops:
 * Juan Daniel Romero, jdromerom
 * Aniket Singh Rawat, dikwickley
 * Rohan Babbar, rohanbabbar04
+* Wei Zhang, ZhangWeiGeo
+* Fedor Goncharov, fedor-goncharov

azure-pipelines.yml

@@ -55,7 +55,7 @@ jobs:
   steps:
   - task: UsePythonVersion@0
     inputs:
-      versionSpec: '3.7'
+      versionSpec: '3.8'
       architecture: 'x64'
 
   - script: |
@@ -65,7 +65,7 @@ jobs:
     displayName: 'Install prerequisites and library'
 
   - script: |
-      python setup.py test
+      pytest
     condition: succeededOrFailed()
     displayName: 'Run tests'
 
@@ -84,7 +84,7 @@ jobs:
   steps:
   - task: UsePythonVersion@0
     inputs:
-      versionSpec: '3.7'
+      versionSpec: '3.8'
       architecture: 'x64'
 
   - script: |
@@ -94,6 +94,6 @@ jobs:
     displayName: 'Install prerequisites and library'
 
   - script: |
-      python setup.py test
+      pytest
     condition: succeededOrFailed()
     displayName: 'Run tests'

docs/source/api/index.rst

@@ -91,6 +91,7 @@ Signal processing
     ConvolveND
     NonStationaryConvolve1D
     NonStationaryConvolve2D
+    NonStationaryConvolve3D
     NonStationaryFilters1D
     NonStationaryFilters2D
     Interp

docs/source/changelog.rst

@@ -3,6 +3,24 @@
 Changelog
 =========
 
+Version 2.2.0
+-------------
+
+*Released on: 11/11/2023*
+
+* Added :class:`pylops.signalprocessing.NonStationaryConvolve3D` operator
+* Added nd-array capabilities to :class:`pylops.basicoperators.Identity` and :class:`pylops.basicoperators.Zero`
+* Added second implementation in :class:`pylops.waveeqprocessing.BlendingContinuous` which is more
+  performant when dealing with small number of receivers
+* Added `forceflat` property to operators with ambiguous `rmatvec` (:class:`pylops.basicoperators.Block`,
+  :class:`pylops.basicoperators.Bilinear`, :class:`pylops.basicoperators.BlockDiag`, :class:`pylops.basicoperators.HStack`,
+  :class:`pylops.basicoperators.MatrixMult`, :class:`pylops.basicoperators.VStack`, and :class:`pylops.basicoperators.Zero`)
+* Improved `dynamic` mode of :class:`pylops.waveeqprocessing.Kirchhoff` operator
+* Modified :class:`pylops.signalprocessing.Convolve1D` to allow both filters that are both shorter and longer of the
+  input vector
+* Modified all solvers to use `matvec/rmatvec` instead of `@/.H @` to improve performance
+
+
 Version 2.1.0
 -------------
 

docs/source/credits.rst

@@ -18,4 +18,6 @@ Contributors
 *  `Muhammad Izzatullah <https://github.com/izzatum>`_, izzatum
 *  `Juan Daniel Romero <https://github.com/jdromerom>`_, jdromerom
 *  `Aniket Singh Rawat <https://github.com/dikwickley>`_, dikwickley
-*  `Rohan Babbar <https://github.com/rohanbabbar04>`_, rohanbabbar04
+*  `Rohan Babbar <https://github.com/rohanbabbar04>`_, rohanbabbar04
+*  `Wei Zhang <https://github.com/ZhangWeiGeo>`_, ZhangWeiGeo
+*  `Fedor Goncharov <https://github.com/fedor-goncharov>`_, fedor-goncharov

docs/source/extensions.rst

@@ -15,7 +15,8 @@ for academic purposes.
 
 Spin-off projects that aim at extending the capabilities of PyLops are:
 
-* `PyLops-GPU <https://github.com/PyLops/pylops-gpu>`_ : PyLops for GPU arrays (incorporated into PyLops).
-* `PyLops-Distributed <https://github.com/PyLops/pylops-distributed>`_: PyLops for distributed systems with many computing nodes.
+* `PyLops-MPI <https://github.com/PyLops/pylops-mpi>`_: PyLops for distributed systems with many computing nodes using MPI
 * `PyProximal <https://github.com/PyLops/pyproximal>`_: Proximal solvers which integrate with PyLops Linear Operators.
-* `Curvelops <https://github.com/PyLops/curvelops>`_: Python wrapper for the Curvelab 2D and 3D digital curvelet transforms.
+* `Curvelops <https://github.com/PyLops/curvelops>`_: Python wrapper for the Curvelab 2D and 3D digital curvelet transforms.
+* `PyLops-GPU <https://github.com/PyLops/pylops-gpu>`_ : PyLops for GPU arrays (unmantained! the core features are now incorporated into PyLops).
+* `PyLops-Distributed <https://github.com/PyLops/pylops-distributed>`_ : PyLops for distributed systems with many computing nodes using Dask (unmantained!).

docs/source/installation.rst

@@ -487,7 +487,7 @@ or via ``pip`` with
 
 Torch
 -----
-`Torch <http://pytorch.org>`_ used to allow seamless integration between PyLops and PyTorch operators.
+`Torch <http://pytorch.org>`_ is used to allow seamless integration between PyLops and PyTorch operators.
 
 Install it via ``conda`` with:
 

examples/plot_convolve.py

@@ -22,6 +22,7 @@
 """
 import matplotlib.pyplot as plt
 import numpy as np
+import scipy as sp
 from scipy.sparse.linalg import lsqr
 
 import pylops
@@ -132,7 +133,7 @@
 plt.subplots_adjust(top=0.8)
 
 ###############################################################################
-# Finally we do the same using three dimensional signals and
+# We now do the same using three dimensional signals and
 # filters taking advantage of the
 # :py:class:`pylops.signalprocessing.ConvolveND` operator.
 ny, nx, nz = 13, 10, 7
@@ -149,7 +150,7 @@
 xlsqr = xlsqr.reshape(Cop.dims)
 
 fig, axs = plt.subplots(3, 3, figsize=(10, 12))
-fig.suptitle("Convolve 3d data", y=0.95, fontsize=14, fontweight="bold")
+fig.suptitle("Convolve 3d data", y=0.98, fontsize=14, fontweight="bold")
 axs[0][0].imshow(x[ny // 3], cmap="gray", vmin=-1, vmax=1)
 axs[0][1].imshow(y[ny // 3], cmap="gray", vmin=-1, vmax=1)
 axs[0][2].imshow(xlsqr[ny // 3], cmap="gray", vmin=-1, vmax=1)
@@ -172,3 +173,115 @@
 axs[2][1].axis("tight")
 axs[2][2].axis("tight")
 plt.tight_layout()
+
+###############################################################################
+# Up until now, we have only considered the case where the filter is compact
+# and much shorter of the input data. There are however scenarios where the
+# opposite (i.e., filter is longer than the signal) is desired. This is for
+# example the case when one wants to estimate a filter (:math:`\mathbf{h}`)
+# to match two signals (:math:`\mathbf{x}` and :math:`\mathbf{y}`):
+#
+# .. math::
+#   J = || \mathbf{y} - \mathbf{X} \mathbf{h} ||_2^2
+#
+# Such a scenario is very commonly used in so-called adaptive substraction
+# techniques. We will try now to use :py:class:`pylops.signalprocessing.Convolve1D`
+# to match two signals that have both a phase and amplitude mismatch.
+
+# Define input signal
+nt = 101
+dt = 0.004
+t = np.arange(nt) * dt
+
+x = np.zeros(nt)
+x[[int(nt / 4), int(nt / 2), int(2 * nt / 3)]] = [3, -2, 1]
+h, th, hcenter = ricker(t[:41], f0=20)
+
+Cop = pylops.signalprocessing.Convolve1D(nt, h=h, offset=hcenter, dtype="float32")
+x = Cop * x
+
+# Phase and amplitude corrupt the input
+amp = 0.9
+phase = 40
+
+y = amp * (
+    x * np.cos(np.deg2rad(phase))
+    + np.imag(sp.signal.hilbert(x)) * np.sin(np.deg2rad(phase))
+)
+
+# Define convolution operator
+nh = 21
+th = np.arange(nh) * dt - dt * nh // 2
+Yop = pylops.signalprocessing.Convolve1D(nh, h=y, offset=nh // 2)
+
+# Find filter to match x to y
+h = Yop / x
+ymatched = Yop @ h
+
+# Find sparse filter to match x to y
+hsparse = pylops.optimization.sparsity.fista(Yop, x, niter=100, eps=1e-1)[0]
+ymatchedsparse = Yop @ hsparse
+
+fig, ax = plt.subplots(1, 1, figsize=(10, 3))
+ax.plot(t, x, "k", lw=2, label=r"$x$")
+ax.plot(t, y, "r", lw=2, label=r"$y$")
+ax.plot(t, ymatched, "--g", lw=2, label=r"$y_{matched}$")
+ax.plot(t, x - ymatched, "--k", lw=2, label=r"$x-y_{matched,sparse}$")
+ax.plot(t, ymatchedsparse, "--m", lw=2, label=r"$y_{matched}$")
+ax.plot(t, x - ymatchedsparse, "--k", lw=2, label=r"$x-y_{matched,sparse}$")
+ax.set_title("Signals to match", fontsize=14, fontweight="bold")
+ax.legend(loc="upper right")
+plt.tight_layout()
+
+fig, axs = plt.subplots(1, 2, figsize=(10, 3))
+axs[0].plot(th, h, "k", lw=2)
+axs[0].set_title("Matching filter", fontsize=14, fontweight="bold")
+axs[1].plot(th, hsparse, "k", lw=2)
+axs[1].set_title("Sparse Matching filter", fontsize=14, fontweight="bold")
+plt.tight_layout()
+
+
+###############################################################################
+# Finally, in some cases one wants to convolve (or correlate) two signals of
+# the same size. This can also be obtained using
+# :py:class:`pylops.signalprocessing.Convolve1D`. We will see here a case
+# where the operator is used to trace-wise auto-correlate signals from
+# a 2-dimensional array representing a seismic dataset.
+
+# Create data
+par = {"ox": -140, "dx": 2, "nx": 140, "ot": 0, "dt": 0.004, "nt": 200, "f0": 20}
+
+v = 1500
+t0 = [0.2, 0.4, 0.5]
+px = [0, 0, 0]
+pxx = [1e-5, 5e-6, 1e-20]
+amp = [1.0, -2, 0.5]
+
+t, t2, x, y = pylops.utils.seismicevents.makeaxis(par)
+wav = pylops.utils.wavelets.ricker(t[:41], f0=par["f0"])[0]
+_, data = pylops.utils.seismicevents.parabolic2d(x, t, t0, px, pxx, amp, wav)
+
+# Convolution operator
+Xop = pylops.signalprocessing.Convolve1D(
+    (par["nx"], par["nt"]), h=data, offset=par["nt"] // 2, axis=-1, method="fft"
+)
+
+corr = Xop.H @ data
+
+fig, axs = plt.subplots(1, 2, figsize=(10, 6))
+axs[0].imshow(data.T, cmap="gray", vmin=-1, vmax=1, extent=(x[0], x[-1], t[-1], t[0]))
+axs[0].set_xlabel("Rec (m)", fontsize=14, fontweight="bold")
+axs[0].set_ylabel("T (s)", fontsize=14, fontweight="bold")
+axs[0].set_title("Data", fontsize=14, fontweight="bold")
+axs[0].axis("tight")
+axs[1].imshow(
+    corr.T,
+    cmap="gray",
+    vmin=-10,
+    vmax=10,
+    extent=(x[0], x[-1], t[par["nt"] // 2], -t[par["nt"] // 2]),
+)
+axs[1].set_xlabel("Rec (m)", fontsize=14, fontweight="bold")
+axs[1].set_title("Auto-correlation", fontsize=14, fontweight="bold")
+axs[1].axis("tight")
+plt.tight_layout()

examples/plot_ista.py

@@ -106,7 +106,7 @@
 x[int(nt / 2)] = 1
 x[nt - 20] = 0.5
 
-h, th, hcenter = pylops.utils.wavelets.ricker(t[:101], f0=20)
+h, th, hcenter = pylops.utils.wavelets.ricker(t[:21], f0=20)
 
 Cop = pylops.signalprocessing.Convolve1D(nt, h=h, offset=hcenter, dtype="float32")
 y = Cop * x