diff --git a/.github/workflows/build.yaml b/.github/workflows/build.yaml index f6852e4c..abf7ff7e 100644 --- a/.github/workflows/build.yaml +++ b/.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: | diff --git a/.github/workflows/codacy-coverage-reporter.yaml b/.github/workflows/codacy-coverage-reporter.yaml index b0a8f90d..b16411b4 100644 --- a/.github/workflows/codacy-coverage-reporter.yaml +++ b/.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 diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 14011ef6..a4ca2e16 100644 --- a/.pre-commit-config.yaml +++ b/.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) diff --git a/.readthedocs.yaml b/.readthedocs.yaml new file mode 100644 index 00000000..a74e3480 --- /dev/null +++ b/.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: . diff --git a/CHANGELOG.md b/CHANGELOG.md index d632bd4a..93ce3982 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/MIGRATION_V1_V2.md b/MIGRATION_V1_V2.md index e198965d..d26271b8 100644 --- a/MIGRATION_V1_V2.md +++ b/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. diff --git a/Makefile b/Makefile index 0ac4b3d5..2db298d1 100755 --- a/Makefile +++ b/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 &&\ diff --git a/README.md b/README.md index a9258eef..3d1fc89f 100644 --- a/README.md +++ b/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 diff --git a/azure-pipelines.yml b/azure-pipelines.yml index d778e70a..ef5306bb 100644 --- a/azure-pipelines.yml +++ b/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' diff --git a/docs/source/api/index.rst b/docs/source/api/index.rst index ea25d5bf..9d7c3f66 100755 --- a/docs/source/api/index.rst +++ b/docs/source/api/index.rst @@ -91,6 +91,7 @@ Signal processing ConvolveND NonStationaryConvolve1D NonStationaryConvolve2D + NonStationaryConvolve3D NonStationaryFilters1D NonStationaryFilters2D Interp diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index a56dd9aa..ad3d94b1 100644 --- a/docs/source/changelog.rst +++ b/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 ------------- diff --git a/docs/source/credits.rst b/docs/source/credits.rst index 2a8558d1..2a5daf61 100755 --- a/docs/source/credits.rst +++ b/docs/source/credits.rst @@ -18,4 +18,6 @@ Contributors * `Muhammad Izzatullah `_, izzatum * `Juan Daniel Romero `_, jdromerom * `Aniket Singh Rawat `_, dikwickley -* `Rohan Babbar `_, rohanbabbar04 \ No newline at end of file +* `Rohan Babbar `_, rohanbabbar04 +* `Wei Zhang `_, ZhangWeiGeo +* `Fedor Goncharov `_, fedor-goncharov \ No newline at end of file diff --git a/docs/source/extensions.rst b/docs/source/extensions.rst index 57c5c8d0..25220df5 100755 --- a/docs/source/extensions.rst +++ b/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 `_ : PyLops for GPU arrays (incorporated into PyLops). -* `PyLops-Distributed `_: PyLops for distributed systems with many computing nodes. +* `PyLops-MPI `_: PyLops for distributed systems with many computing nodes using MPI * `PyProximal `_: Proximal solvers which integrate with PyLops Linear Operators. -* `Curvelops `_: Python wrapper for the Curvelab 2D and 3D digital curvelet transforms. \ No newline at end of file +* `Curvelops `_: Python wrapper for the Curvelab 2D and 3D digital curvelet transforms. +* `PyLops-GPU `_ : PyLops for GPU arrays (unmantained! the core features are now incorporated into PyLops). +* `PyLops-Distributed `_ : PyLops for distributed systems with many computing nodes using Dask (unmantained!). diff --git a/docs/source/installation.rst b/docs/source/installation.rst index 2e22eb2d..3b237f23 100755 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -487,7 +487,7 @@ or via ``pip`` with Torch ----- -`Torch `_ used to allow seamless integration between PyLops and PyTorch operators. +`Torch `_ is used to allow seamless integration between PyLops and PyTorch operators. Install it via ``conda`` with: diff --git a/examples/plot_convolve.py b/examples/plot_convolve.py index 19431ff3..4afa9604 100644 --- a/examples/plot_convolve.py +++ b/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() diff --git a/examples/plot_ista.py b/examples/plot_ista.py index e37f2247..d6359905 100755 --- a/examples/plot_ista.py +++ b/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 diff --git a/examples/plot_sliding.py b/examples/plot_sliding.py index 2d4be7f8..ad553487 100644 --- a/examples/plot_sliding.py +++ b/examples/plot_sliding.py @@ -29,12 +29,14 @@ ############################################################################### # Let's start by creating a 1-dimensional array of size :math:`n_t` and create # a sliding operator to compute its transformed representation. -nwins = 4 + +# sliding window parameters nwin = 26 nover = 3 nop = 64 -dimd = nwin * nwins - 3 * nover +# length of input signal (chosen to ensure perfect match with sliding windows) +dimd = 95 t = np.arange(dimd) * 0.004 data = np.sin(2 * np.pi * 20 * t) diff --git a/pylops/avo/poststack.py b/pylops/avo/poststack.py index 11c49682..7e514707 100644 --- a/pylops/avo/poststack.py +++ b/pylops/avo/poststack.py @@ -102,7 +102,7 @@ def _PoststackLinearModelling( # Create wavelet operator if len(wav.shape) == 1: - C = convmtx(wav, nt0)[:, len(wav) // 2 : -len(wav) // 2 + 1] + C = convmtx(wav, nt0, len(wav) // 2)[:nt0] else: C = nonstationary_convmtx(wav, nt0, hc=wav.shape[1] // 2, pad=(nt0, nt0)) # Combine operators diff --git a/pylops/avo/prestack.py b/pylops/avo/prestack.py index 491dfb9c..8630bc9a 100644 --- a/pylops/avo/prestack.py +++ b/pylops/avo/prestack.py @@ -191,7 +191,7 @@ def PrestackLinearModelling( D = get_block_diag(theta)(*([D] * nG)) # Create wavelet operator - C = ncp.asarray(convmtx(wav, nt0))[:, len(wav) // 2 : -len(wav) // 2 + 1] + C = ncp.asarray(convmtx(wav, nt0, len(wav) // 2)[:nt0]) C = [C] * ntheta C = get_block_diag(theta)(*C) @@ -346,7 +346,7 @@ def PrestackWaveletModelling( M = ncp.dot(G, ncp.dot(D, m.T.ravel())).reshape(ntheta, nt0) Mconv = VStack( [ - MatrixMult(convmtx(M[itheta], nwav)[wavc : -nwav + wavc + 1], dtype=dtype) + MatrixMult(convmtx(M[itheta], nwav, wavc)[:nt0], dtype=dtype) for itheta in range(ntheta) ] ) diff --git a/pylops/basicoperators/block.py b/pylops/basicoperators/block.py index 88e989bd..15eea4a8 100644 --- a/pylops/basicoperators/block.py +++ b/pylops/basicoperators/block.py @@ -1,7 +1,6 @@ __all__ = ["Block"] import multiprocessing as mp - from typing import Iterable, Optional from pylops import LinearOperator @@ -15,12 +14,18 @@ class _Block(LinearOperator): Used to be able to provide operators from different libraries to Block. """ - def __init__(self, ops: Iterable[Iterable[LinearOperator]], - dtype: Optional[DTypeLike] = None, - _HStack=HStack, - _VStack=VStack, - args_HStack: Optional[dict] = None, - args_VStack: Optional[dict] = None, name: str = 'B'): + + def __init__( + self, + ops: Iterable[Iterable[LinearOperator]], + forceflat: bool = None, + dtype: Optional[DTypeLike] = None, + _HStack=HStack, + _VStack=VStack, + args_HStack: Optional[dict] = None, + args_VStack: Optional[dict] = None, + name: str = "B", + ): if args_HStack is None: self.args_HStack = {} else: @@ -30,7 +35,12 @@ def __init__(self, ops: Iterable[Iterable[LinearOperator]], else: self.args_VStack = args_VStack hblocks = [_HStack(hblock, dtype=dtype, **self.args_HStack) for hblock in ops] - super().__init__(Op=_VStack(ops=hblocks, dtype=dtype, **self.args_VStack), name=name) + super().__init__( + Op=_VStack( + ops=hblocks, forceflat=forceflat, dtype=dtype, **self.args_VStack + ), + name=name, + ) def _matvec(self, x: NDArray) -> NDArray: return super()._matvec(x) @@ -53,6 +63,10 @@ class Block(_Block): nproc : :obj:`int`, optional Number of processes used to evaluate the N operators in parallel using ``multiprocessing``. If ``nproc=1``, work in serial mode. + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after rmatvec. dtype : :obj:`str`, optional Type of elements in input array. @@ -123,9 +137,16 @@ class Block(_Block): \end{bmatrix} """ - def __init__(self, ops: Iterable[Iterable[LinearOperator]], - nproc: int = 1, - dtype: Optional[DTypeLike] = None): + + def __init__( + self, + ops: Iterable[Iterable[LinearOperator]], + nproc: int = 1, + forceflat: bool = None, + dtype: Optional[DTypeLike] = None, + ): if nproc > 1: self.pool = mp.Pool(processes=nproc) - super().__init__(ops=ops, dtype=dtype, args_VStack={"nproc": nproc}) + super().__init__( + ops=ops, forceflat=forceflat, dtype=dtype, args_VStack={"nproc": nproc} + ) diff --git a/pylops/basicoperators/blockdiag.py b/pylops/basicoperators/blockdiag.py index 2447ffcb..e13ed026 100644 --- a/pylops/basicoperators/blockdiag.py +++ b/pylops/basicoperators/blockdiag.py @@ -44,6 +44,10 @@ class BlockDiag(LinearOperator): nproc : :obj:`int`, optional Number of processes used to evaluate the N operators in parallel using ``multiprocessing``. If ``nproc=1``, work in serial mode. + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after matvec and rmatvec. dtype : :obj:`str`, optional Type of elements in input array. @@ -108,6 +112,7 @@ def __init__( self, ops: Sequence[LinearOperator], nproc: int = 1, + forceflat: bool = None, dtype: Optional[DTypeLike] = None, ) -> None: self.ops = ops @@ -122,6 +127,22 @@ def __init__( self.mops = int(mops.sum()) self.nnops = np.insert(np.cumsum(nops), 0, 0) self.mmops = np.insert(np.cumsum(mops), 0, 0) + # define dims (check if all operators have the same, + # otherwise make same as self.mops and forceflat=True) + dims = [op.dims for op in self.ops] + if len(set(dims)) == 1: + dims = (len(ops), *dims[0]) + else: + dims = (self.mops,) + forceflat = True + # define dimsd (check if all operators have the same, + # otherwise make same as self.nops and forceflat=True) + dimsd = [op.dimsd for op in self.ops] + if len(set(dimsd)) == 1: + dimsd = (len(ops), *dimsd[0]) + else: + dimsd = (self.nops,) + forceflat = True # create pool for multiprocessing self._nproc = nproc self.pool: Optional[mp.pool.Pool] = None @@ -130,7 +151,13 @@ def __init__( dtype = _get_dtype(ops) if dtype is None else np.dtype(dtype) clinear = all([getattr(oper, "clinear", True) for oper in self.ops]) - super().__init__(dtype=dtype, shape=(self.nops, self.mops), clinear=clinear) + super().__init__( + dtype=dtype, + dims=dims, + dimsd=dimsd, + clinear=clinear, + forceflat=forceflat, + ) @property def nproc(self) -> int: diff --git a/pylops/basicoperators/hstack.py b/pylops/basicoperators/hstack.py index c97fac1c..5cfbbec0 100644 --- a/pylops/basicoperators/hstack.py +++ b/pylops/basicoperators/hstack.py @@ -44,6 +44,10 @@ class HStack(LinearOperator): nproc : :obj:`int`, optional Number of processes used to evaluate the N operators in parallel using ``multiprocessing``. If ``nproc=1``, work in serial mode. + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after matvec. dtype : :obj:`str`, optional Type of elements in input array. @@ -107,6 +111,7 @@ def __init__( self, ops: Sequence[LinearOperator], nproc: int = 1, + forceflat: bool = None, dtype: Optional[str] = None, ) -> None: self.ops = ops @@ -121,15 +126,28 @@ def __init__( raise ValueError("operators have different number of rows") self.nops = int(nops[0]) self.mmops = np.insert(np.cumsum(mops), 0, 0) + # define dimsd (check if all operators have the same, + # otherwise make same as self.nops and forceflat=True) + dimsd = [op.dimsd for op in self.ops] + if len(set(dimsd)) == 1: + dimsd = dimsd[0] + else: + dimsd = (self.nops,) + forceflat = True # create pool for multiprocessing self._nproc = nproc self.pool = None if self.nproc > 1: self.pool = mp.Pool(processes=nproc) - dtype = _get_dtype(self.ops) if dtype is None else np.dtype(dtype) clinear = all([getattr(oper, "clinear", True) for oper in self.ops]) - super().__init__(dtype=dtype, shape=(self.nops, self.mops), clinear=clinear) + super().__init__( + dtype=dtype, + shape=(self.nops, self.mops), + dimsd=dimsd, + clinear=clinear, + forceflat=forceflat, + ) @property def nproc(self) -> int: diff --git a/pylops/basicoperators/identity.py b/pylops/basicoperators/identity.py index a41f5ce9..c2d05a30 100644 --- a/pylops/basicoperators/identity.py +++ b/pylops/basicoperators/identity.py @@ -1,13 +1,14 @@ __all__ = ["Identity"] -from typing import Optional +from typing import Optional, Union import numpy as np from pylops import LinearOperator from pylops.utils.backend import get_array_module -from pylops.utils.typing import DTypeLike, NDArray +from pylops.utils.decorators import reshaped +from pylops.utils.typing import DTypeLike, InputDimsLike, NDArray class Identity(LinearOperator): @@ -19,16 +20,29 @@ class Identity(LinearOperator): removes last :math:`N - M` elements from data in adjoint and pads with :math:`0` in forward. + Note that the identity operator can handle both 1d and nd arrays; in the + case of nd arrays, all elements of N must be larger or equal than those of M + (or all elements of M must be larger or equal than those of N). + Parameters ---------- - N : :obj:`int` + N : :obj:`int` or :obj:`tuple` Number of samples in data (and model, if ``M`` is not provided). - M : :obj:`int`, optional - Number of samples in model. + If a tuple is provided, this is interpreted as the data (and model) + are nd-arrays. + M : :obj:`int` or :obj:`tuple`, optional + Number of samples in model. If a tuple is provided, this is interpreted + as the model is an nd-array. Note that when `M` is a tuple, `N` must be + also a tuple with the same number of elements. inplace : :obj:`bool`, optional Work inplace (``True``) or make a new copy (``False``). By default, data is a reference to the model (in forward) and model is a reference to the data (in adjoint). + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after matvec and rmatvec. Note that this is only + required when `N` and `M` are tuples (input and output arrays are nd-arrays). dtype : :obj:`str`, optional Type of elements in input array. name : :obj:`str`, optional @@ -44,6 +58,15 @@ class Identity(LinearOperator): Operator contains a matrix that can be solved explicitly (``True``) or not (``False``) + Raises + ------ + ValueError + - If ``M`` is a tuple with different number of elements of ``N`` + - If ``N`` ``M`` are non-identical tuples and some values are largers + in ``N`` and some in ``M`` + NotImplementedError + If ``N`` or ``M`` have type different from int or tuple/list + Notes ----- For :math:`M = N`, an *Identity* operator simply moves the model @@ -87,38 +110,90 @@ class Identity(LinearOperator): def __init__( self, - N: int, - M: Optional[int] = None, + N: Union[int, InputDimsLike], + M: Optional[Union[int, InputDimsLike]] = None, inplace: bool = True, + forceflat: bool = None, dtype: DTypeLike = "float64", name: str = "I", ) -> None: M = N if M is None else M - super().__init__(dtype=np.dtype(dtype), shape=(N, M), name=name) + if isinstance(N, int) and isinstance(M, int): + # N and M are scalars (1d-arrays) + super().__init__( + dtype=np.dtype(dtype), + dims=(M,), + dimsd=(N,), + forceflat=forceflat, + name=name, + ) + # identify behaviour for matvec/rmatvec: 'same' for N=M, + # 'data' for N>M, and 'model' for M>N + if N == M: + self.mode = "same" + elif N < M: + self.mode = "model" + self.sliceN = slice(0, N) + self.sliceM = slice(0, M) + else: + self.mode = "data" + self.sliceN = slice(0, N) + self.sliceM = slice(0, M) + elif isinstance(N, (tuple, list)) and isinstance(M, (tuple, list)): + # N and M are tuples (nd-arrays) + # First check that all elements in N and M are the same or that + # all elements of either N or M are bigger than the other one and + # raise error is not the case + if np.array_equal(N, M): + self.mode = "same" + elif np.array_equal(M, np.maximum(N, M)): + self.mode = "model" + self.sliceN = tuple([slice(0, n) for n in N]) + self.sliceM = tuple([slice(0, m) for m in M]) + elif np.array_equal(N, np.maximum(N, M)): + self.mode = "data" + self.sliceN = tuple([slice(0, n) for n in N]) + self.sliceM = tuple([slice(0, m) for m in M]) + else: + raise ValueError( + "N and M are not identical, " + "and some values are larger in N and some in M" + ) + super().__init__( + dtype=np.dtype(dtype), dims=M, dimsd=N, forceflat=forceflat, name=name + ) + else: + raise NotImplementedError( + f"N and M must have same type and equal to " + f"int, tuple, or list, instead their types" + f" are type(N)={type(N)} and type(M)={type(M)}" + ) self.inplace = inplace + @reshaped def _matvec(self, x: NDArray) -> NDArray: ncp = get_array_module(x) if not self.inplace: x = x.copy() - if self.shape[0] == self.shape[1]: + if self.mode == "same": y = x - elif self.shape[0] < self.shape[1]: - y = x[: self.shape[0]] + elif self.mode == "model": + y = x[self.sliceN] else: - y = ncp.zeros(self.shape[0], dtype=self.dtype) - y[: self.shape[1]] = x + y = ncp.zeros(self.dimsd, dtype=self.dtype) + y[self.sliceM] = x return y + @reshaped def _rmatvec(self, x: NDArray) -> NDArray: ncp = get_array_module(x) if not self.inplace: x = x.copy() - if self.shape[0] == self.shape[1]: + if self.mode == "same": y = x - elif self.shape[0] < self.shape[1]: - y = ncp.zeros(self.shape[1], dtype=self.dtype) - y[: self.shape[0]] = x + elif self.mode == "model": + y = ncp.zeros(self.dims, dtype=self.dtype) + y[self.sliceN] = x else: - y = x[: self.shape[1]] + y = x[self.sliceM] return y diff --git a/pylops/basicoperators/matrixmult.py b/pylops/basicoperators/matrixmult.py index 65a62d4f..a5f713ab 100644 --- a/pylops/basicoperators/matrixmult.py +++ b/pylops/basicoperators/matrixmult.py @@ -29,6 +29,12 @@ class MatrixMult(LinearOperator): Number of samples for each other dimension of model (model/data will be reshaped and ``A`` applied multiple times to each column of the model/data). + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after matvec and rmatvec. Note that this is only + required when `otherdims=None`, otherwise pylops will detect whether to + return a 1d or nd array. dtype : :obj:`str`, optional Type of elements in input array. name : :obj:`str`, optional @@ -56,6 +62,7 @@ def __init__( self, A: NDArray, otherdims: Optional[Union[int, InputDimsLike]] = None, + forceflat: bool = None, dtype: DTypeLike = "float64", name: str = "M", ) -> None: @@ -81,12 +88,26 @@ def __init__( self.reshape = True explicit = False + # Check if forceflat is needed and set it back to None otherwise + if otherdims is not None and forceflat is not None: + logging.warning( + "setting forceflat=None since otherdims!=None. " + "PyLops will automatically detect whether to return " + "a 1d or nd array based on the shape of the input " + "array." + ) + forceflat = None # Check dtype for correctness (upcast to complex when A is complex) if np.iscomplexobj(A) and not np.iscomplexobj(np.ones(1, dtype=dtype)): dtype = A.dtype logging.warning("Matrix A is a complex object, dtype cast to %s" % dtype) super().__init__( - dtype=np.dtype(dtype), dims=dims, dimsd=dimsd, explicit=explicit, name=name + dtype=np.dtype(dtype), + dims=dims, + dimsd=dimsd, + explicit=explicit, + forceflat=forceflat, + name=name, ) def _matvec(self, x: NDArray) -> NDArray: diff --git a/pylops/basicoperators/restriction.py b/pylops/basicoperators/restriction.py index 4aaf10ea..c2e51a31 100644 --- a/pylops/basicoperators/restriction.py +++ b/pylops/basicoperators/restriction.py @@ -1,5 +1,7 @@ __all__ = ["Restriction"] +import logging + from typing import Sequence, Union import numpy as np @@ -11,6 +13,8 @@ from pylops.utils.backend import get_array_module, to_cupy_conditional from pylops.utils.typing import DTypeLike, InputDimsLike, IntNDArray, NDArray +logging.basicConfig(format="%(levelname)s: %(message)s", level=logging.WARNING) + def _compute_iavamask(dims, axis, iava, ncp): """Compute restriction mask when using cupy arrays""" @@ -46,6 +50,12 @@ class Restriction(LinearOperator): Work inplace (``True``) or make a new copy (``False``). By default, data is a reference to the model (in forward) and model is a reference to the data (in adjoint). + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after rmatvec. Note that this is only + required when `len(dims)=2`, otherwise pylops will detect whether to + return a 1d or nd array. dtype : :obj:`str`, optional Type of elements in input array. name : :obj:`str`, optional @@ -97,6 +107,7 @@ def __init__( iava: Union[IntNDArray, Sequence[int]], axis: int = -1, inplace: bool = True, + forceflat: bool = None, dtype: DTypeLike = "float64", name: str = "R", ) -> None: @@ -105,7 +116,20 @@ def __init__( axis = normalize_axis_index(axis, len(dims)) dimsd = list(dims) # data dimensions dimsd[axis] = len(iava) - super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dimsd, name=name) + + # check if forceflat is needed and set it back to None otherwise + if len(dims) > 2: + if forceflat is not None: + logging.warning( + f"setting forceflat=None since len(dims)={len(dims)}>2. " + f"PyLops will automatically detect whether to return " + f"a 1d or nd array based on the shape of the input" + f"array." + ) + forceflat = None + + super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dimsd, + forceflat=forceflat, name=name) iavareshape = np.ones(len(self.dims), dtype=int) iavareshape[axis] = len(iava) diff --git a/pylops/basicoperators/sum.py b/pylops/basicoperators/sum.py index ec42cabe..855cf521 100644 --- a/pylops/basicoperators/sum.py +++ b/pylops/basicoperators/sum.py @@ -1,5 +1,7 @@ __all__ = ["Sum"] +import logging + import numpy as np from pylops import LinearOperator @@ -8,6 +10,8 @@ from pylops.utils.decorators import reshaped from pylops.utils.typing import DTypeLike, InputDimsLike, NDArray +logging.basicConfig(format="%(levelname)s: %(message)s", level=logging.WARNING) + class Sum(LinearOperator): r"""Sum operator. @@ -24,6 +28,12 @@ class Sum(LinearOperator): .. versionadded:: 2.0.0 Axis along which model is summed. + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after rmatvec. Note that this is only + required when `len(dims)=2`, otherwise pylops will detect whether to + return a 1d or nd array. dtype : :obj:`str`, optional Type of elements in input array. name : :obj:`str`, optional @@ -61,6 +71,7 @@ def __init__( self, dims: InputDimsLike, axis: int = -1, + forceflat: bool = None, dtype: DTypeLike = "float64", name: str = "S", ) -> None: @@ -71,7 +82,22 @@ def __init__( # data dimensions dimsd = list(dims).copy() dimsd.pop(self.axis) - super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dimsd, name=name) + # check if forceflat is needed and set it back to None otherwise + if len(dims) > 2 and forceflat is not None: + logging.warning( + f"setting forceflat=None since len(dims)={len(dims)}>2. " + f"PyLops will automatically detect whether to return " + f"a 1d or nd array based on the shape of the input" + f"array." + ) + forceflat = None + super().__init__( + dtype=np.dtype(dtype), + dims=dims, + dimsd=dimsd, + forceflat=forceflat, + name=name, + ) # array of ones with dims of model in self.axis for np.tile in adjoint mode self.tile = np.ones(len(self.dims), dtype=int) diff --git a/pylops/basicoperators/vstack.py b/pylops/basicoperators/vstack.py index dc6734f8..812b1a7e 100644 --- a/pylops/basicoperators/vstack.py +++ b/pylops/basicoperators/vstack.py @@ -44,6 +44,10 @@ class VStack(LinearOperator): nproc : :obj:`int`, optional Number of processes used to evaluate the N operators in parallel using ``multiprocessing``. If ``nproc=1``, work in serial mode. + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after rmatvec. dtype : :obj:`str`, optional Type of elements in input array. @@ -107,6 +111,7 @@ def __init__( self, ops: Sequence[LinearOperator], nproc: int = 1, + forceflat: bool = None, dtype: Optional[DTypeLike] = None, ) -> None: self.ops = ops @@ -121,15 +126,28 @@ def __init__( raise ValueError("operators have different number of columns") self.mops = int(mops[0]) self.nnops = np.insert(np.cumsum(nops), 0, 0) + # define dims (check if all operators have the same, + # otherwise make same as self.mops and forceflat=True) + dims = [op.dims for op in self.ops] + if len(set(dims)) == 1: + dims = dims[0] + else: + dims = (self.mops,) + forceflat = True # create pool for multiprocessing self._nproc = nproc self.pool = None if self.nproc > 1: self.pool = mp.Pool(processes=nproc) - dtype = _get_dtype(self.ops) if dtype is None else np.dtype(dtype) clinear = all([getattr(oper, "clinear", True) for oper in self.ops]) - super().__init__(dtype=dtype, shape=(self.nops, self.mops), clinear=clinear) + super().__init__( + dtype=dtype, + shape=(self.nops, self.mops), + dims=dims, + clinear=clinear, + forceflat=forceflat, + ) @property def nproc(self) -> int: diff --git a/pylops/basicoperators/zero.py b/pylops/basicoperators/zero.py index 04e5f12c..eb375d9a 100644 --- a/pylops/basicoperators/zero.py +++ b/pylops/basicoperators/zero.py @@ -1,12 +1,12 @@ __all__ = ["Zero"] -from typing import Optional +from typing import Optional, Union import numpy as np from pylops import LinearOperator from pylops.utils.backend import get_array_module -from pylops.utils.typing import DTypeLike, NDArray +from pylops.utils.typing import DTypeLike, InputDimsLike, NDArray class Zero(LinearOperator): @@ -17,10 +17,19 @@ class Zero(LinearOperator): Parameters ---------- - N : :obj:`int` - Number of samples in data (and model in M is not provided). - M : :obj:`int`, optional - Number of samples in model. + N : :obj:`int` or :obj:`tuple` + Number of samples in data (and model, if ``M`` is not provided). + If a tuple is provided, this is interpreted as the data (and model) + are nd-arrays. + M : :obj:`int` or :obj:`tuple`, optional + Number of samples in model. If a tuple is provided, this is interpreted + as the model is an nd-array. Note that when `M` is a tuple, `N` must be + also a tuple with the same number of elements. + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after matvec and rmatvec. Note that this is only + required when `N` and `M` are tuples (input and output arrays are nd-arrays). dtype : :obj:`str`, optional Type of elements in input array. name : :obj:`str`, optional @@ -55,13 +64,32 @@ class Zero(LinearOperator): def __init__( self, - N: int, - M: Optional[int] = None, + N: Union[int, InputDimsLike], + M: Optional[Union[int, InputDimsLike]] = None, + forceflat: bool = None, dtype: DTypeLike = "float64", name: str = "Z", ) -> None: M = N if M is None else M - super().__init__(dtype=np.dtype(dtype), shape=(N, M), name=name) + if isinstance(N, int) and isinstance(M, int): + # N and M are scalars (1d-arrays) + dims, dimsd = (M,), (N,) + elif isinstance(N, (tuple, list)) and isinstance(M, (tuple, list)): + # N and M are tuples (nd-arrays) + dims, dimsd = M, N + else: + raise NotImplementedError( + f"N and M must have same type and equal to " + f"int, tuple, or list, instead their types" + f" are type(N)={type(N)} and type(M)={type(M)}" + ) + super().__init__( + dtype=np.dtype(dtype), + dims=dims, + dimsd=dimsd, + forceflat=forceflat, + name=name, + ) def _matvec(self, x: NDArray) -> NDArray: ncp = get_array_module(x) diff --git a/pylops/linearoperator.py b/pylops/linearoperator.py index b16358cb..0a719cf3 100644 --- a/pylops/linearoperator.py +++ b/pylops/linearoperator.py @@ -85,13 +85,20 @@ class LinearOperator(_LinearOperator): .. versionadded:: 2.0.0 Dimensions of data. If not provided, ``(self.shape[0],)`` is used. - explicit : :obj:`bool`, optional - Operator contains a matrix that can be solved explicitly - (``True``) or not (``False``). Defaults to ``False``. clinear : :obj:`bool`, optional .. versionadded:: 1.17.0 Operator is complex-linear. Defaults to ``True``. + explicit : :obj:`bool`, optional + Operator contains a matrix that can be solved explicitly + (``True``) or not (``False``). Defaults to ``False``. + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after matvec/rmatvec if the input is ambiguous + (i.e., is a 1D array both when operating with ND arrays and with 1D arrays). + Defaults to ``None`` for operators that have no ambiguity or to ``True`` + for those with ambiguity. name : :obj:`str`, optional .. versionadded:: 2.0.0 @@ -108,6 +115,7 @@ def __init__( dimsd: Optional[ShapeLike] = None, clinear: Optional[bool] = None, explicit: Optional[bool] = None, + forceflat: Optional[bool] = None, name: Optional[str] = None, ) -> None: if Op is not None: @@ -124,6 +132,9 @@ def __init__( ) if explicit and hasattr(Op, "A"): self.A = Op.A + forceflat = ( + getattr(self.Op, "forceflat", None) if forceflat is None else forceflat + ) name = getattr(Op, "name", None) if name is None else name if dtype is not None: @@ -138,6 +149,8 @@ def __init__( self.clinear = clinear if explicit is not None: self.explicit = explicit + if forceflat is not None: + self.forceflat = forceflat self.name = name # counters @@ -263,6 +276,21 @@ def explicit(self, new_explicit: bool) -> None: def explicit(self): del self._explicit + @property + def forceflat(self): + return getattr(self, "_forceflat", None) + + @forceflat.setter + def forceflat(self, new_forceflat: bool) -> None: + # note that this can also be None so we check before forcing bool + self._forceflat = ( + new_forceflat if new_forceflat is None else bool(new_forceflat) + ) + + @forceflat.deleter + def forceflat(self): + del self._forceflat + @property def name(self): return getattr(self, "_name", None) @@ -328,16 +356,32 @@ def __add__(self, x: LinearOperator) -> LinearOperator: Op, exclude=[ "explicit", + "forceflat", "name", ], ) Op.clinear = Op.clinear and Opx.clinear Op.explicit = False + if self.forceflat is None and Opx.forceflat is None: + Op.forceflat = None + elif self.forceflat is not None and Opx.forceflat is not None: + # Define forceflat only if differing, otherwise raise error + if self.forceflat != Opx.forceflat: + raise ValueError( + f"Operators have conflicting forceflat {Op.forceflat} != {Opx.forceflat}" + ) + Op.forceflat = self.forceflat + else: # Only one of them is None + Op.forceflat = ( + self.forceflat if self.forceflat is not None else Opx.forceflat + ) + # Replace if shape-like if len(self.dims) == 1: Op.dims = Opx.dims if len(self.dimsd) == 1: Op.dimsd = Opx.dimsd + return Op else: return NotImplemented @@ -374,7 +418,7 @@ def _copy_attributes( """Copy attributes from one LinearOperator to another""" if exclude is None: exclude = ["name"] - attrs = ["dims", "dimsd", "clinear", "explicit", "name"] + attrs = ["dims", "dimsd", "clinear", "explicit", "forceflat", "name"] if exclude is not None: for item in exclude: attrs.remove(item) @@ -582,9 +626,22 @@ def dot(self, x: NDArray) -> NDArray: # to allow mixing pylops and scipy operators) Opx = aslinearoperator(x) Op = _ProductLinearOperator(self, Opx) - self._copy_attributes(Op, exclude=["dims", "explicit", "name"]) + self._copy_attributes(Op, exclude=["dims", "explicit", "forceflat", "name"]) Op.clinear = Op.clinear and Opx.clinear Op.explicit = False + if self.forceflat is None and Opx.forceflat is None: + Op.forceflat = None + elif self.forceflat is not None and Opx.forceflat is not None: + # Define forceflat only if differing, otherwise raise error + if self.forceflat != Opx.forceflat: + raise ValueError( + f"Operators have conflicting forceflat {Op.forceflat} != {Opx.forceflat}" + ) + Op.forceflat = self.forceflat + else: # Only one of them is None + Op.forceflat = ( + self.forceflat if self.forceflat is not None else Opx.forceflat + ) Op.dims = Opx.dims return Op elif np.isscalar(x): @@ -608,17 +665,25 @@ def dot(self, x: NDArray) -> NDArray: if is_dims_shaped: # (dims1, ..., dimsK) => (dims1 * ... * dimsK,) == self.shape x = x.ravel() - if is_dims_shaped_matrix: + if is_dims_shaped_matrix and not self.forceflat: # (dims1, ..., dimsK, P) => (dims1 * ... * dimsK, P) x = x.reshape((-1, x.shape[-1])) if x.ndim == 1: y = self.matvec(x) - if is_dims_shaped and get_ndarray_multiplication(): + if ( + is_dims_shaped + and not self.forceflat + and get_ndarray_multiplication() + ): y = y.reshape(self.dimsd) return y elif x.ndim == 2: y = self.matmat(x) - if is_dims_shaped_matrix and get_ndarray_multiplication(): + if ( + is_dims_shaped_matrix + and not self.forceflat + and get_ndarray_multiplication() + ): y = y.reshape((*self.dimsd, -1)) return y else: @@ -776,7 +841,6 @@ def tosparse(self) -> NDArray: # loop through columns of self for i in range(n): - # make a unit vector for the ith column unit_i = np.zeros(n) unit_i[i] = 1 diff --git a/pylops/optimization/cls_leastsquares.py b/pylops/optimization/cls_leastsquares.py index 8382a2fc..a9e20fec 100644 --- a/pylops/optimization/cls_leastsquares.py +++ b/pylops/optimization/cls_leastsquares.py @@ -16,7 +16,6 @@ from pylops.optimization.basesolver import Solver from pylops.optimization.basic import cg, cgls from pylops.utils.backend import get_array_module -from pylops.utils.decorators import disable_ndarray_multiplication from pylops.utils.typing import NDArray if TYPE_CHECKING: @@ -89,7 +88,6 @@ def _print_finalize(self) -> None: print(f"\nTotal time (s) = {self.telapsed:.2f}") print("-" * 55 + "\n") - @disable_ndarray_multiplication def setup( self, y: NDArray, @@ -157,13 +155,13 @@ def setup( # normal equations if Weight is not None: - self.y_normal = self.OpH * Weight * y + self.y_normal = self.Op.rmatvec(Weight.matvec(y)) else: - self.y_normal = self.OpH * y + self.y_normal = self.Op.rmatvec(y) if Weight is not None: - self.Op_normal = self.OpH * Weight * self.Op + self.Op_normal = self.OpH @ Weight @ self.Op else: - self.Op_normal = self.OpH * self.Op + self.Op_normal = self.OpH @ self.Op # add regularization terms if epsI > 0: @@ -178,9 +176,8 @@ def setup( and self.dataregs is not None ): for epsR, Reg, datareg in zip(self.epsRs, self.Regs, self.dataregs): - self.RegH = Reg.H - self.y_normal += epsR**2 * self.RegH * datareg - self.Op_normal += epsR**2 * self.RegH * Reg + self.y_normal += epsR**2 * Reg.rmatvec(datareg) + self.Op_normal += epsR**2 * Reg.H @ Reg if epsNRs is not None and NRegs is not None: for epsNR, NReg in zip(epsNRs, NRegs): @@ -197,7 +194,6 @@ def step(self) -> None: "step method is not implemented. Use directly run or solve." ) - @disable_ndarray_multiplication def run( self, x: NDArray, @@ -240,7 +236,7 @@ def run( """ if x is not None: - self.y_normal = self.y_normal - self.Op_normal * x + self.y_normal = self.y_normal - self.Op_normal.matvec(x) if engine == "scipy" and self.ncp == np: if "atol" not in kwargs_solver: kwargs_solver["atol"] = "legacy" @@ -451,7 +447,6 @@ def _print_finalize(self) -> None: print(f"\nTotal time (s) = {self.telapsed:.2f}") print("-" * 65 + "\n") - @disable_ndarray_multiplication def setup( self, y: NDArray, @@ -504,10 +499,10 @@ def setup( self.RegOp: LinearOperator if Weight is not None: if Regs is None: - self.RegOp = Weight * self.Op + self.RegOp = Weight @ self.Op else: self.RegOp = RegularizedOperator( - Weight * self.Op, Regs, epsRs=self.epsRs + Weight @ self.Op, Regs, epsRs=self.epsRs ) else: if Regs is None: @@ -517,7 +512,7 @@ def setup( # augumented data if Weight is not None: - self.datatot: NDArray = Weight * self.y.copy() + self.datatot: NDArray = Weight.matvec(self.y.copy()) else: self.datatot = self.y.copy() @@ -537,7 +532,6 @@ def step(self) -> None: "step method is not implemented. Use directly run or solve." ) - @disable_ndarray_multiplication def run( self, x: NDArray, @@ -586,7 +580,7 @@ def run( """ if x is not None: - self.datatot = self.datatot - self.RegOp * x + self.datatot = self.datatot - self.RegOp.matvec(x) if engine == "scipy" and self.ncp == np: if show: kwargs_solver["show"] = 1 @@ -722,7 +716,6 @@ def _print_finalize(self) -> None: print(f"\nTotal time (s) = {self.telapsed:.2f}") print("-" * 65 + "\n") - @disable_ndarray_multiplication def setup( self, y: NDArray, @@ -746,7 +739,7 @@ def setup( self.ncp = get_array_module(y) # preconditioned operator - self.POp = self.Op * P + self.POp = self.Op @ P # print setup if show: @@ -759,7 +752,6 @@ def step(self) -> None: "step method is not implemented. Use directly run or solve." ) - @disable_ndarray_multiplication def run( self, x: NDArray, @@ -808,7 +800,7 @@ def run( """ if x is not None: - self.y = self.y - self.Op * x + self.y = self.y - self.Op.matvec(x) if engine == "scipy" and self.ncp == np: if show: kwargs_solver["show"] = 1 @@ -830,7 +822,7 @@ def run( pinv = pinv.ravel() else: raise NotImplementedError("Engine must be scipy or pylops") - xinv = self.P * pinv + xinv = self.P.matvec(pinv) if x is not None: xinv = x + xinv return xinv, istop, itn, r1norm, r2norm diff --git a/pylops/optimization/cls_sparsity.py b/pylops/optimization/cls_sparsity.py index 11f38aee..6ab4a30d 100644 --- a/pylops/optimization/cls_sparsity.py +++ b/pylops/optimization/cls_sparsity.py @@ -14,7 +14,6 @@ from pylops.optimization.leastsquares import regularized_inversion from pylops.utils import deps from pylops.utils.backend import get_array_module, get_module_name -from pylops.utils.decorators import disable_ndarray_multiplication from pylops.utils.typing import InputDimsLike, NDArray, SamplingLike spgl1_message = deps.spgl1_import("the spgl1 solver") @@ -422,18 +421,16 @@ def _step_model(self, x: NDArray, **kwargs_solver) -> NDArray: if self.iiter == 0: # first iteration (unweighted least-squares) if self.ncp == np: - x = ( - self.Op.H - @ lsqr( + x = self.Op.rmatvec( + lsqr( self.Op @ self.Op.H + (self.epsI**2) * self.Iop, self.y, **kwargs_solver, )[0] ) else: - x = ( - self.Op.H - @ cgls( + x = self.Op.rmatvec( + cgls( self.Op @ self.Op.H + (self.epsI**2) * self.Iop, self.y, self.ncp.zeros(int(self.Op.shape[0]), dtype=self.Op.dtype), @@ -446,25 +443,25 @@ def _step_model(self, x: NDArray, **kwargs_solver) -> NDArray: self.rw = self.rw / self.rw.max() R = Diagonal(self.rw, dtype=self.rw.dtype) if self.ncp == np: - x = ( - R - @ self.Op.H - @ lsqr( - self.Op @ R @ self.Op.H + self.epsI**2 * self.Iop, - self.y, - **kwargs_solver, - )[0] + x = R.matvec( + self.Op.rmatvec( + lsqr( + self.Op @ R @ self.Op.H + self.epsI**2 * self.Iop, + self.y, + **kwargs_solver, + )[0] + ) ) else: - x = ( - R - @ self.Op.H - @ cgls( - self.Op @ R @ self.Op.H + self.epsI**2 * self.Iop, - self.y, - self.ncp.zeros(int(self.Op.shape[0]), dtype=self.Op.dtype), - **kwargs_solver, - )[0] + x = R.matvec( + self.Op.rmatvec( + cgls( + self.Op @ R @ self.Op.H + self.epsI**2 * self.Iop, + self.y, + self.ncp.zeros(int(self.Op.shape[0]), dtype=self.Op.dtype), + **kwargs_solver, + )[0] + ) ) return x @@ -494,7 +491,7 @@ def step(self, x: NDArray, show: bool = False, **kwargs_solver) -> NDArray: x = self._step(x, **kwargs_solver) # compute residual - self.r: NDArray = self.y - self.Op * x + self.r: NDArray = self.y - self.Op.matvec(x) self.rnorm = self.ncp.linalg.norm(self.r) self.iiter += 1 @@ -540,7 +537,7 @@ def run( nouter = nouter if self.nouter is None else self.nouter if x is not None: self.x0 = x.copy() - self.y = self.y - self.Op * x + self.y = self.y - self.Op.matvec(x) # choose xold to ensure tolerance test is passed initially xold = x.copy() + np.inf while self.iiter < nouter and self.ncp.linalg.norm(x - xold) >= self.tolIRLS: @@ -1140,7 +1137,8 @@ def setup( Parameters ---------- y : :obj:`np.ndarray` - Data of size :math:`[N \times 1]` + Data of size :math:`[N \times 1]` or :math:`[N \times R]` where + a solution for multiple right-hand-side is found when ``R>1``. x0: :obj:`numpy.ndarray`, optional Initial guess niter : :obj:`int` @@ -1194,6 +1192,21 @@ def setup( self.ncp = get_array_module(y) + # choose matvec/rmatvec or matmat/rmatmat based on R + if y.ndim > 1 and y.shape[1] > 1: + self.Opmatvec = self.Op.matmat + self.Oprmatvec = self.Op.rmatmat + if self.SOp is not None: + self.SOpmatvec = self.SOp.matmat + self.SOprmatvec = self.SOp.rmatmat + else: + self.Opmatvec = self.Op.matvec + self.Oprmatvec = self.Op.rmatvec + if self.SOp is not None: + self.SOpmatvec = self.SOp.matvec + self.SOprmatvec = self.SOp.rmatvec + + # choose thresholding function if threshkind not in [ "hard", "soft", @@ -1216,7 +1229,6 @@ def setup( "soft-percentile, or half-percentile thresholding" ) - # choose thresholding function self.threshf: Callable[[NDArray, float], NDArray] if threshkind == "soft": self.threshf = _softthreshold @@ -1240,7 +1252,7 @@ def setup( self.alpha = alpha elif not hasattr(self, "alpha"): # compute largest eigenvalues of Op^H * Op - Op1 = self.Op.H * self.Op + Op1 = self.Op.H @ self.Op if get_module_name(self.ncp) == "numpy": maxeig: float = np.abs( Op1.eigs( @@ -1318,7 +1330,7 @@ def step(self, x: NDArray, show: bool = False) -> Tuple[NDArray, float]: # store old vector xold = x.copy() # compute residual - res: NDArray = self.y - self.Op @ x + res: NDArray = self.y - self.Opmatvec(x) if self.monitorres: self.normres = np.linalg.norm(res) if self.normres > self.normresold: @@ -1331,18 +1343,18 @@ def step(self, x: NDArray, show: bool = False) -> Tuple[NDArray, float]: self.normresold = self.normres # compute gradient - grad: NDArray = self.alpha * (self.Op.H @ res) + grad: NDArray = self.alpha * (self.Oprmatvec(res)) # update inverted model x_unthesh: NDArray = x + grad if self.SOp is not None: - x_unthesh = self.SOp.H @ x_unthesh + x_unthesh = self.SOprmatvec(x_unthesh) if self.perc is None and self.decay is not None: x = self.threshf(x_unthesh, self.decay[self.iiter] * self.thresh) elif self.perc is not None: x = self.threshf(x_unthesh, 100 - self.perc) if self.SOp is not None: - x = self.SOp @ x + x = self.SOpmatvec(x) # model update xupdate = np.linalg.norm(x - xold) @@ -1578,7 +1590,7 @@ def step(self, x: NDArray, z: NDArray, show: bool = False) -> NDArray: ---------- x : :obj:`np.ndarray` Current model vector to be updated by a step of ISTA - x : :obj:`np.ndarray` + z : :obj:`np.ndarray` Current auxiliary model vector to be updated by a step of ISTA show : :obj:`bool`, optional Display iteration log @@ -1596,7 +1608,7 @@ def step(self, x: NDArray, z: NDArray, show: bool = False) -> NDArray: # store old vector xold = x.copy() # compute residual - resz: NDArray = self.y - self.Op @ z + resz: NDArray = self.y - self.Opmatvec(z) if self.monitorres: self.normres = np.linalg.norm(resz) if self.normres > self.normresold: @@ -1609,18 +1621,18 @@ def step(self, x: NDArray, z: NDArray, show: bool = False) -> NDArray: self.normresold = self.normres # compute gradient - grad: NDArray = self.alpha * (self.Op.H @ resz) + grad: NDArray = self.alpha * (self.Oprmatvec(resz)) # update inverted model x_unthesh: NDArray = z + grad if self.SOp is not None: - x_unthesh = self.SOp.H @ x_unthesh + x_unthesh = self.SOprmatvec(x_unthesh) if self.perc is None and self.decay is not None: x = self.threshf(x_unthesh, self.decay[self.iiter] * self.thresh) elif self.perc is not None: x = self.threshf(x_unthesh, 100 - self.perc) if self.SOp is not None: - x = self.SOp @ x + x = self.SOpmatvec(x) # update auxiliary coefficients told = self.t @@ -1753,7 +1765,6 @@ def _print_finalize(self) -> None: print(f"\nTotal time (s) = {self.telapsed:.2f}") print("-" * 80 + "\n") - @disable_ndarray_multiplication def setup( self, y: NDArray, @@ -1800,7 +1811,6 @@ def step(self) -> None: "step method is not implemented. Use directly run or solve." ) - @disable_ndarray_multiplication def run( self, x: NDArray, @@ -1871,7 +1881,7 @@ def run( """ pinv, _, _, info = ext_spgl1( - self.Op if self.SOp is None else self.Op * self.SOp.H, + self.Op if self.SOp is None else self.Op @ self.SOp.H, self.y, tau=self.tau, sigma=self.sigma, @@ -1879,7 +1889,7 @@ def run( **kwargs_spgl1, ) - xinv = pinv.copy() if self.SOp is None else self.SOp.H * pinv + xinv = pinv.copy() if self.SOp is None else self.SOp.rmatvec(pinv) return xinv, pinv, info def solve( @@ -2238,13 +2248,15 @@ def step( )[0] # Shrinkage self.d = [ - _softthreshold(self.RegsL1[ireg] * x + self.b[ireg], self.epsRL1s[ireg]) + _softthreshold( + self.RegsL1[ireg].matvec(x) + self.b[ireg], self.epsRL1s[ireg] + ) for ireg in range(self.nregsL1) ] # Bregman update self.b = [ - self.b[ireg] + self.tau * (self.RegsL1[ireg] * x - self.d[ireg]) + self.b[ireg] + self.tau * (self.RegsL1[ireg].matvec(x) - self.d[ireg]) for ireg in range(self.nregsL1) ] diff --git a/pylops/signalprocessing/__init__.py b/pylops/signalprocessing/__init__.py index 2ce1fed1..8efa532e 100755 --- a/pylops/signalprocessing/__init__.py +++ b/pylops/signalprocessing/__init__.py @@ -12,6 +12,7 @@ ConvolveND ND convolution operator. NonStationaryConvolve1D 1D nonstationary convolution operator. NonStationaryConvolve2D 2D nonstationary convolution operator. + NonStationaryConvolve3D 3D nonstationary convolution operator. NonStationaryFilters1D 1D nonstationary filter estimation operator. NonStationaryFilters2D 2D nonstationary filter estimation operator. Interp Interpolation operator. @@ -43,6 +44,7 @@ from .convolve2d import * from .nonstatconvolve1d import * from .nonstatconvolve2d import * +from .nonstatconvolve3d import * from .shift import * from .interp import * from .bilinear import * @@ -71,6 +73,7 @@ "Convolve2D", "NonStationaryConvolve1D", "NonStationaryConvolve2D", + "NonStationaryConvolve3D", "NonStationaryFilters1D", "NonStationaryFilters2D", "Interp", diff --git a/pylops/signalprocessing/_nonstatconvolve3d_cuda.py b/pylops/signalprocessing/_nonstatconvolve3d_cuda.py new file mode 100644 index 00000000..5b309050 --- /dev/null +++ b/pylops/signalprocessing/_nonstatconvolve3d_cuda.py @@ -0,0 +1,154 @@ +from math import floor + +from numba import cuda + + +@cuda.jit(max_registers=40) +def _matvec_rmatvec( + x, y, hs, hshape, xdims, ohx, ohy, ohz, dhx, dhy, dhz, nhx, nhy, nhz, rmatvec +): + """Cuda kernels for NonStationaryConvolve3D operator + + Cuda implementation of matvec and rmatvec for NonStationaryConvolve3D operator. See + :class:`pylops.signalprocessing.NonStationaryConvolve3D` for details about input parameters. + + """ + ix, iy, iz = cuda.grid(3) + + if ix < xdims[0] and iy < xdims[1] and iz < xdims[2]: + + # find closest filters and interpolate h + ihx_l = int(floor((ix - ohx) / dhx)) # id number of left for hs_arr + ihy_b = int(floor((iy - ohy) / dhy)) # id number of back for hs_arr + ihz_t = int(floor((iz - ohz) / dhz)) # id number of top for hs_arr + + dhx_r = (ix - ohx) / dhx - ihx_l # weight for right psfs, left 1-dhx_r + dhy_f = (iy - ohy) / dhy - ihy_b # weight for front psfs, back 1-dhy_f + dhz_d = (iz - ohz) / dhz - ihz_t # weight for down psfs, top 1-dhz_d + + if ihx_l < 0: + ihx_l = ihx_r = 0 + dhx_l = dhx_r = 0.5 + elif ihx_l >= nhx - 1: + ihx_l = ihx_r = nhx - 1 + dhx_l = dhx_r = 0.5 + else: + ihx_r = ihx_l + 1 + dhx_l = 1.0 - dhx_r + + if ihy_b < 0: + ihy_b = ihy_f = 0 + dhy_b = dhy_f = 0.5 + elif ihy_b >= nhy - 1: + ihy_b = ihy_f = nhy - 1 + dhy_b = dhy_f = 0.5 + else: + ihy_f = ihy_b + 1 + dhy_b = 1.0 - dhy_f + + if ihz_t < 0: + ihz_t = ihz_d = 0 + dhz_t = dhz_d = 0.5 + elif ihz_t >= nhz - 1: + ihz_t = ihz_d = nhz - 1 + dhz_t = dhz_d = 0.5 + else: + ihz_d = ihz_t + 1 + dhz_t = 1.0 - dhz_d + + h_lbt = hs[ihx_l, ihy_b, ihz_t] + h_lbd = hs[ihx_l, ihy_b, ihz_d] + h_lft = hs[ihx_l, ihy_f, ihz_t] + h_lfd = hs[ihx_l, ihy_f, ihz_d] + + h_rbt = hs[ihx_r, ihy_b, ihz_t] + h_rbd = hs[ihx_r, ihy_b, ihz_d] + h_rft = hs[ihx_r, ihy_f, ihz_t] + h_rfd = hs[ihx_r, ihy_f, ihz_d] + + # find extremes of model where to apply h (in case h is going out of model) + xextremes = ( + max(0, ix - hshape[0] // 2), + min(ix + hshape[0] // 2 + 1, xdims[0]), + ) + yextremes = ( + max(0, iy - hshape[1] // 2), + min(iy + hshape[1] // 2 + 1, xdims[1]), + ) + zextremes = ( + max(0, iz - hshape[2] // 2), + min(iz + hshape[2] // 2 + 1, xdims[2]), + ) + + # find extremes of h (in case h is going out of model) + hxextremes = ( + max(0, -ix + hshape[0] // 2), + min(hshape[0], hshape[0] // 2 + (xdims[0] - ix)), + ) + hyextremes = ( + max(0, -iy + hshape[1] // 2), + min(hshape[1], hshape[1] // 2 + (xdims[1] - iy)), + ) + hzextremes = ( + max(0, -iz + hshape[2] // 2), + min(hshape[2], hshape[2] // 2 + (xdims[2] - iz)), + ) + + # place filter in output + for ixx, hxx in zip( + range(xextremes[0], xextremes[1]), range(hxextremes[0], hxextremes[1]) + ): + for iyy, hyy in zip( + range(yextremes[0], yextremes[1]), range(hyextremes[0], hyextremes[1]) + ): + for izz, hzz in zip( + range(zextremes[0], zextremes[1]), + range(hzextremes[0], hzextremes[1]), + ): + h = ( + dhx_l * dhy_b * dhz_t * h_lbt[hxx, hyy, hzz] + + dhx_l * dhy_b * dhz_d * h_lbd[hxx, hyy, hzz] + + dhx_l * dhy_f * dhz_t * h_lft[hxx, hyy, hzz] + + dhx_l * dhy_f * dhz_d * h_lfd[hxx, hyy, hzz] + + dhx_r * dhy_b * dhz_t * h_rbt[hxx, hyy, hzz] + + dhx_r * dhy_b * dhz_d * h_rbd[hxx, hyy, hzz] + + dhx_r * dhy_f * dhz_t * h_rft[hxx, hyy, hzz] + + dhx_r * dhy_f * dhz_d * h_rfd[hxx, hyy, hzz] + ) + + if rmatvec: + y[ix, iy, iz] += h * x[ixx, iyy, izz] + else: + cuda.atomic.add(y, (ixx, iyy, izz), x[ix, iy, iz] * h) + + +def _matvec_rmatvec_call( + x, + y, + hs, + hshape, + xdims, + ohx, + ohy, + ohz, + dhx, + dhy, + dhz, + nhx, + nhy, + nhz, + rmatvec=False, + num_blocks=(1, 32, 16), + num_threads_per_blocks=(24, 24, 24), +): + """Caller for NonStationaryConvolve3D operator + + Caller for cuda implementation of matvec and rmatvec for NonStationaryConvolve3D operator, with same signature + as numpy/numba counterparts. See :class:`pylops.signalprocessing.NonStationaryConvolve3D` for details about + input parameters. + + """ + _matvec_rmatvec[num_blocks, num_threads_per_blocks]( + x, y, hs, hshape, xdims, ohx, ohy, ohz, dhx, dhy, dhz, nhx, nhy, nhz, rmatvec + ) + return y diff --git a/pylops/signalprocessing/bilinear.py b/pylops/signalprocessing/bilinear.py index 7924e602..6bc62e58 100644 --- a/pylops/signalprocessing/bilinear.py +++ b/pylops/signalprocessing/bilinear.py @@ -37,6 +37,12 @@ class Bilinear(LinearOperator): for interpolation. dims : :obj:`list` Number of samples for each dimension + forceflat : :obj:`bool`, optional + .. versionadded:: 2.2.0 + + Force an array to be flattened after rmatvec. Note that this is only + required when `len(dims)=2`, otherwise pylops will detect whether to + return a 1d or nd array. dtype : :obj:`str`, optional Type of elements in input array. name : :obj:`str`, optional @@ -90,13 +96,30 @@ def __init__( self, iava: IntNDArray, dims: InputDimsLike, + forceflat: bool = None, dtype: DTypeLike = "float64", name: str = "B", ) -> None: # define dimension of data ndims = len(dims) dimsd = [len(iava[1])] + list(dims[2:]) - super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dimsd, name=name) + # check if forceflat is needed and set it back to None otherwise + if ndims > 2: + if forceflat is not None: + logging.warning( + f"setting forceflat=None since len(dims)={len(dims)}>2. " + f"PyLops will automatically detect whether to return " + f"a 1d or nd array based on the shape of the input" + f"array." + ) + forceflat = None + super().__init__( + dtype=np.dtype(dtype), + dims=dims, + dimsd=dimsd, + forceflat=forceflat, + name=name, + ) ncp = get_array_module(iava) # check non-unique pairs (works only with numpy arrays) diff --git a/pylops/signalprocessing/convolve1d.py b/pylops/signalprocessing/convolve1d.py index 9921ccf9..74057aa8 100644 --- a/pylops/signalprocessing/convolve1d.py +++ b/pylops/signalprocessing/convolve1d.py @@ -1,5 +1,6 @@ __all__ = ["Convolve1D"] +from functools import partial from typing import Callable, Tuple, Union import numpy as np @@ -8,6 +9,7 @@ from pylops import LinearOperator from pylops.utils._internal import _value_or_sized_to_tuple from pylops.utils.backend import ( + get_array_module, get_convolve, get_fftconvolve, get_oaconvolve, @@ -18,7 +20,10 @@ def _choose_convfunc( - x: npt.ArrayLike, method: Union[None, str], dims + x: npt.ArrayLike, + method: Union[None, str], + dims: Union[int, InputDimsLike], + axis: int = -1, ) -> Tuple[Callable, str]: """Choose convolution function @@ -34,28 +39,184 @@ def _choose_convfunc( if method is None: method = "fft" if method == "fft": - convfunc = get_fftconvolve(x) + convfunc = partial(get_fftconvolve(x), axes=axis) elif method == "overlapadd": - convfunc = get_oaconvolve(x) + convfunc = partial(get_oaconvolve(x), axes=axis)(x) else: raise NotImplementedError("method must be fft or overlapadd") return convfunc, method +def _pad_along_axis(array: np.ndarray, pad_size: tuple, axis: int = 0) -> np.ndarray: + + npad = [(0, 0)] * array.ndim + npad[axis] = pad_size + return np.pad(array, pad_width=npad) + + +class _Convolve1Dshort(LinearOperator): + r"""1D convolution operator with compact filter (shorter than input signal)""" + + def __init__( + self, + dims: Union[int, InputDimsLike], + h: NDArray, + offset: int = 0, + axis: int = -1, + method: str = None, + dtype: DTypeLike = "float64", + name: str = "C", + ) -> None: + dims = _value_or_sized_to_tuple(dims) + super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dims, name=name) + self.axis = axis + self.nh = h.size if h.ndim == 1 else h.shape[axis] + if offset > self.nh - 1: + raise ValueError("offset must be smaller than h.shape[axis] - 1") + self.h = h + self.offset = 2 * (self.nh // 2 - int(offset)) + if self.nh % 2 == 0: + self.offset -= 1 + if self.offset != 0: + self.h = _pad_along_axis( + self.h, + (max(self.offset, 0), -min(self.offset, 0)), + axis=-1 if h.ndim == 1 else axis, + ) + self.hstar = np.flip(self.h, axis=-1) + + # add dimensions to filter to match dimensions of model and data + if self.h.ndim == 1: + hdims = np.ones(len(self.dims), dtype=int) + hdims[self.axis] = len(self.h) + self.h = self.h.reshape(hdims) + self.hstar = self.hstar.reshape(hdims) + + # choose method and function handle + self.convfunc, self.method = _choose_convfunc(h, method, self.dims, self.axis) + + @reshaped + def _matvec(self, x: NDArray) -> NDArray: + if type(self.h) is not type(x): + self.h = to_cupy_conditional(x, self.h) + self.convfunc, self.method = _choose_convfunc( + self.h, self.method, self.dims + ) + return self.convfunc(x, self.h, mode="same") + + @reshaped + def _rmatvec(self, x: NDArray) -> NDArray: + if type(self.hstar) is not type(x): + self.hstar = to_cupy_conditional(x, self.hstar) + self.convfunc, self.method = _choose_convfunc( + self.hstar, self.method, self.dims + ) + return self.convfunc(x, self.hstar, mode="same") + + +class _Convolve1Dlong(LinearOperator): + """1D convolution operator with extended filter (larger than input signal)""" + + def __init__( + self, + dims: Union[int, InputDimsLike], + h: NDArray, + offset: int = 0, + axis: int = -1, + method: str = None, + dtype: DTypeLike = "float64", + name: str = "C", + ) -> None: + dims = _value_or_sized_to_tuple(dims) + dimsd = h.shape + super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dimsd, name=name) + + # create filter + self.axis = axis + if offset > self.dims[self.axis] - 1: + raise ValueError("offset must be smaller than self.dims[self.axis] - 1") + self.nh = h.size if h.ndim == 1 else h.shape[axis] + self.h = h + self.offset = 2 * (self.dims[self.axis] // 2 - int(offset)) + if self.dims[self.axis] % 2 == 0: + self.offset -= 1 + self.hstar = np.flip(self.h, axis=-1) + + self.pad = np.zeros((len(dims), 2), dtype=int) + self.pad[self.axis, 0] = max(self.offset, 0) + self.pad[self.axis, 1] = -min(self.offset, 0) + + self.padd = np.zeros((len(dims), 2), dtype=int) + self.padd[self.axis, 1] = max(self.offset, 0) + self.padd[self.axis, 0] = -min(self.offset, 0) + + # add dimensions to filter to match dimensions of model and data + if self.h.ndim == 1: + hdims = np.ones(len(self.dims), dtype=int) + hdims[self.axis] = len(self.h) + self.h = self.h.reshape(hdims) + self.hstar = self.hstar.reshape(hdims) + + # choose method and function handle + self.convfunc, self.method = _choose_convfunc(h, method, self.dims, self.axis) + + @reshaped + def _matvec(self, x: NDArray) -> NDArray: + if type(self.h) is not type(x): + self.h = to_cupy_conditional(x, self.h) + self.convfunc, self.method = _choose_convfunc( + self.h, self.method, self.dims + ) + x = np.pad(x, self.pad) + y = self.convfunc(self.h, x, mode="same") + return y + + @reshaped + def _rmatvec(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + if type(self.h) is not type(x): + self.hstar = to_cupy_conditional(x, self.hstar) + self.convfunc, self.method = _choose_convfunc( + self.hstar, self.method, self.dims + ) + x = np.pad(x, self.padd) + y = self.convfunc(self.hstar, x) + if self.dims[self.axis] % 2 == 0: + y = ncp.take( + y, + range( + len(y) // 2 - self.dims[self.axis] // 2, + len(y) // 2 + self.dims[self.axis] // 2, + ), + axis=self.axis, + ) + else: + y = ncp.take( + y, + range( + len(y) // 2 - self.dims[self.axis] // 2, + len(y) // 2 + self.dims[self.axis] // 2 + 1, + ), + axis=self.axis, + ) + return y + + class Convolve1D(LinearOperator): r"""1D convolution operator. - Apply one-dimensional convolution with a compact filter to model (and data) - along an ``axis`` of a multi-dimensional array. + Apply one-dimensional convolution with i) a compact filter (shorter than input signal) or + ii) an extended filter (larger than input signal) to model (and data) along an ``axis`` + of a multi-dimensional array. Parameters ---------- dims : :obj:`list` or :obj:`int` - Number of samples for each dimension + Number of samples for each dimension of the model h : :obj:`numpy.ndarray` - 1d compact filter to be convolved to input signal + 1d filter to be convolved to input signal offset : :obj:`int` - Index of the center of the compact filter + Index of the center of the filter axis : :obj:`int`, optional .. versionadded:: 2.0.0 @@ -90,7 +251,7 @@ class Convolve1D(LinearOperator): Notes ----- The Convolve1D operator applies convolution between the input signal - :math:`x(t)` and a compact filter kernel :math:`h(t)` in forward model: + :math:`x(t)` and a filter kernel :math:`h(t)` in forward model: .. math:: y(t) = \int\limits_{-\infty}^{\infty} h(t-\tau) x(\tau) \,\mathrm{d}\tau @@ -137,52 +298,23 @@ def __init__( dtype: DTypeLike = "float64", name: str = "C", ) -> None: - dims = _value_or_sized_to_tuple(dims) - super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dims, name=name) - - self.axis = axis - if offset > len(h) - 1: - raise ValueError("offset must be smaller than len(h) - 1") - self.h = h - self.hstar = np.flip(self.h, axis=-1) - self.nh = len(h) - self.offset = 2 * (self.nh // 2 - int(offset)) - if self.nh % 2 == 0: - self.offset -= 1 - if self.offset != 0: - self.h = np.pad( - self.h, - ( - self.offset if self.offset > 0 else 0, - -self.offset if self.offset < 0 else 0, - ), - mode="constant", - ) - self.hstar = np.flip(self.h, axis=-1) - - # add dimensions to filter to match dimensions of model and data - hdims = np.ones(len(self.dims), dtype=int) - hdims[self.axis] = len(self.h) - self.h = self.h.reshape(hdims) - self.hstar = self.hstar.reshape(hdims) - - # choose method and function handle - self.convfunc, self.method = _choose_convfunc(h, method, self.dims) + nh = h.size if h.ndim == 1 else h.shape[axis] + if nh <= _value_or_sized_to_tuple(dims)[axis]: + convop = _Convolve1Dshort + else: + convop = _Convolve1Dlong + Op = convop( + dims=dims, + h=h, + offset=offset, + axis=axis, + method=method, + dtype=dtype, + ) + super().__init__(Op=Op, name=name) - @reshaped def _matvec(self, x: NDArray) -> NDArray: - if type(self.h) is not type(x): - self.h = to_cupy_conditional(x, self.h) - self.convfunc, self.method = _choose_convfunc( - self.h, self.method, self.dims - ) - return self.convfunc(x, self.h, mode="same") + return super()._matvec(x) - @reshaped def _rmatvec(self, x: NDArray) -> NDArray: - if type(self.hstar) is not type(x): - self.hstar = to_cupy_conditional(x, self.hstar) - self.convfunc, self.method = _choose_convfunc( - self.hstar, self.method, self.dims - ) - return self.convfunc(x, self.hstar, mode="same") + return super()._rmatvec(x) diff --git a/pylops/signalprocessing/nonstatconvolve2d.py b/pylops/signalprocessing/nonstatconvolve2d.py index 569736d6..a7056986 100644 --- a/pylops/signalprocessing/nonstatconvolve2d.py +++ b/pylops/signalprocessing/nonstatconvolve2d.py @@ -36,15 +36,15 @@ class NonStationaryConvolve2D(LinearOperator): Apply non-stationary two-dimensional convolution. A varying compact filter is provided on a coarser grid and on-the-fly interpolation is applied in forward and adjoint modes. - Both input and output have size :math`n_x \time n_z`. + Both input and output have size :math:`n_x \times n_z`. Parameters ---------- dims : :obj:`list` or :obj:`int` - Number of samples for each dimension (which we refer to as :math`n_x \time n_z`). + Number of samples for each dimension (which we refer to as :math:`n_x \times n_z`). hs : :obj:`numpy.ndarray` Bank of 2d compact filters of size - :math:`n_{\text{filts},x} \times n_{\text{filts},z} \times n_h \times n_{h,x} \times n_{h,z}`. + :math:`n_{\text{filts},x} \times n_{\text{filts},z} \times n_{h,x} \times n_{h,z}`. Filters must have odd number of samples and are assumed to be centered in the middle of the filter support. ihx : :obj:`tuple` @@ -424,7 +424,7 @@ def _register_multiplications(self, engine: str) -> None: self._mv = self.__matvec self._rmv = self.__rmatvec - # use same matvec method as inNonStationaryConvolve2D + # use same matvec method as in NonStationaryConvolve2D __matvec = staticmethod(NonStationaryConvolve2D._matvec_rmatvec) @staticmethod diff --git a/pylops/signalprocessing/nonstatconvolve3d.py b/pylops/signalprocessing/nonstatconvolve3d.py new file mode 100644 index 00000000..236851f8 --- /dev/null +++ b/pylops/signalprocessing/nonstatconvolve3d.py @@ -0,0 +1,355 @@ +__all__ = ["NonStationaryConvolve3D"] + +import os +from typing import Tuple, Union + +import numpy as np + +from pylops import LinearOperator +from pylops.utils import deps +from pylops.utils._internal import _value_or_sized_to_tuple +from pylops.utils.backend import get_array_module +from pylops.utils.decorators import reshaped +from pylops.utils.typing import DTypeLike, InputDimsLike, NDArray + +jit_message = deps.numba_import("the nonstatconvolve3d module") + +if jit_message is None: + from numba import jit, prange + + from ._nonstatconvolve3d_cuda import ( + _matvec_rmatvec_call as _matvec_rmatvec_cuda_call, + ) + + # detect whether to use parallel or not + numba_threads = int(os.getenv("NUMBA_NUM_THREADS", "1")) + parallel = True if numba_threads != 1 else False +else: + prange = range + + +class NonStationaryConvolve3D(LinearOperator): + r"""3D non-stationary convolution operator. + + Apply non-stationary three-dimensional convolution. A varying compact filter + is provided on a coarser grid and on-the-fly interpolation is applied + in forward and adjoint modes. Both input and output have size :math:`n_x \times n_y \times n_z`. + + Parameters + ---------- + dims : :obj:`list` or :obj:`int` + Number of samples for each dimension (which we refer to as :math:`n_x \times n_y \times n_z`). + hs : :obj:`numpy.ndarray` + Bank of 3d compact filters of size + :math:`n_{\text{filts},x} \times n_{\text{filts},y} \times + n_{\text{filts},z} \times n_{h,x} \times n_{h,y} \times n_{h,z}`. + Filters must have odd number of samples and are assumed to be + centered in the middle of the filter support. + ihx : :obj:`tuple` + Indices of the x locations of the filters ``hs`` in the model (and data). Note + that the filters must be regularly sampled, i.e. :math:`dh_x=\text{diff}(ihx)=\text{const.}` + ihy : :obj:`tuple` + Indices of the y locations of the filters ``hs`` in the model (and data). Note + that the filters must be regularly sampled, i.e. :math:`dh_y=\text{diff}(ihy)=\text{const.}` + ihz : :obj:`tuple` + Indices of the z locations of the filters ``hs`` in the model (and data). Note + that the filters must be regularly sampled, i.e. :math:`dh_z=\text{diff}(ihz)=\text{const.}` + engine : :obj:`str`, optional + Engine used for spread computation (``numpy``, ``numba``, or ``cuda``) + num_threads_per_blocks : :obj:`tuple`, optional + Number of threads in each block (only when ``engine=cuda``) + dtype : :obj:`str`, optional + Type of elements in input array. + name : :obj:`str`, optional + Name of operator (to be used by :func:`pylops.utils.describe.describe`) + + Attributes + ---------- + shape : :obj:`tuple` + Operator shape + explicit : :obj:`bool` + Operator contains a matrix that can be solved explicitly (``True``) or + not (``False``) + + Raises + ------ + ValueError + If filters ``hs`` have even size + ValueError + If ``ihx``, ``ihy`` or ``ihz`` is not regularly sampled + NotImplementedError + If ``engine`` is neither ``numpy``, ``fftw``, nor ``scipy``. + + Notes + ----- + See :class:`pylops.signalprocessing.NonStationaryConvolve2D`. + + """ + + def __init__( + self, + dims: Union[int, InputDimsLike], + hs: NDArray, + ihx: InputDimsLike, + ihy: InputDimsLike, + ihz: InputDimsLike, + engine: str = "numpy", + num_threads_per_blocks: Tuple[int, int, int] = (2, 16, 16), + dtype: DTypeLike = "float64", + name: str = "C", + ) -> None: + if engine not in ["numpy", "numba", "cuda"]: + raise NotImplementedError("engine must be numpy or numba or cuda") + if hs.shape[3] % 2 == 0 or hs.shape[4] % 2 == 0 or hs.shape[5] % 2 == 0: + raise ValueError("filters hs must have odd length") + if ( + len(np.unique(np.diff(ihx))) > 1 + or len(np.unique(np.diff(ihy))) > 1 + or len(np.unique(np.diff(ihz))) > 1 + ): + raise ValueError( + "the indices of filters 'ih' are must be regularly sampled" + ) + + if ( + min(ihx) < 0 + or min(ihy) < 0 + or min(ihz) < 0 + or max(ihx) >= dims[0] + or max(ihy) >= dims[1] + or max(ihz) >= dims[2] + ): + raise ValueError( + "the indices of filters 'ih' must be larger than 0 and smaller than `dims`" + ) + self.hs = hs + self.hshape = hs.shape[3:] + self.ohx, self.dhx, self.nhx = ihx[0], ihx[1] - ihx[0], len(ihx) + self.ohy, self.dhy, self.nhy = ihy[0], ihy[1] - ihy[0], len(ihy) + self.ohz, self.dhz, self.nhz = ihz[0], ihz[1] - ihz[0], len(ihz) + self.ehx, self.ehx, self.ehz = ihx[-1], ihy[-1], ihz[-1] + self.dims = _value_or_sized_to_tuple(dims) + self.engine = engine + super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dims, name=name) + + # create additional input parameters for engine=cuda + self.kwargs_cuda = {} + if engine == "cuda": + self.kwargs_cuda["num_threads_per_blocks"] = num_threads_per_blocks + num_blocks_x = ( + self.dims[0] + num_threads_per_blocks[0] - 1 + ) // num_threads_per_blocks[0] + num_blocks_y = ( + self.dims[1] + num_threads_per_blocks[1] - 1 + ) // num_threads_per_blocks[1] + num_blocks_z = ( + self.dims[2] + num_threads_per_blocks[2] - 1 + ) // num_threads_per_blocks[2] + self.kwargs_cuda["num_blocks"] = (num_blocks_x, num_blocks_y, num_blocks_z) + self._register_multiplications(engine) + + def _register_multiplications(self, engine: str) -> None: + if engine == "numba": + numba_opts = dict(nopython=True, fastmath=True, nogil=True, parallel=True) + self._mvrmv = jit(**numba_opts)(self._matvec_rmatvec) + elif engine == "cuda": + self._mvrmv = _matvec_rmatvec_cuda_call + else: + self._mvrmv = self._matvec_rmatvec + + @staticmethod + def _matvec_rmatvec( + x: NDArray, + y: NDArray, + hs: NDArray, + hshape: Tuple[int, int, int], + xdims: Tuple[int, int, int], + ohx: float, + ohy: float, + ohz: float, + dhx: float, + dhy: float, + dhz: float, + nhx: int, + nhy: int, + nhz: int, + rmatvec: bool = False, + ) -> NDArray: + for ix in prange(xdims[0]): + for iy in range(xdims[1]): + for iz in range(xdims[2]): + # find closest filters and interpolate h + ihx_l = int( + np.floor((ix - ohx) / dhx) + ) # id number of left for hs_arr + ihy_b = int( + np.floor((iy - ohy) / dhy) + ) # id number of back for hs_arr + ihz_t = int( + np.floor((iz - ohz) / dhz) + ) # id number of top for hs_arr + + dhx_r = ( + ix - ohx + ) / dhx - ihx_l # weight for right psfs, left 1-ihz_t + dhy_f = ( + iy - ohy + ) / dhy - ihy_b # weight for front psfs, left 1-ihz_t + dhz_d = ( + iz - ohz + ) / dhz - ihz_t # weight for down psfs, top 1-dhz_d + + if ihx_l < 0: + ihx_l = ihx_r = 0 + dhx_l = dhx_r = 0.5 + elif ihx_l >= nhx - 1: + ihx_l = ihx_r = nhx - 1 + dhx_l = dhx_r = 0.5 + else: + ihx_r = ihx_l + 1 + dhx_l = 1.0 - dhx_r + + if ihy_b < 0: + ihy_b = ihy_f = 0 + dhy_b = dhy_f = 0.5 + elif ihy_b >= nhy - 1: + ihy_b = ihy_f = nhy - 1 + dhy_b = dhy_f = 0.5 + else: + ihy_f = ihy_b + 1 + dhy_b = 1.0 - dhy_f + + if ihz_t < 0: + ihz_t = ihz_d = 0 + dhz_t = dhz_d = 0.5 + elif ihz_t >= nhz - 1: + ihz_t = ihz_d = nhz - 1 + dhz_t = dhz_d = 0.5 + else: + ihz_d = ihz_t + 1 + dhz_t = 1.0 - dhz_d + + h_lbt = hs[ihx_l, ihy_b, ihz_t] + h_lbd = hs[ihx_l, ihy_b, ihz_d] + h_lft = hs[ihx_l, ihy_f, ihz_t] + h_lfd = hs[ihx_l, ihy_f, ihz_d] + + h_rbt = hs[ihx_r, ihy_b, ihz_t] + h_rbd = hs[ihx_r, ihy_b, ihz_d] + h_rft = hs[ihx_r, ihy_f, ihz_t] + h_rfd = hs[ihx_r, ihy_f, ihz_d] + + h = ( + dhx_l * dhy_b * dhz_t * h_lbt + + dhx_l * dhy_b * dhz_d * h_lbd + + dhx_l * dhy_f * dhz_t * h_lft + + dhx_l * dhy_f * dhz_d * h_lfd + + dhx_r * dhy_b * dhz_t * h_rbt + + dhx_r * dhy_b * dhz_d * h_rbd + + dhx_r * dhy_f * dhz_t * h_rft + + dhx_r * dhy_f * dhz_d * h_rfd + ) + + # find extremes of model where to apply h (in case h is going out of model) + xextremes = ( + max(0, ix - hshape[0] // 2), + min(ix + hshape[0] // 2 + 1, xdims[0]), + ) + yextremes = ( + max(0, iy - hshape[1] // 2), + min(iy + hshape[1] // 2 + 1, xdims[1]), + ) + zextremes = ( + max(0, iz - hshape[2] // 2), + min(iz + hshape[2] // 2 + 1, xdims[2]), + ) + + # find extremes of h (in case h is going out of model) + hxextremes = ( + max(0, -ix + hshape[0] // 2), + min(hshape[0], hshape[0] // 2 + (xdims[0] - ix)), + ) + hyextremes = ( + max(0, -iy + hshape[1] // 2), + min(hshape[1], hshape[1] // 2 + (xdims[1] - iy)), + ) + hzextremes = ( + max(0, -iz + hshape[2] // 2), + min(hshape[2], hshape[2] // 2 + (xdims[2] - iz)), + ) + + if not rmatvec: + y[ + xextremes[0] : xextremes[1], + yextremes[0] : yextremes[1], + zextremes[0] : zextremes[1], + ] += ( + x[ix, iy, iz] + * h[ + hxextremes[0] : hxextremes[1], + hyextremes[0] : hyextremes[1], + hzextremes[0] : hzextremes[1], + ] + ) + else: + y[ix, iy, iz] = np.sum( + h[ + hxextremes[0] : hxextremes[1], + hyextremes[0] : hyextremes[1], + hzextremes[0] : hzextremes[1], + ] + * x[ + xextremes[0] : xextremes[1], + yextremes[0] : yextremes[1], + zextremes[0] : zextremes[1], + ] + ) + return y + + @reshaped + def _matvec(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + y = ncp.zeros(self.dims, dtype=self.dtype) + y = self._mvrmv( + x, + y, + self.hs, + self.hshape, + self.dims, + self.ohx, + self.ohy, + self.ohz, + self.dhx, + self.dhy, + self.dhz, + self.nhx, + self.nhy, + self.nhz, + rmatvec=False, + **self.kwargs_cuda + ) + return y + + @reshaped + def _rmatvec(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + y = ncp.zeros(self.dims, dtype=self.dtype) + y = self._mvrmv( + x, + y, + self.hs, + self.hshape, + self.dims, + self.ohx, + self.ohy, + self.ohz, + self.dhx, + self.dhy, + self.dhz, + self.nhx, + self.nhy, + self.nhz, + rmatvec=True, + **self.kwargs_cuda + ) + return y diff --git a/pylops/utils/decorators.py b/pylops/utils/decorators.py index f8767b28..887d2d53 100644 --- a/pylops/utils/decorators.py +++ b/pylops/utils/decorators.py @@ -54,10 +54,16 @@ def add_ndarray_support_to_solver(func: Callable) -> Callable: @wraps(func) def wrapper(A, b, *args, **kwargs): # SciPy-type signature + x0flat = False if "x0" in kwargs and kwargs["x0"] is not None: + if kwargs["x0"].ndim == 1: + x0flat = True kwargs["x0"] = kwargs["x0"].ravel() with disabled_ndarray_multiplication(): res = list(func(A, b.ravel(), *args, **kwargs)) + # reshape if x0 was provided unflattened + # (unless the operator has forceflat=True) + if not x0flat and not getattr(A, "forceflat", None): res[0] = res[0].reshape(getattr(A, "dims", (A.shape[1],))) return tuple(res) diff --git a/pylops/utils/metrics.py b/pylops/utils/metrics.py index fa71967b..e1e7a55c 100644 --- a/pylops/utils/metrics.py +++ b/pylops/utils/metrics.py @@ -75,10 +75,10 @@ def snr(xref, xcmp): return snr -def psnr(xref, xcmp, xmax=None): +def psnr(xref, xcmp, xmax=None, xmin=0.0): """Peak Signal to Noise Ratio (PSNR) - Compute Peak Signal to Noise Ratio between two vectors. + Compute Peak Signal to Noise Ratio between two vectors Parameters ---------- @@ -89,6 +89,10 @@ def psnr(xref, xcmp, xmax=None): xmax : :obj:`float`, optional Maximum value to use. If ``None``, the actual maximum of the reference vector is used + xmin : :obj:`float`, optional + Minimum value to use. If ``None``, the actual minimum of + the reference vector is used (``0`` is default for + backward compatibility) Returns ------- @@ -98,5 +102,8 @@ def psnr(xref, xcmp, xmax=None): """ if xmax is None: xmax = xref.max() - psrn = 10.0 * np.log10(xmax**2 / mse(xref, xcmp)) - return psrn + if xmin is None: + xmin = xref.min() + xrange = xmax - xmin + psnr = 10.0 * np.log10(xrange**2 / mse(xref, xcmp)) + return psnr diff --git a/pylops/utils/signalprocessing.py b/pylops/utils/signalprocessing.py index 17917021..4aa277fa 100644 --- a/pylops/utils/signalprocessing.py +++ b/pylops/utils/signalprocessing.py @@ -5,6 +5,7 @@ "dip_estimate", ] +import warnings from typing import Tuple import numpy as np @@ -15,39 +16,49 @@ from pylops.utils.typing import NDArray -def convmtx(h: npt.ArrayLike, n: int) -> NDArray: +def convmtx(h: npt.ArrayLike, n: int, offset: int = 0) -> NDArray: r"""Convolution matrix - Equivalent of `MATLAB's convmtx function - `_ . Makes a dense convolution matrix :math:`\mathbf{C}` such that the dot product ``np.dot(C, x)`` is the convolution of - the filter :math:`h` and the input signal :math:`x`. + the filter :math:`h` centered on `offset` and the input signal :math:`x`. + + Equivalent of `MATLAB's convmtx function + `_ for: + - ``mode='full'`` when used with ``offset=0``. + - ``mode='same'`` when used with ``offset=len(h)//2`` (after truncating the rows as ``C[:n]``) Parameters ---------- h : :obj:`np.ndarray` Convolution filter (1D array) n : :obj:`int` - Number of columns (if :math:`\text{len}(h) < n`) or rows - (if :math:`\text{len}(h) \geq n`) of convolution matrix + Number of columns of convolution matrix + offset : :obj:`int` + Index of the center of the filter Returns ------- C : :obj:`np.ndarray` Convolution matrix of size :math:`\text{len}(h)+n-1 \times n` - (if :math:`\text{len}(h) < n`) or :math:`n \times \text{len}(h)+n-1` - (if :math:`\text{len}(h) \geq n`) """ + warnings.warn( + "A new implementation of convmtx is provided in v2.2.0 to match " + "MATLAB's convmtx method as stated in the docstring. The implementation " + "of convmtx provided prior to v2.2.0 was instead not consistent " + "with the documentation. Users are highly encouraged " + "to modify their codes accordingly.", + FutureWarning, + ) + ncp = get_array_module(h) - if len(h) < n: - col_1 = ncp.r_[h[0], ncp.zeros(n - 1, dtype=h.dtype)] - row_1 = ncp.r_[h, ncp.zeros(n - 1, dtype=h.dtype)] - else: - row_1 = ncp.r_[h[0], ncp.zeros(n - 1, dtype=h.dtype)] - col_1 = ncp.r_[h, ncp.zeros(n - 1, dtype=h.dtype)] + nh = len(h) + col_1 = ncp.r_[h, ncp.zeros(n + nh - 2, dtype=h.dtype)] + row_1 = ncp.r_[h[0], ncp.zeros(n - 1, dtype=h.dtype)] C = get_toeplitz(h)(col_1, row_1) + # apply offset + C = C[offset : offset + nh + n - 1] return C diff --git a/pylops/waveeqprocessing/blending.py b/pylops/waveeqprocessing/blending.py index 70242cb6..ca18ccb3 100644 --- a/pylops/waveeqprocessing/blending.py +++ b/pylops/waveeqprocessing/blending.py @@ -33,8 +33,9 @@ class BlendingContinuous(LinearOperator): Time sampling in seconds times : :obj:`np.ndarray` Absolute ignition times for each source - nproc : :obj:`int`, optional - Number of processors used when applying operator + shiftall : :obj:`bool`, optional + Shift all shots together (``True``) or one at the time (``False``). Defaults to ``shiftall=False`` (original + implementation), however ``shiftall=True`` should be preferred when ``nr`` is small. dtype : :obj:`str`, optional Operator dtype name : :obj:`str`, optional @@ -64,6 +65,7 @@ def __init__( ns: int, dt: float, times: NDArray, + shiftall: bool = False, dtype: DTypeLike = "float64", name: str = "B", ) -> None: @@ -73,40 +75,85 @@ def __init__( self.ns = ns self.dt = dt self.times = times + self.shiftall = shiftall self.nttot = int(np.max(self.times) / self.dt + self.nt + 1) - self.PadOp = Pad((self.nr, self.nt), ((0, 0), (0, 1)), dtype=self.dtype) - # Define shift operators - self.shifts = [] - self.ShiftOps = [] - for i in range(self.ns): - shift = self.times[i] - # This is the part that fits on the grid - shift_int = int(shift // self.dt) - self.shifts.append(shift_int) - # This is the fractional part - diff = (shift / self.dt - shift_int) * self.dt - if diff == 0: - self.ShiftOps.append(None) - else: - self.ShiftOps.append( - Shift( - (self.nr, self.nt + 1), - diff, - axis=1, - sampling=self.dt, - real=False, - dtype=self.dtype, + if not self.shiftall: + # original implementation, where each source is shifted indipendently + self.PadOp = Pad((self.nr, self.nt), ((0, 0), (0, 1)), dtype=self.dtype) + # Define shift operators + self.shifts = [] + self.ShiftOps = [] + for i in range(self.ns): + shift = self.times[i] + # This is the part that fits on the grid + shift_int = int(shift // self.dt) + self.shifts.append(shift_int) + # This is the fractional part + diff = (shift / self.dt - shift_int) * self.dt + if diff == 0: + self.ShiftOps.append(None) + else: + self.ShiftOps.append( + Shift( + (self.nr, self.nt + 1), + diff, + axis=1, + sampling=self.dt, + real=True, + dtype=self.dtype, + ) ) - ) + else: + # alternative implementation, where all sources are shifted at the same time + self.PadOp = Pad( + (self.ns, self.nr, self.nt), ((0, 0), (0, 0), (0, 1)), dtype=self.dtype + ) + # Define shift operator + self.shifts = (times // self.dt).astype(np.int32) + diff = (times / self.dt - self.shifts) * self.dt + diff = np.repeat(diff[:, np.newaxis], self.nr, axis=1) + self.ShiftOp = Shift( + (self.ns, self.nr, self.nt + 1), + diff, + axis=-1, + sampling=self.dt, + real=True, + dtype=self.dtype, + ) + self.diff = diff + super().__init__( dtype=np.dtype(dtype), dims=(self.ns, self.nr, self.nt), dimsd=(self.nr, self.nttot), name=name, ) + self._register_multiplications() + + @reshaped + def _matvec_smallrecs(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + blended_data = ncp.zeros((self.nr, self.nttot), dtype=self.dtype) + shifted_data = self.ShiftOp._matvec(self.PadOp._matvec(x.ravel())).reshape( + self.ns, self.nr, self.nt + 1 + ) + for i, shift_int in enumerate(self.shifts): + blended_data[:, shift_int : shift_int + self.nt + 1] += shifted_data[i] + return blended_data + + @reshaped + def _rmatvec_smallrecs(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + shifted_data = ncp.zeros((self.ns, self.nr, self.nt + 1), dtype=self.dtype) + for i, shift_int in enumerate(self.shifts): + shifted_data[i, :, :] = x[:, shift_int : shift_int + self.nt + 1] + deblended_data = self.PadOp._rmatvec( + self.ShiftOp._rmatvec(shifted_data.ravel()) + ).reshape(self.dims) + return deblended_data @reshaped - def _matvec(self, x: NDArray) -> NDArray: + def _matvec_largerecs(self, x: NDArray) -> NDArray: ncp = get_array_module(x) blended_data = ncp.zeros((self.nr, self.nttot), dtype=self.dtype) for i, shift_int in enumerate(self.shifts): @@ -118,7 +165,7 @@ def _matvec(self, x: NDArray) -> NDArray: return blended_data @reshaped - def _rmatvec(self, x: NDArray) -> NDArray: + def _rmatvec_largerecs(self, x: NDArray) -> NDArray: ncp = get_array_module(x) deblended_data = ncp.zeros((self.ns, self.nr, self.nt), dtype=self.dtype) for i, shift_int in enumerate(self.shifts): @@ -133,6 +180,14 @@ def _rmatvec(self, x: NDArray) -> NDArray: deblended_data[i, :, :] = shifted_data return deblended_data + def _register_multiplications(self) -> None: + if self.shiftall: + self._matvec = self._matvec_smallrecs + self._rmatvec = self._rmatvec_smallrecs + else: + self._matvec = self._matvec_largerecs + self._rmatvec = self._rmatvec_largerecs + def BlendingGroup( nt: int, diff --git a/pylops/waveeqprocessing/kirchhoff.py b/pylops/waveeqprocessing/kirchhoff.py index 595664df..06c29251 100644 --- a/pylops/waveeqprocessing/kirchhoff.py +++ b/pylops/waveeqprocessing/kirchhoff.py @@ -35,10 +35,12 @@ class Kirchhoff(LinearOperator): - r"""Kirchhoff Demigration operator. + r"""Kirchhoff demigration operator. Kirchhoff-based demigration/migration operator. Uses a high-frequency - approximation of Green's function propagators based on ``trav``. + approximation of the Green's function propagators based on traveltimes + and amplitudes that are either computed internally by solving the Eikonal equation, + or passed directly by the user (which can use any other propagation engine of choice). Parameters ---------- @@ -73,23 +75,21 @@ class Kirchhoff(LinearOperator): dynamic : :obj:`bool`, optional .. versionadded:: 2.0.0 - Include dynamic weights in computations (``True``) or not (``False``). - Note that when ``mode=byot``, the user is required to provide such weights - in ``amp``. + Apply dynamic weights in computations (``True``) or not (``False``). This includes both the amplitude + terms of the Green's function and the reflectivity-related scaling term (see equations below). trav : :obj:`numpy.ndarray` or :obj:`tuple`, optional Traveltime table of size :math:`\lbrack (n_y) n_x n_z \times n_s n_r \rbrack` or pair of traveltime tables of size :math:`\lbrack (n_y) n_x n_z \times n_s \rbrack` and :math:`\lbrack (n_y) n_x n_z \times n_r \rbrack` (to be provided if ``mode='byot'``). Note that the latter approach is recommended as less memory demanding - than the former. + than the former. Moreover, only ``mode='dynamic'`` is only possible when traveltimes are provided in + the latter form. amp : :obj:`numpy.ndarray`, optional .. versionadded:: 2.0.0 - Amplitude table of size - :math:`\lbrack (n_y) n_x n_z \times n_s n_r \rbrack` or pair of amplitude tables - of size :math:`\lbrack (n_y) n_x n_z \times n_s \rbrack` and :math:`\lbrack (n_y) n_x n_z \times n_r \rbrack` - (to be provided if ``mode='byot'``). Note that the latter approach is recommended as less memory demanding - than the former. + Pair of amplitude tables of size :math:`\lbrack (n_y) n_x n_z \times n_s \rbrack` and + :math:`\lbrack (n_y) n_x n_z \times n_r \rbrack` (to be provided if ``mode='byot'``). Note that this parameter + is only required when ``mode='dynamic'`` is chosen. aperture : :obj:`float` or :obj:`tuple`, optional .. versionadded:: 2.0.0 @@ -104,18 +104,9 @@ class Kirchhoff(LinearOperator): Maximum allowed angle (either source or receiver side) in degrees. If ``None``, angle aperture limitations are not introduced. See ``aperture`` for implementation details regarding scalar and tuple cases. - - anglerefl : :obj:`np.ndarray`, optional - .. versionadded:: 2.0.0 - - Angle between the normal of the reflectors and the vertical axis in degrees snell : :obj:`float` or :obj:`tuple`, optional - .. versionadded:: 2.0.0 - - Threshold on Snell's law evaluation. If larger, the source-receiver-image - point is discarded. If ``None``, no check on the validity of the Snell's - law is performed. See ``aperture`` for implementation details regarding - scalar and tuple cases. + Deprecated, will be removed in v3.0.0. Simply kept for back-compatibility with previous implementation, + but effectively not affecting the behaviour of the operator. engine : :obj:`str`, optional Engine used for computations (``numpy`` or ``numba``). dtype : :obj:`str`, optional @@ -141,29 +132,34 @@ class Kirchhoff(LinearOperator): Notes ----- The Kirchhoff demigration operator synthesizes seismic data given a - propagation velocity model :math:`v` and a reflectivity model :math:`m`. - In forward mode [1]_, [2]_: + propagation velocity model :math:`v(\mathbf{x})` and a + reflectivity model :math:`m(\mathbf{x})`. In forward mode [1]_, [2]_, [3]_: .. math:: d(\mathbf{x_r}, \mathbf{x_s}, t) = - \widetilde{w}(t) * \int_V G(\mathbf{x_r}, \mathbf{x}, t) - m(\mathbf{x}) G(\mathbf{x}, \mathbf{x_s}, t)\,\mathrm{d}\mathbf{x} - - where :math:`m(\mathbf{x})` represents the reflectivity - at every location in the subsurface, :math:`G(\mathbf{x}, \mathbf{x_s}, t)` - and :math:`G(\mathbf{x_r}, \mathbf{x}, t)` are the Green's functions - from source-to-subsurface-to-receiver and finally :math:`\widetilde{w}(t)` is - a filtered version of the wavelet :math:`w(t)` [3]_ (or the wavelet itself when - ``wavfilter=False``). In our implementation, the following high-frequency - approximation of the Green's functions is adopted: + \widetilde{w}(t) * \int_V \frac{2 \cos\theta} {v(\mathbf{x})} + G(\mathbf{x_r}, \mathbf{x}, t) G(\mathbf{x}, \mathbf{x_s}, t) + m(\mathbf{x}) \,\mathrm{d}\mathbf{x} + + where :math:`G(\mathbf{x}, \mathbf{x_s}, t)` and :math:`G(\mathbf{x_r}, \mathbf{x}, t)` + are the Green's functions from source-to-subsurface-to-receiver and finally + :math:`\widetilde{w}(t)` is either a filtered version of the wavelet :math:`w(t)` + as explained below (``wavfilter=True``) or the wavelet itself when (``wavfilter=False``). + Moreover, an angle scaling is included in the modelling operator, + where the reflection angle :math:`\theta=(\theta_s-\theta_r)/2` is half of the opening angle, + with :math:`\theta_s` and :math:`\theta_r` representing the angles between the source-side + and receiver-side rays and the vertical at the image point, respectively. + + In our implementation, the following high-frequency approximation of the + Green's functions is adopted: .. math:: G(\mathbf{x_r}, \mathbf{x}, \omega) = a(\mathbf{x_r}, \mathbf{x}) e^{j \omega t(\mathbf{x_r}, \mathbf{x})} where :math:`a(\mathbf{x_r}, \mathbf{x})` is the amplitude and - :math:`t(\mathbf{x_r}, \mathbf{x})` is the traveltime. When ``dynamic=False`` the - amplitude is disregarded leading to a kinematic-only Kirchhoff operator. + :math:`t(\mathbf{x_r}, \mathbf{x})` is the traveltime. When ``dynamic=False``, the + amplitude correction terms are disregarded leading to a kinematic-only Kirchhoff operator. .. math:: d(\mathbf{x_r}, \mathbf{x_s}, t) = @@ -171,36 +167,32 @@ class Kirchhoff(LinearOperator): t(\mathbf{x}, \mathbf{x_s}))} m(\mathbf{x}) \,\mathrm{d}\mathbf{x} On the other hand, when ``dynamic=True``, the amplitude scaling is defined as - :math:`a(\mathbf{x}, \mathbf{y})=\frac{1}{\|\mathbf{x} - \mathbf{y}\|}`, - that is, the reciprocal of the distance between the two points, - approximating the geometrical spreading of the wavefront. - Moreover an angle scaling is included in the modelling operator - added as follows: - .. math:: - d(\mathbf{x_r}, \mathbf{x_s}, t) = - \tilde{w}(t) * \int_V a(\mathbf{x}, \mathbf{x_s}) a(\mathbf{x}, \mathbf{x_r}) - \frac{|cos \theta_s + cos \theta_r|} {v(\mathbf{x})} e^{j \omega (t(\mathbf{x_r}, \mathbf{x}) + - t(\mathbf{x}, \mathbf{x_s}))} m(\mathbf{x}) \,\mathrm{d}\mathbf{x} + * ``2D``: :math:`a(\mathbf{x}, \mathbf{y})=\frac{1}{\sqrt{\text{dist}(\mathbf{x}, \mathbf{y})}}` + * ``3D``: :math:`a(\mathbf{x}, \mathbf{y})=\frac{1}{\text{dist}(\mathbf{x}, \mathbf{y})}` + + approximating the geometrical spreading of the wavefront. For ``mode=analytic``, + :math:`\text{dist}(\mathbf{x}, \mathbf{y})=\|\mathbf{x} - \mathbf{y}\|`, whilst for + ``mode=eikonal``, this is computed internally by the Eikonal solver. - where :math:`\theta_s` and :math:`\theta_r` are the angles between the source-side - and receiver-side rays and the normal to the reflector at the image point (or - the vertical axis at the image point when ``reflslope=None``), respectively. + The wavelet filtering is applied as follows [4]_: + + * ``2D``: :math:`\tilde{W}(f)=\sqrt{j\omega} \cdot W(f)` + * ``3D``: :math:`\tilde{W}(f)=-j\omega \cdot W(f)` Depending on the choice of ``mode`` the traveltime and amplitude of the Green's function will be also computed differently: - * ``mode=analytic`` or ``mode=eikonal``: traveltimes, geometrical spreading, and angles - are computed for every source-image point-receiver triplets and the - Green's functions are implemented from traveltime look-up tables, placing - scaled reflectivity values at corresponding source-to-receiver time in the - data. - * ``byot``: bring your own tables. Traveltime table are provided + * ``mode=analytic`` or ``mode=eikonal``: traveltimes, amplitudes, and angles + are computed for every source-image point-receiver triplets upfront and the + Green's functions are implemented from traveltime and amplitude look-up tables, + placing scaled reflectivity values at corresponding source-to-receiver time + in the data. + * ``mode=byot``: bring your own tables. Traveltime table are provided directly by user using ``trav`` input parameter. Similarly, in this case one - can provide their own amplitude scaling ``amp`` (which should include the angle - scaling too). + can also provide their own amplitude scaling ``amp``. - Three aperture limitations have been also implemented as defined by: + Two aperture limitations have been also implemented as defined by: * ``aperture``: the maximum allowed aperture is expressed as the ratio of offset over depth. This aperture limitation avoid including grazing angles @@ -208,28 +200,27 @@ class Kirchhoff(LinearOperator): edges of the aperture; * ``angleaperture``: the maximum allowed angle aperture is expressed as the difference between the incident or emerging angle at every image point and - the vertical axis (or the normal to the reflector if ``anglerefl`` is provided. - This aperture limitation also avoid including grazing angles whose contributions - can introduce aliasing effects. Note that for a homogenous medium and slowly varying - heterogenous medium the offset and angle aperture limits may work in the same way; - * ``snell``: the maximum allowed snell's angle is expressed as the absolute value of - the sum between incident and emerging angles defined as in the ``angleaperture`` case. - This aperture limitation is introduced to turn a scattering-based Kirchhoff engine into - a reflection-based Kirchhoff engine where each image point is not considered as scatter - but as a local horizontal (or straight) reflector. + the vertical axis. This aperture limitation also avoid including grazing angles + whose contributions can introduce aliasing effects. Note that for a homogenous + medium and slowly varying heterogeneous medium the offset and angle aperture limits + may work in the same way. Finally, the adjoint of the demigration operator is a *migration* operator which projects data in the model domain creating an image of the subsurface reflectivity. - .. [1] Bleistein, N., Cohen, J.K., and Stockwell, J.W.. + .. [1] Bleistein, N., Cohen, J.K., and Stockwell, J.W. "Mathematics of Multidimensional Seismic Imaging, Migration and Inversion", 2000. .. [2] Santos, L.T., Schleicher, J., Tygel, M., and Hubral, P. "Seismic modeling by demigration", Geophysics, 65(4), pp. 1281-1289, 2000. - .. [3] Safron, L. "Multicomponent least-squares Kirchhoff depth migration", + .. [3] Yang, K., and Zhang, J. "Comparison between Born and Kirchhoff operators for + least-squares reverse time migration and the constraint of the propagation of the + background wavefield", Geophysics, 84(5), pp. R725-R739, 2019. + + .. [4] Safron, L. "Multicomponent least-squares Kirchhoff depth migration", MSc Thesis, 2018. """ @@ -252,7 +243,6 @@ def __init__( amp: Optional[NDArray] = None, aperture: Optional[Tuple[float, float]] = None, angleaperture: Union[float, Tuple[float, float]] = 90.0, - anglerefl: Optional[NDArray] = None, snell: Optional[Tuple[float, float]] = None, engine: str = "numpy", dtype: DTypeLike = "float64", @@ -288,17 +278,20 @@ def __init__( self.dt = t[1] - t[0] self.nt = len(t) - # store ix-iy locations of sources and receivers - dx = x[1] - x[0] - if self.ndims == 2: - self.six = np.tile((srcs[0] - x[0]) // dx, (nr, 1)).T.astype(int).ravel() - self.rix = np.tile((recs[0] - x[0]) // dx, (ns, 1)).astype(int).ravel() - elif self.ndims == 3: - # TODO: 3D normalized distances - pass - - # compute traveltime + # store ix-iy locations of sources and receivers for aperture filter self.dynamic = dynamic + if self.dynamic: + dx = x[1] - x[0] + if self.ndims == 2: + self.six = ( + np.tile((srcs[0] - x[0]) // dx, (nr, 1)).T.astype(int).ravel() + ) + self.rix = np.tile((recs[0] - x[0]) // dx, (ns, 1)).astype(int).ravel() + elif self.ndims == 3: + # TODO: 3D normalized distances + raise NotImplementedError("dynamic=True currently not available in 3D") + + # compute traveltime and distances self.travsrcrec = True # use separate tables for src and rec traveltimes if mode in ["analytic", "eikonal", "byot"]: if mode in ["analytic", "eikonal"]: @@ -306,63 +299,72 @@ def __init__( ( self.trav_srcs, self.trav_recs, - self.dist_srcs, - self.dist_recs, + dist_srcs, + dist_recs, trav_srcs_grad, trav_recs_grad, ) = Kirchhoff._traveltime_table(z, x, srcs, recs, vel, y=y, mode=mode) if self.dynamic: # need to add a scalar in the denominator of amplitude term to avoid - # division by 0, currently set to 1/100 of max distance to avoid having - # to large scaling around the source. This number may change in future + # division by 0, currently set to 1e-2 of max distance to avoid having + # too large scaling around the source. This number may change in future # or left to the user to define epsdist = 1e-2 - self.maxdist = epsdist * ( - np.max(self.dist_srcs) + np.max(self.dist_recs) - ) - # compute angles + self.maxdist = epsdist * (np.max(dist_srcs) + np.max(dist_recs)) if self.ndims == 2: - # 2d with vertical - if anglerefl is None: - self.angle_srcs = np.arctan2( - trav_srcs_grad[0], trav_srcs_grad[1] - ).reshape(np.prod(dims), ns) - self.angle_recs = np.arctan2( - trav_recs_grad[0], trav_recs_grad[1] - ).reshape(np.prod(dims), nr) - self.cosangle_srcs = np.cos(self.angle_srcs) - self.cosangle_recs = np.cos(self.angle_recs) - else: - # TODO: 2D with normal - raise NotImplementedError( - "angle scaling with anglerefl currently not available" - ) + self.amp_srcs, self.amp_recs = 1.0 / np.sqrt( + dist_srcs + self.maxdist + ), 1.0 / np.sqrt(dist_recs + self.maxdist) else: - # TODO: 3D - raise NotImplementedError( - "dynamic=True currently not available in 3D" - ) + self.amp_srcs, self.amp_recs = 1.0 / ( + dist_srcs + self.maxdist + ), 1.0 / (dist_recs + self.maxdist) else: if isinstance(trav, tuple): self.trav_srcs, self.trav_recs = trav else: self.travsrcrec = False self.trav = trav + + if self.dynamic and not self.travsrcrec: + raise NotImplementedError( + "separate traveltime tables must be provided " + "when selecting mode=dynamic" + ) if self.dynamic: if isinstance(amp, tuple): self.amp_srcs, self.amp_recs = amp else: - self.amp = amp - # in byot mode, angleaperture and snell checks are not performed - self.angle_srcs = np.ones( - (self.ny * self.nx * self.nz, ns), dtype=dtype - ) - self.angle_recs = np.ones( - (self.ny * self.nx * self.nz, nr), dtype=dtype - ) + raise NotImplementedError( + "separate amplitude tables must be provided " + ) + + if self.travsrcrec: + # compute traveltime gradients at image points + trav_srcs_grad = np.gradient( + self.trav_srcs.reshape(*dims, ns), + axis=np.arange(self.ndims), + ) + trav_recs_grad = np.gradient( + self.trav_recs.reshape(*dims, nr), + axis=np.arange(self.ndims), + ) else: raise NotImplementedError("method must be analytic, eikonal or byot") + # compute angles with vertical + if self.dynamic: + if self.ndims == 2: + self.angle_srcs = np.arctan2( + trav_srcs_grad[0], trav_srcs_grad[1] + ).reshape(np.prod(dims), ns) + self.angle_recs = np.arctan2( + trav_recs_grad[0], trav_recs_grad[1] + ).reshape(np.prod(dims), nr) + else: + # TODO: 3D + raise NotImplementedError("dynamic=True currently not available in 3D") + # pre-compute traveltime indices if total traveltime is used if not self.travsrcrec: self.itrav = (self.trav / self.dt).astype("int32") @@ -383,31 +385,25 @@ def __init__( self.aperturetap = taper(41, 20, "hanning")[20:] # define aperture + # if aperture=None, we want to ensure the check is always matched (no aperture limits...) if aperture is not None: warnings.warn( "Aperture is currently defined as ratio of offset over depth, " - "and may be not ideal for highly heterogenous media" + "and may be not ideal for highly heterogeneous media" ) self.aperture = ( - (2 * self.nx / self.nz,) - if aperture is None - else _value_or_sized_to_array(aperture) + (self.nx + 1,) if aperture is None else _value_or_sized_to_array(aperture) ) if len(self.aperture) == 1: self.aperture = np.array([0.8 * self.aperture[0], self.aperture[0]]) - # define angle aperture and snell law + # define angle aperture angleaperture = [0.0, 1000.0] if angleaperture is None else angleaperture self.angleaperture = np.deg2rad(_value_or_sized_to_array(angleaperture)) if len(self.angleaperture) == 1: self.angleaperture = np.array( [0.8 * self.angleaperture[0], self.angleaperture[0]] ) - self.snell = ( - (np.pi,) if snell is None else np.deg2rad(_value_or_sized_to_array(snell)) - ) - if len(self.snell) == 1: - self.snell = np.array([0.8 * self.snell[0], self.snell[0]]) # dimensions self.ns, self.nr = ns, nr @@ -578,7 +574,7 @@ def _traveltime_table( dist_recs = trav_recs * vel elif mode == "eikonal": - if skfmm is not None: + if skfmm_message is None: dist_srcs = np.zeros((ny * nx * nz, ns)) dist_recs = np.zeros((ny * nx * nz, nr)) trav_srcs = np.zeros((ny * nx * nz, ns)) @@ -635,12 +631,13 @@ def _wavelet_reshaping( dimrec: int, dimv: int, ) -> NDArray: - """Apply wavelet reshaping as from theory in [1]_""" + """Apply wavelet reshaping to account for omega scaling factor + originating from the wave equation""" f = np.fft.rfftfreq(len(wav), dt) W = np.fft.rfft(wav, n=len(wav)) if dimsrc == 2 and dimv == 2: # 2D - Wfilt = W * (2 * np.pi * f) + Wfilt = W * np.sqrt(1j * 2 * np.pi * f) elif (dimsrc == 2 or dimrec == 2) and dimv == 3: # 2.5D raise NotImplementedError("2.D wavelet currently not available") @@ -750,225 +747,6 @@ def _travsrcrec_kirch_rmatvec( ) return y - @staticmethod - def _amp_kirch_matvec( - x: NDArray, - y: NDArray, - nsnr: int, - nt: int, - ni: int, - itrav: NDArray, - travd: NDArray, - amp: NDArray, - aperturemin: float, - aperturemax: float, - aperturetap: NDArray, - nz: int, - six: NDArray, - rix: NDArray, - angleaperturemin: float, - angleaperturemax: float, - angles_srcs: NDArray, - angles_recs: NDArray, - snellmin: float, - snellmax: float, - ) -> NDArray: - nr = angles_recs.shape[-1] - daperture = aperturemax - aperturemin - dangleaperture = angleaperturemax - angleaperturemin - dsnell = snellmax - snellmin - for isrcrec in prange(nsnr): - # extract traveltime, amplitude, src/rec coordinates at given src/pair - itravisrcrec = itrav[:, isrcrec] - travdisrcrec = travd[:, isrcrec] - ampisrcrec = amp[:, isrcrec] - sixisrcrec = six[isrcrec] - rixisrcrec = rix[isrcrec] - # extract source and receiver angles - angles_src = angles_srcs[:, isrcrec // nr] - angles_rec = angles_recs[:, isrcrec % nr] - for ii in range(ni): - # extract traveltime, amplitude at given image point - itravii = itravisrcrec[ii] - travdii = travdisrcrec[ii] - damp = ampisrcrec[ii] - # extract source and receiver angle - angle_src = angles_src[ii] - angle_rec = angles_rec[ii] - abs_angle_src = abs(angle_src) - abs_angle_rec = abs(angle_rec) - abs_angle_src_rec = abs(angle_src + angle_rec) - aptap = 1.0 - # angle apertures checks - if ( - abs_angle_src < angleaperturemax - and abs_angle_rec < angleaperturemax - and abs_angle_src_rec < snellmax - ): - if abs_angle_src >= angleaperturemin: - # extract source angle aperture taper value - aptap = ( - aptap - * aperturetap[ - int( - 20 - * (abs_angle_src - angleaperturemin) - // dangleaperture - ) - ] - ) - if abs_angle_rec >= angleaperturemin: - # extract receiver angle aperture taper value - aptap = ( - aptap - * aperturetap[ - int( - 20 - * (abs_angle_rec - angleaperturemin) - // dangleaperture - ) - ] - ) - if abs_angle_src_rec >= snellmin: - # extract snell taper value - aptap = ( - aptap - * aperturetap[ - int(20 * (abs_angle_src_rec - snellmin) // dsnell) - ] - ) - - # identify x-index of image point - iz = ii % nz - # aperture check - aperture = abs(sixisrcrec - rixisrcrec) / iz - if aperture < aperturemax: - if aperture >= aperturemin: - # extract aperture taper value - aptap = ( - aptap - * aperturetap[ - int(20 * ((aperture - aperturemin) // daperture)) - ] - ) - # time limit check - if 0 <= itravii < nt - 1: - # assign values - y[isrcrec, itravii] += x[ii] * (1 - travdii) * damp * aptap - y[isrcrec, itravii + 1] += x[ii] * travdii * damp * aptap - return y - - @staticmethod - def _amp_kirch_rmatvec( - x: NDArray, - y: NDArray, - nsnr: int, - nt: int, - ni: int, - itrav: NDArray, - travd: NDArray, - amp: NDArray, - aperturemin: float, - aperturemax: float, - aperturetap: NDArray, - nz: int, - six: NDArray, - rix: NDArray, - angleaperturemin: float, - angleaperturemax: float, - angles_srcs: NDArray, - angles_recs: NDArray, - snellmin: float, - snellmax: float, - ) -> NDArray: - nr = angles_recs.shape[-1] - daperture = aperturemax - aperturemin - dangleaperture = angleaperturemax - angleaperturemin - dsnell = snellmax - snellmin - for ii in prange(ni): - itravii = itrav[ii] - travdii = travd[ii] - ampii = amp[ii] - # extract source and receiver angles - angle_srcs = angles_srcs[ii] - angle_recs = angles_recs[ii] - # identify x-index of image point - iz = ii % nz - for isrcrec in range(nsnr): - itravisrcrecii = itravii[isrcrec] - travdisrcrecii = travdii[isrcrec] - sixisrcrec = six[isrcrec] - rixisrcrec = rix[isrcrec] - # extract source and receiver angle - angle_src = angle_srcs[isrcrec // nr] - angle_rec = angle_recs[isrcrec % nr] - abs_angle_src = abs(angle_src) - abs_angle_rec = abs(angle_rec) - abs_angle_src_rec = abs(angle_src + angle_rec) - aptap = 1.0 - # angle apertures checks - if ( - abs_angle_src < angleaperturemax - and abs_angle_rec < angleaperturemax - and abs_angle_src_rec < snellmax - ): - if abs_angle_src >= angleaperturemin: - # extract source angle aperture taper value - aptap = ( - aptap - * aperturetap[ - int( - 20 - * (abs_angle_src - angleaperturemin) - // dangleaperture - ) - ] - ) - if abs_angle_rec >= angleaperturemin: - # extract receiver angle aperture taper value - aptap = ( - aptap - * aperturetap[ - int( - 20 - * (abs_angle_rec - angleaperturemin) - // dangleaperture - ) - ] - ) - if abs_angle_src_rec >= snellmin: - # extract snell taper value - aptap = ( - aptap - * aperturetap[ - int(20 * (abs_angle_src_rec - snellmin) // dsnell) - ] - ) - - # aperture check - aperture = abs(sixisrcrec - rixisrcrec) / iz - if aperture < aperturemax: - if aperture >= aperturemin: - # extract aperture taper value - aptap = ( - aptap - * aperturetap[ - int(20 * ((aperture - aperturemin) // daperture)) - ] - ) - # time limit check - if 0 <= itravisrcrecii < nt - 1: - # assign values - y[ii] += ( - ( - x[isrcrec, itravisrcrecii] * (1 - travdisrcrecii) - + x[isrcrec, itravisrcrecii + 1] * travdisrcrecii - ) - * ampii[isrcrec] - * aptap - ) - return y - @staticmethod def _ampsrcrec_kirch_matvec( x: NDArray, @@ -981,8 +759,8 @@ def _ampsrcrec_kirch_matvec( vel: NDArray, trav_srcs: NDArray, trav_recs: NDArray, - dist_srcs: NDArray, - dist_recs: NDArray, + amp_srcs: NDArray, + amp_recs: NDArray, aperturemin: float, aperturemax: float, aperturetap: NDArray, @@ -993,44 +771,39 @@ def _ampsrcrec_kirch_matvec( angleaperturemax: float, angles_srcs: NDArray, angles_recs: NDArray, - snellmin: float, - snellmax: float, - maxdist: float, ) -> NDArray: daperture = aperturemax - aperturemin dangleaperture = angleaperturemax - angleaperturemin - dsnell = snellmax - snellmin for isrc in prange(ns): travisrc = trav_srcs[:, isrc] - distisrc = dist_srcs[:, isrc] + ampisrc = amp_srcs[:, isrc] angleisrc = angles_srcs[:, isrc] for irec in range(nr): travirec = trav_recs[:, irec] trav = travisrc + travirec itrav = (trav / dt).astype("int32") travd = trav / dt - itrav - distirec = dist_recs[:, irec] + ampirec = amp_recs[:, irec] angleirec = angles_recs[:, irec] - dist = distisrc + distirec - amp = np.abs(np.cos(angleisrc) + np.cos(angleirec)) / (dist + maxdist) sixisrcrec = six[isrc * nr + irec] rixisrcrec = rix[isrc * nr + irec] + # compute cosine of half opening angle and total amplitude scaling + cosangle = np.cos((angleisrc - angleirec) / 2.0) + amp = 2.0 * cosangle * ampisrc * ampirec / vel for ii in range(ni): itravii = itrav[ii] travdii = travd[ii] - damp = amp[ii] / vel[ii] + damp = amp[ii] # extract source and receiver angle at given image point angle_src = angleisrc[ii] angle_rec = angleirec[ii] abs_angle_src = abs(angle_src) abs_angle_rec = abs(angle_rec) - abs_angle_src_rec = abs(angle_src + angle_rec) - aptap = 1.0 # angle apertures checks + aptap = 1.0 if ( abs_angle_src < angleaperturemax and abs_angle_rec < angleaperturemax - and abs_angle_src_rec < snellmax ): if abs_angle_src >= angleaperturemin: # extract source angle aperture taper value @@ -1056,19 +829,11 @@ def _ampsrcrec_kirch_matvec( ) ] ) - if abs_angle_src_rec >= snellmin: - # extract snell taper value - aptap = ( - aptap - * aperturetap[ - int(20 * (abs_angle_src_rec - snellmin) // dsnell) - ] - ) # identify x-index of image point iz = ii % nz # aperture check - aperture = abs(sixisrcrec - rixisrcrec) / iz + aperture = abs(sixisrcrec - rixisrcrec) / (iz + 1) if aperture < aperturemax: if aperture >= aperturemin: # extract aperture taper value @@ -1102,8 +867,8 @@ def _ampsrcrec_kirch_rmatvec( vel: NDArray, trav_srcs: NDArray, trav_recs: NDArray, - dist_srcs: NDArray, - dist_recs: NDArray, + amp_srcs: NDArray, + amp_recs: NDArray, aperturemin: float, aperturemax: float, aperturetap: NDArray, @@ -1114,18 +879,14 @@ def _ampsrcrec_kirch_rmatvec( angleaperturemax: float, angles_srcs: NDArray, angles_recs: NDArray, - snellmin: float, - snellmax: float, - maxdist: float, ) -> NDArray: daperture = aperturemax - aperturemin dangleaperture = angleaperturemax - angleaperturemin - dsnell = snellmax - snellmin for ii in prange(ni): trav_srcsii = trav_srcs[ii] trav_recsii = trav_recs[ii] - dist_srcsii = dist_srcs[ii] - dist_recsii = dist_recs[ii] + amp_srcsii = amp_srcs[ii] + amp_recsii = amp_recs[ii] velii = vel[ii] angle_srcsii = angles_srcs[ii] angle_recsii = angles_recs[ii] @@ -1133,31 +894,27 @@ def _ampsrcrec_kirch_rmatvec( iz = ii % nz for isrc in range(ns): trav_srcii = trav_srcsii[isrc] - dist_srcii = dist_srcsii[isrc] + amp_srcii = amp_srcsii[isrc] angle_src = angle_srcsii[isrc] for irec in range(nr): trav_recii = trav_recsii[irec] travii = trav_srcii + trav_recii itravii = int(travii / dt) travdii = travii / dt - itravii - dist_recii = dist_recsii[irec] + amp_recii = amp_recsii[irec] angle_rec = angle_recsii[irec] - dist = dist_srcii + dist_recii - ampii = np.abs(np.cos(angle_src) + np.cos(angle_rec)) / ( - dist + maxdist - ) - damp = ampii / velii sixisrcrec = six[isrc * nr + irec] rixisrcrec = rix[isrc * nr + irec] abs_angle_src = abs(angle_src) abs_angle_rec = abs(angle_rec) - abs_angle_src_rec = abs(angle_src + angle_rec) - aptap = 1.0 + # compute cosine of half opening angle and total amplitude scaling + cosangle = np.cos((angle_src - angle_rec) / 2.0) + damp = 2.0 * cosangle * amp_srcii * amp_recii / velii # angle apertures checks + aptap = 1.0 if ( abs_angle_src < angleaperturemax and abs_angle_rec < angleaperturemax - and abs_angle_src_rec < snellmax ): if abs_angle_src >= angleaperturemin: # extract source angle aperture taper value @@ -1183,17 +940,9 @@ def _ampsrcrec_kirch_rmatvec( ) ] ) - if abs_angle_src_rec >= snellmin: - # extract snell taper value - aptap = ( - aptap - * aperturetap[ - int(20 * (abs_angle_src_rec - snellmin) // dsnell) - ] - ) # aperture check - aperture = abs(sixisrcrec - rixisrcrec) / iz + aperture = abs(sixisrcrec - rixisrcrec) / (iz + 1) if aperture < aperturemax: if aperture >= aperturemin: # extract aperture taper value @@ -1221,7 +970,7 @@ def _ampsrcrec_kirch_rmatvec( def _register_multiplications(self, engine: str) -> None: if engine not in ["numpy", "numba"]: raise KeyError("engine must be numpy or numba") - if engine == "numba" and jit is not None: + if engine == "numba" and jit_message is None: # numba numba_opts = dict( nopython=True, nogil=True, parallel=parallel @@ -1229,9 +978,6 @@ def _register_multiplications(self, engine: str) -> None: if self.dynamic and self.travsrcrec: self._kirch_matvec = jit(**numba_opts)(self._ampsrcrec_kirch_matvec) self._kirch_rmatvec = jit(**numba_opts)(self._ampsrcrec_kirch_rmatvec) - elif self.dynamic and not self.travsrcrec: - self._kirch_matvec = jit(**numba_opts)(self._amp_kirch_matvec) - self._kirch_rmatvec = jit(**numba_opts)(self._amp_kirch_rmatvec) elif self.travsrcrec: self._kirch_matvec = jit(**numba_opts)(self._travsrcrec_kirch_matvec) self._kirch_rmatvec = jit(**numba_opts)(self._travsrcrec_kirch_rmatvec) @@ -1240,14 +986,11 @@ def _register_multiplications(self, engine: str) -> None: self._kirch_rmatvec = jit(**numba_opts)(self._trav_kirch_rmatvec) else: - if engine == "numba" and jit is None: + if engine == "numba" and jit_message is not None: logging.warning(jit_message) if self.dynamic and self.travsrcrec: self._kirch_matvec = self._ampsrcrec_kirch_matvec self._kirch_rmatvec = self._ampsrcrec_kirch_rmatvec - elif self.dynamic and not self.travsrcrec: - self._kirch_matvec = self._amp_kirch_matvec - self._kirch_rmatvec = self._amp_kirch_rmatvec elif self.travsrcrec: self._kirch_matvec = self._travsrcrec_kirch_matvec self._kirch_rmatvec = self._travsrcrec_kirch_rmatvec @@ -1270,32 +1013,8 @@ def _matvec(self, x: NDArray) -> NDArray: self.vel, self.trav_srcs, self.trav_recs, - self.dist_srcs, - self.dist_recs, - self.aperture[0], - self.aperture[1], - self.aperturetap, - self.nz, - self.six, - self.rix, - self.angleaperture[0], - self.angleaperture[1], - self.angle_srcs, - self.angle_recs, - self.snell[0], - self.snell[1], - self.maxdist, - ) - elif self.dynamic and not self.travsrcrec: - inputs = ( - x.ravel(), - y, - self.nsnr, - self.nt, - self.ni, - self.itrav, - self.travd, - self.amp, + self.amp_srcs, + self.amp_recs, self.aperture[0], self.aperture[1], self.aperturetap, @@ -1306,10 +1025,8 @@ def _matvec(self, x: NDArray) -> NDArray: self.angleaperture[1], self.angle_srcs, self.angle_recs, - self.snell[0], - self.snell[1], ) - elif not self.dynamic and self.travsrcrec: + elif self.travsrcrec: inputs = ( x.ravel(), y, @@ -1321,7 +1038,7 @@ def _matvec(self, x: NDArray) -> NDArray: self.trav_srcs, self.trav_recs, ) - elif not self.dynamic and not self.travsrcrec: + elif not self.travsrcrec: inputs = (x.ravel(), y, self.nsnr, self.nt, self.ni, self.itrav, self.travd) y = self._kirch_matvec(*inputs) @@ -1345,32 +1062,8 @@ def _rmatvec(self, x: NDArray) -> NDArray: self.vel, self.trav_srcs, self.trav_recs, - self.dist_srcs, - self.dist_recs, - self.aperture[0], - self.aperture[1], - self.aperturetap, - self.nz, - self.six, - self.rix, - self.angleaperture[0], - self.angleaperture[1], - self.angle_srcs, - self.angle_recs, - self.snell[0], - self.snell[1], - self.maxdist, - ) - elif self.dynamic and not self.travsrcrec: - inputs = ( - x, - y, - self.nsnr, - self.nt, - self.ni, - self.itrav, - self.travd, - self.amp, + self.amp_srcs, + self.amp_recs, self.aperture[0], self.aperture[1], self.aperturetap, @@ -1381,10 +1074,8 @@ def _rmatvec(self, x: NDArray) -> NDArray: self.angleaperture[1], self.angle_srcs, self.angle_recs, - self.snell[0], - self.snell[1], ) - elif not self.dynamic and self.travsrcrec: + elif self.travsrcrec: inputs = ( x, y, @@ -1396,7 +1087,7 @@ def _rmatvec(self, x: NDArray) -> NDArray: self.trav_srcs, self.trav_recs, ) - elif not self.dynamic and not self.travsrcrec: + elif not self.travsrcrec: inputs = (x, y, self.nsnr, self.nt, self.ni, self.itrav, self.travd) y = self._kirch_rmatvec(*inputs) diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 00000000..d2f6c854 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,52 @@ +[build-system] +requires = [ + "setuptools >= 65", + "setuptools_scm[toml]", + "wheel", +] +build-backend = "setuptools.build_meta" + +[project] +name = "pylops" +description = "Python library implementing linear operators to allow solving large-scale optimization problems" +readme = "README.md" +authors = [ + {name = "Matteo Ravasi", email = "matteoravasi@gmail.com"}, +] +license = {file = "LICENSE.md"} +keywords = ["algebra", "inverse problems", "large-scale optimization"] +classifiers = [ + "Development Status :: 5 - Production/Stable", + "Intended Audience :: Developers", + "Intended Audience :: Science/Research", + "Intended Audience :: Education", + "License :: OSI Approved :: GNU Lesser General Public License v3 (LGPLv3)", + "Natural Language :: English", + "Operating System :: OS Independent", + "Programming Language :: Python :: 3 :: Only", + "Programming Language :: Python :: 3.8", + "Programming Language :: Python :: 3.9", + "Programming Language :: Python :: 3.10", + "Topic :: Scientific/Engineering :: Mathematics", +] +dependencies = [ + "numpy >= 1.21.0", + "scipy >= 1.4.0", +] +dynamic = ["version"] + +[project.optional-dependencies] +advanced = [ + "llvmlite", + "numba", + "pyfftw", + "PyWavelets", + "scikit-fmm", + "spgl1", +] + +[tool.setuptools.packages.find] +exclude = ["pytests"] + +[tool.setuptools_scm] +version_file = "pylops/version.py" diff --git a/pytests/test_basicoperators.py b/pytests/test_basicoperators.py index 02ae6946..8f3f0528 100755 --- a/pytests/test_basicoperators.py +++ b/pytests/test_basicoperators.py @@ -63,7 +63,9 @@ def test_LinearRegression(par): assert dottest(LRop, par["ny"], 2) x = np.array([1.0, 2.0], dtype=np.float32) - xlsqr = lsqr(LRop, LRop * x, damp=1e-10, iter_lim=300, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr( + LRop, LRop * x, damp=1e-10, iter_lim=300, atol=1e-8, btol=1e-8, show=0 + )[0] assert_array_almost_equal(x, xlsqr, decimal=3) y = LRop * x @@ -82,7 +84,9 @@ def test_MatrixMult(par): assert dottest(Gop, par["ny"], par["nx"], complexflag=0 if par["imag"] == 0 else 3) x = np.ones(par["nx"]) + par["imag"] * np.ones(par["nx"]) - xlsqr = lsqr(Gop, Gop * x, damp=1e-20, iter_lim=300, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr(Gop, Gop * x, damp=1e-20, iter_lim=300, atol=1e-8, btol=1e-8, show=0)[ + 0 + ] assert_array_almost_equal(x, xlsqr, decimal=4) @@ -100,7 +104,9 @@ def test_MatrixMult_sparse(par): assert dottest(Gop, par["ny"], par["nx"], complexflag=0 if par["imag"] == 1 else 3) x = np.ones(par["nx"]) + par["imag"] * np.ones(par["nx"]) - xlsqr = lsqr(Gop, Gop * x, damp=1e-20, iter_lim=300, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr(Gop, Gop * x, damp=1e-20, iter_lim=300, atol=1e-8, btol=1e-8, show=0)[ + 0 + ] assert_array_almost_equal(x, xlsqr, decimal=4) @@ -133,7 +139,9 @@ def test_MatrixMult_repeated(par): ) x = (np.ones((par["nx"], 5)) + par["imag"] * np.ones((par["nx"], 5))).ravel() - xlsqr = lsqr(Gop, Gop * x, damp=1e-20, iter_lim=300, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr(Gop, Gop * x, damp=1e-20, iter_lim=300, atol=1e-8, btol=1e-8, show=0)[ + 0 + ] assert_array_almost_equal(x, xlsqr, decimal=4) @@ -470,6 +478,45 @@ def test_Sum2D(par): assert dottest(Sop, np.prod(dim_d), par["ny"] * par["nx"]) +@pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j), (par3)]) +def test_Sum2D_forceflat(par): + """Dot-test for Sum operator on 2d signal with forceflat""" + np.random.seed(10) + flat_dimsd = par["ny"] + flat_dims = par["ny"] * par["nx"] + x = np.random.randn(flat_dims) + par["imag"] * np.random.randn(flat_dims) + + Sop_True = Sum((par["ny"], par["nx"]), axis=-1, forceflat=True) + y = Sop_True @ x + xadj = Sop_True.H @ y + assert y.shape == (flat_dimsd,) + assert xadj.shape == (flat_dims,) + + Sop_None = Sum((par["ny"], par["nx"]), axis=-1) + y = Sop_None @ x + xadj = Sop_None.H @ y + assert y.shape == (par["ny"],) + assert xadj.shape == (par["ny"], par["nx"]) + + Sop_False = Sum((par["ny"], par["nx"]), axis=-1, forceflat=False) + y = Sop_False @ x + xadj = Sop_False.H @ y + assert y.shape == (par["ny"],) + assert xadj.shape == (par["ny"], par["nx"]) + + with pytest.raises(ValueError): + Sop_True * Sop_False.H + + Sop = Sop_True * Sop_None.H + assert Sop.forceflat is True + + Sop = Sop_False * Sop_None.H + assert Sop.forceflat is False + + Sop = Sop_None * Sop_None.H + assert Sop.forceflat is None + + @pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j), (par3)]) def test_Sum3D(par): """Dot-test, forward and adjoint for Sum operator on 3d signal""" diff --git a/pytests/test_convolve.py b/pytests/test_convolve.py index 4cddec63..d26938ba 100755 --- a/pytests/test_convolve.py +++ b/pytests/test_convolve.py @@ -133,7 +133,9 @@ def test_Convolve1D(par): x = np.zeros((par["nx"])) x[par["nx"] // 2] = 1.0 - xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr( + Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0 + )[0] assert_array_almost_equal(x, xlsqr, decimal=1) # 1D on 2D @@ -153,7 +155,9 @@ def test_Convolve1D(par): int(par["nx"] / 2 - 3) : int(par["nx"] / 2 + 3), ] = 1.0 x = x.ravel() - xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr( + Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0 + )[0] assert_array_almost_equal(x, xlsqr, decimal=1) # 1D on 3D @@ -175,10 +179,31 @@ def test_Convolve1D(par): int(par["nx"] / 2 - 3) : int(par["nx"] / 2 + 3), ] = 1.0 x = x.ravel() - xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0)[ + 0 + ] assert_array_almost_equal(x, xlsqr, decimal=1) +@pytest.mark.parametrize( + "par", [(par1_1d), (par2_1d), (par3_1d), (par4_1d), (par5_1d), (par6_1d)] +) +def test_Convolve1D_long(par): + """Dot-test and inversion for Convolve1D operator with long filter""" + np.random.seed(10) + # 1D + if par["axis"] == 0: + x = np.zeros((par["nx"])) + x[par["nx"] // 2] = 1.0 + Xop = Convolve1D(nfilt[0], h=x, offset=nfilt[0] // 2, dtype="float64") + assert dottest(Xop, par["nx"], nfilt[0]) + + h1lsqr = lsqr( + Xop, Xop * h1, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0 + )[0] + assert_array_almost_equal(h1, h1lsqr, decimal=1) + + @pytest.mark.parametrize( "par", [(par1_2d), (par2_2d), (par3_2d), (par4_2d), (par5_2d), (par6_2d)] ) @@ -200,7 +225,9 @@ def test_Convolve2D(par): int(par["nx"] / 2 - 3) : int(par["nx"] / 2 + 3), ] = 1.0 x = x.ravel() - xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr( + Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0 + )[0] assert_array_almost_equal(x, xlsqr, decimal=1) # 2D on 3D @@ -224,7 +251,9 @@ def test_Convolve2D(par): int(par["nx"] / 2 - 3) : int(par["nx"] / 2 + 3), ] = 1.0 x = x.ravel() - xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0)[ + 0 + ] # due to ringing in solution we cannot use assert_array_almost_equal assert np.linalg.norm(xlsqr - x) / np.linalg.norm(xlsqr) < 2e-1 diff --git a/pytests/test_interpolation.py b/pytests/test_interpolation.py index a691ba16..a22c04ef 100755 --- a/pytests/test_interpolation.py +++ b/pytests/test_interpolation.py @@ -405,6 +405,30 @@ def test_Bilinear_2dsignal(par): assert_array_almost_equal(y, x[iava[0], iava[1]]) +@pytest.mark.parametrize("par", [(par1), (par2)]) +def test_Bilinear_2dsignal_flatten(par): + """Dot-test and forward for Interp operator for 2d signal with forceflat""" + np.random.seed(1) + dims = (par["nx"], par["nt"]) + flat_dims = par["nx"] * par["nt"] + dimsd = 10 + + x = np.random.normal(0, 1, dims) + par["imag"] * np.random.normal(0, 1, dims) + + iava = np.vstack((np.arange(0, dimsd), np.arange(0, dimsd))) + Iop_True = Bilinear(iava, dims=dims, dtype=par["dtype"], forceflat=True) + y = Iop_True @ x + xadj = Iop_True.H @ y + assert y.shape == (dimsd,) + assert xadj.shape == (flat_dims,) + + Iop_None = Bilinear(iava, dims=dims, dtype=par["dtype"]) + y = Iop_None @ x + xadj = Iop_None.H @ y + assert y.shape == (dimsd,) + assert xadj.shape == dims + + @pytest.mark.parametrize("par", [(par1), (par2)]) def test_Bilinear_3dsignal(par): """Dot-test and forward for Interp operator for 3d signal""" diff --git a/pytests/test_kirchhoff.py b/pytests/test_kirchhoff.py index 02deeafd..d601d91b 100755 --- a/pytests/test_kirchhoff.py +++ b/pytests/test_kirchhoff.py @@ -49,7 +49,7 @@ r2d = np.vstack((rx, 2 * np.ones(PAR["nrx"]))) r3d = np.vstack((ryy.ravel(), rxx.ravel(), 2 * np.ones(PAR["nrx"] * PAR["nry"]))) -wav, _, wavc = ricker(t[:41], f0=40) +wav, _, wavc = ricker(t[:21], f0=40) par1 = {"mode": "analytic", "dynamic": False} par2 = {"mode": "eikonal", "dynamic": False} @@ -287,7 +287,7 @@ def test_kirchhoff3d(par): ) def test_kirchhoff2d_trav_vs_travsrcrec(par): """Compare 2D Kirchhoff operator forward and adjoint when using trav (original behavior) - or trav_src and trav_rec (new reccomended behaviour)""" + or trav_src and trav_rec (new recommended behaviour)""" # new behaviour Dop = Kirchhoff( @@ -311,21 +311,7 @@ def test_kirchhoff2d_trav_vs_travsrcrec(par): ) + Dop.trav_recs.reshape(PAR["nx"] * PAR["nz"], 1, PAR["nrx"]) trav = trav.reshape(PAR["nx"] * PAR["nz"], PAR["nsx"] * PAR["nrx"]) if par["dynamic"]: - dist = Dop.dist_srcs.reshape( - PAR["nx"] * PAR["nz"], PAR["nsx"], 1 - ) + Dop.dist_recs.reshape(PAR["nx"] * PAR["nz"], 1, PAR["nrx"]) - dist = dist.reshape(PAR["nx"] * PAR["nz"], PAR["nsx"] * PAR["nrx"]) - - cosangle = np.cos(Dop.angle_srcs).reshape( - PAR["nx"] * PAR["nz"], PAR["nsx"], 1 - ) + np.cos(Dop.angle_recs).reshape(PAR["nx"] * PAR["nz"], 1, PAR["nrx"]) - cosangle = cosangle.reshape(PAR["nx"] * PAR["nz"], PAR["nsx"] * PAR["nrx"]) - - epsdist = 1e-2 - amp = 1 / (dist + epsdist * np.max(dist)) - - amp *= np.abs(cosangle) - amp /= v0 + amp = (Dop.amp_srcs, Dop.amp_recs) D1op = Kirchhoff( z, @@ -361,7 +347,7 @@ def test_kirchhoff2d_trav_vs_travsrcrec(par): ) def test_kirchhoff3d_trav_vs_travsrcrec(par): """Compare 3D Kirchhoff operator forward and adjoint when using trav (original behavior) - or trav_src and trav_rec (new reccomended behaviour)""" + or trav_src and trav_rec (new recommended behaviour)""" # new behaviour Dop = Kirchhoff( diff --git a/pytests/test_lsm.py b/pytests/test_lsm.py index 33024710..7f39cefd 100755 --- a/pytests/test_lsm.py +++ b/pytests/test_lsm.py @@ -48,7 +48,7 @@ r2d = np.vstack((rx, 2 * np.ones(PAR["nrx"]))) r3d = np.vstack((ryy.ravel(), rxx.ravel(), 2 * np.ones(PAR["nrx"] * PAR["nry"]))) -wav, _, wavc = ricker(t[:41], f0=40) +wav, _, wavc = ricker(t[:21], f0=40) par1 = {"mode": "analytic", "dynamic": False} par2 = {"mode": "eikonal", "dynamic": False} diff --git a/pytests/test_metrics.py b/pytests/test_metrics.py index 7ca20f67..420434af 100755 --- a/pytests/test_metrics.py +++ b/pytests/test_metrics.py @@ -49,7 +49,7 @@ def test_psnr(par): xref = np.ones(par["nx"]) xcmp = np.zeros(par["nx"]) - psnrsame = psnr(xref, xref, xmax=1.0) - psnrcmp = psnr(xref, xcmp, xmax=1.0) + psnrsame = psnr(xref, xref, xmax=1.0, xmin=0.0) + psnrcmp = psnr(xref, xcmp, xmax=1.0, xmin=0.0) assert psnrsame == np.inf assert psnrcmp == 0.0 diff --git a/pytests/test_nonstatconvolve.py b/pytests/test_nonstatconvolve.py index 0a90c489..688fe39d 100755 --- a/pytests/test_nonstatconvolve.py +++ b/pytests/test_nonstatconvolve.py @@ -7,8 +7,10 @@ from pylops.signalprocessing import ( Convolve1D, Convolve2D, + ConvolveND, NonStationaryConvolve1D, NonStationaryConvolve2D, + NonStationaryConvolve3D, NonStationaryFilters1D, NonStationaryFilters2D, ) @@ -16,10 +18,18 @@ # filters nfilts = (5, 7) +nfilts3 = (5, 5, 7) + h1 = triang(nfilts[0], sym=True) h2 = np.outer(triang(nfilts[0], sym=True), triang(nfilts[1], sym=True)) +h3 = np.outer( + triang(nfilts[0], sym=True), + np.outer(triang(nfilts[0], sym=True), triang(nfilts[1], sym=True)[np.newaxis]).T, +).reshape(nfilts3) + h1stat = np.vstack([h1, h1, h1]) h1ns = np.vstack([h1, -h1, 2 * h1]) + h2stat = np.vstack( [ h2.ravel(), @@ -41,6 +51,26 @@ ] ).reshape(3, 2, nfilts[0], nfilts[1]) +h3stat = np.vstack( + [ + h3.ravel(), + ] + * 8 +).reshape(2, 2, 2, nfilts[0], nfilts[0], nfilts[1]) +h3ns = np.vstack( + [ + 2 * h3.ravel(), + h3.ravel(), + h3.ravel(), + h3.ravel(), + h3.ravel(), + h3.ravel(), + -h3.ravel(), + 2 * h3.ravel(), + ] +).reshape(2, 2, 2, nfilts[0], nfilts[0], nfilts[1]) + + par1_1d = { "nz": 21, "nx": 31, @@ -144,7 +174,9 @@ def test_NonStationaryConvolve1D(par): x = np.zeros((par["nx"])) x[par["nx"] // 2] = 1.0 - xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr( + Cop, Cop * x, damp=1e-20, iter_lim=200, atol=1e-8, btol=1e-8, show=0 + )[0] assert_array_almost_equal(x, xlsqr, decimal=1) # 1D on 2D @@ -164,7 +196,9 @@ def test_NonStationaryConvolve1D(par): int(par["nz"] / 2 - 3) : int(par["nz"] / 2 + 3), ] = 1.0 x = x.ravel() - xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=400, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=400, atol=1e-8, btol=1e-8, show=0)[ + 0 + ] assert_array_almost_equal(x, xlsqr, decimal=1) @@ -208,7 +242,9 @@ def test_NonStationaryConvolve2D(par): int(par["nz"] / 2 - 3) : int(par["nz"] / 2 + 3), ] = 1.0 x = x.ravel() - xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=400, atol=1e-8, btol=1e-8, show=0)[0] + xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=400, atol=1e-8, btol=1e-8, show=0)[ + 0 + ] assert_array_almost_equal(x, xlsqr, decimal=1) @@ -240,7 +276,7 @@ def test_StationaryConvolve2D(par): (par1_1d), ], ) -def test_NonStationaryFilters2D(par): +def test_NonStationaryFilters1D(par): """Dot-test and inversion for NonStationaryFilters2D operator""" x = np.zeros((par["nx"])) x[par["nx"] // 4], x[par["nx"] // 2], x[3 * par["nx"] // 4] = 1.0, 1.0, 1.0 @@ -275,3 +311,53 @@ def test_NonStationaryFilters2D(par): h2lsqr = lsqr(Cop, Cop * h2ns.ravel(), damp=1e-20, iter_lim=400, show=0)[0] assert_array_almost_equal(h2ns.ravel(), h2lsqr, decimal=1) + + +@pytest.mark.parametrize("par", [(par_2d)]) +def test_NonStationaryConvolve3D(par): + """Dot-test and inversion for NonStationaryConvolve3D operator""" + Cop = NonStationaryConvolve3D( + dims=(par["nx"], par["nx"], par["nz"]), + hs=h3ns, + ihx=(int(par["nx"] // 4), int(3 * par["nx"] // 4)), + ihy=(int(par["nx"] // 4), int(3 * par["nx"] // 4)), + ihz=(int(par["nz"] // 4), int(3 * par["nz"] // 4)), + dtype="float64", + ) + assert dottest( + Cop, par["nx"] * par["nx"] * par["nz"], par["nx"] * par["nx"] * par["nz"] + ) + + x = np.zeros((par["nx"], par["nx"], par["nz"])) + x[ + int(par["nx"] / 2 - 3) : int(par["nx"] / 2 + 3), + int(par["nx"] / 2 - 3) : int(par["nx"] / 2 + 3), + int(par["nz"] / 2 - 3) : int(par["nz"] / 2 + 3), + ] = 1.0 + x = x.ravel() + xlsqr = lsqr(Cop, Cop * x, damp=1e-20, iter_lim=40, atol=1e-8, btol=1e-8, show=0)[0] + # given the size of the problem, we can only run few iterations and test accuracy up to 30% + assert np.linalg.norm(x - xlsqr) / np.linalg.norm(x) < 0.3 + + +@pytest.mark.parametrize("par", [(par_2d)]) +def test_StationaryConvolve3D(par): + """Check that Convolve3D and NonStationaryConvolve3D return same result for + stationary filter""" + Cop = NonStationaryConvolve3D( + dims=(par["nx"], par["nx"], par["nz"]), + hs=h3stat, + ihx=(int(par["nx"] // 4), int(3 * par["nx"] // 4)), + ihy=(int(par["nx"] // 4), int(3 * par["nx"] // 4)), + ihz=(int(par["nz"] // 4), int(3 * par["nz"] // 4)), + dtype="float64", + ) + Cop_stat = ConvolveND( + dims=(par["nx"], par["nx"], par["nz"]), + h=h3, + offset=(nfilts[0] // 2, nfilts[0] // 2, nfilts[1] // 2), + dtype="float64", + ) + x = np.random.normal(0, 1, (par["nx"], par["nx"], par["nz"])) + + assert_array_almost_equal(Cop_stat * x, Cop * x, decimal=10) diff --git a/pytests/test_signalutils.py b/pytests/test_signalutils.py index 79c6d4e8..d7551a2e 100755 --- a/pytests/test_signalutils.py +++ b/pytests/test_signalutils.py @@ -4,47 +4,89 @@ from pylops.utils.signalprocessing import convmtx, nonstationary_convmtx -par1 = {"nt": 51, "imag": 0, "dtype": "float32"} # real -par2 = {"nt": 51, "imag": 1j, "dtype": "complex64"} # complex +par1 = {"nt": 51, "nh": 7, "imag": 0, "dtype": "float32"} # odd sign, odd filt, real +par1j = { + "nt": 51, + "nh": 7, + "imag": 1j, + "dtype": "complex64", +} # odd sign, odd filt, complex +par2 = {"nt": 50, "nh": 7, "imag": 0, "dtype": "float32"} # even sign, odd filt, real +par2j = { + "nt": 50, + "nh": 7, + "imag": 1j, + "dtype": "complex64", +} # even sign, odd filt, complex +par3 = {"nt": 51, "nh": 6, "imag": 0, "dtype": "float32"} # odd sign, even filt, real +par3j = { + "nt": 51, + "nh": 6, + "imag": 1j, + "dtype": "complex64", +} # odd sign, even filt, complex +par4 = {"nt": 50, "nh": 6, "imag": 0, "dtype": "float32"} # even sign, even filt, real +par4j = { + "nt": 50, + "nh": 6, + "imag": 1j, + "dtype": "complex64", +} # even sign, even filt, complex + np.random.seed(10) -@pytest.mark.parametrize("par", [(par1), (par2)]) +@pytest.mark.parametrize("par", [(par1), (par1j), (par2), (par2j)]) def test_convmtx(par): - """Compare convmtx with np.convolve""" + """Compare convmtx with np.convolve (small filter)""" + x = np.random.normal(0, 1, par["nt"]) + par["imag"] * np.random.normal( + 0, 1, par["nt"] + ) + + h = np.hanning(par["nh"]) + H = convmtx(h, par["nt"], par["nh"] // 2) + + y = np.convolve(x, h, mode="same") + y1 = np.dot(H[: par["nt"]], x) + assert_array_almost_equal(y, y1, decimal=4) + + +@pytest.mark.parametrize("par", [(par1), (par1j), (par2), (par2j)]) +def test_convmtx1(par): + """Compare convmtx with np.convolve (large filter)""" x = np.random.normal(0, 1, par["nt"]) + par["imag"] * np.random.normal( 0, 1, par["nt"] ) - nh = 7 - h = np.hanning(7) - H = convmtx(h, par["nt"]) - H = H[:, nh // 2 : -nh // 2 + 1] + h = np.hanning(par["nh"]) + X = convmtx( + x, par["nh"], par["nh"] // 2 - 1 if par["nh"] % 2 == 0 else par["nh"] // 2 + ) y = np.convolve(x, h, mode="same") - y1 = np.dot(H, x) + y1 = np.dot(X[: par["nt"]], h) assert_array_almost_equal(y, y1, decimal=4) -@pytest.mark.parametrize("par", [(par1), (par2)]) +@pytest.mark.parametrize("par", [(par1), (par1j)]) def test_nonstationary_convmtx(par): """Compare nonstationary_convmtx with convmtx for stationary filter""" x = np.random.normal(0, 1, par["nt"]) + par["imag"] * np.random.normal( 0, 1, par["nt"] ) - nh = 7 - h = np.hanning(7) - H = convmtx(h, par["nt"]) - H = H[:, nh // 2 : -nh // 2 + 1] + h = np.hanning(par["nh"]) + H = convmtx( + h, par["nt"], par["nh"] // 2 - 1 if par["nh"] % 2 == 0 else par["nh"] // 2 + ) H1 = nonstationary_convmtx( np.repeat(h[:, np.newaxis], par["nt"], axis=1).T, par["nt"], - hc=nh // 2, + hc=par["nh"] // 2, pad=(par["nt"], par["nt"]), ) - y = np.dot(H, x) + y = np.dot(H[: par["nt"]], x) y1 = np.dot(H1, x) assert_array_almost_equal(y, y1, decimal=4) diff --git a/pytests/test_solver.py b/pytests/test_solver.py index 6030bf99..ba32f79b 100644 --- a/pytests/test_solver.py +++ b/pytests/test_solver.py @@ -94,7 +94,7 @@ def test_cg(par): "par", [(par1), (par2), (par3), (par4), (par1j), (par2j), (par3j), (par3j)] ) def test_cg_ndarray(par): - """CG with linear operator""" + """CG with linear operator (and ndarray as input and output)""" np.random.seed(10) dims = dimsd = (par["nx"], par["ny"]) @@ -119,6 +119,35 @@ def test_cg_ndarray(par): assert_array_almost_equal(x, xinv, decimal=4) +@pytest.mark.parametrize( + "par", [(par1), (par2), (par3), (par4), (par1j), (par2j), (par3j), (par3j)] +) +def test_cg_forceflat(par): + """CG with linear operator (and forced 1darray as input and output)""" + np.random.seed(10) + + dims = dimsd = (par["nx"], par["ny"]) + x = np.ones(dims) + par["imag"] * np.ones(dims) + + A = np.random.normal(0, 10, (x.size, x.size)) + par["imag"] * np.random.normal( + 0, 10, (x.size, x.size) + ) + A = np.conj(A).T @ A # to ensure definite positive matrix + Aop = MatrixMult(A, dtype=par["dtype"], forceflat=True) + Aop.dims = dims + Aop.dimsd = dimsd + + if par["x0"]: + x0 = np.random.normal(0, 10, dims) + par["imag"] * np.random.normal(0, 10, dims) + else: + x0 = None + + y = Aop * x + xinv = cg(Aop, y, x0=x0, niter=2 * x.size, tol=1e-5, show=True)[0] + assert xinv.shape == x.ravel().shape + assert_array_almost_equal(x.ravel(), xinv, decimal=4) + + @pytest.mark.parametrize( "par", [(par1), (par2), (par3), (par4), (par1j), (par2j), (par3j), (par3j)] ) diff --git a/readthedocs.yml b/readthedocs.yml deleted file mode 100644 index a4cf3385..00000000 --- a/readthedocs.yml +++ /dev/null @@ -1,14 +0,0 @@ -# .readthedocs.yml - -version: 2 - -build: - os: ubuntu-20.04 - tools: - python: "3.9" - -python: - install: - - requirements: requirements-dev.txt - - method: setuptools - path: . diff --git a/setup.cfg b/setup.cfg index 1691b3bc..4c13f09c 100755 --- a/setup.cfg +++ b/setup.cfg @@ -11,7 +11,6 @@ per-file-ignores = __init__.py: F401, F403, F405 max-line-length = 88 - # mypy global options [mypy] plugins = numpy.typing.mypy_plugin diff --git a/setup.py b/setup.py deleted file mode 100755 index 8b82afa2..00000000 --- a/setup.py +++ /dev/null @@ -1,54 +0,0 @@ -import os - -from setuptools import find_packages, setup - - -def src(pth): - return os.path.join(os.path.dirname(__file__), pth) - - -# Project description -descr = ( - "Python library implementing linear operators to allow solving large-scale optimization " - "problems without requiring to explicitly create a dense (or sparse) matrix." -) - -# Setup -setup( - name="pylops", - description=descr, - long_description=open(src("README.md")).read(), - long_description_content_type="text/markdown", - keywords=["algebra", "inverse problems", "large-scale optimization"], - classifiers=[ - "Development Status :: 5 - Production/Stable", - "Intended Audience :: Developers", - "Intended Audience :: Science/Research", - "License :: OSI Approved :: GNU Lesser General Public License v2 or later (LGPLv2+)", - "Natural Language :: English", - "Programming Language :: Python :: 3.8", - "Programming Language :: Python :: 3.9", - "Topic :: Scientific/Engineering :: Mathematics", - ], - author="mrava", - author_email="matteoravasi@gmail.com", - install_requires=["numpy >= 1.21.0", "scipy >= 1.4.0"], - extras_require={ - "advanced": [ - "llvmlite", - "numba", - "pyfftw", - "PyWavelets", - "scikit-fmm", - "spgl1", - ] - }, - packages=find_packages(exclude=["pytests"]), - use_scm_version=dict( - root=".", relative_to=__file__, write_to=src("pylops/version.py") - ), - setup_requires=["pytest-runner", "setuptools_scm"], - test_suite="pytests", - tests_require=["pytest"], - zip_safe=True, -) diff --git a/tutorials/ctscan.py b/tutorials/ctscan.py index 7869593d..9017499f 100755 --- a/tutorials/ctscan.py +++ b/tutorials/ctscan.py @@ -28,11 +28,17 @@ # such a type: # # .. math:: -# t(r,\theta; x) = \tan(90°-\theta)x + \frac{r}{\sin(\theta)} +# t\sin(\theta) + x\cos(\theta) = r, # # where :math:`\theta` is the angle between the x-axis (:math:`x`) and # the perpendicular to the summation line and :math:`r` is the distance -# from the origin of the summation line. +# from the origin of the summation line. Radon transform in CT +# corresponds to the integral of the input image along the straight line above. +# To implement the integration in PyLops we simply need to express +# :math:`t(r,\theta;x)` which is given by: +# +# .. math:: +# t(r,\theta; x) = \tan\left(\frac{\pi}{2}-\theta\right)x + \frac{r}{\sin(\theta)}. @jit(nopython=True) @@ -44,6 +50,13 @@ def radoncurve(x, r, theta): ) +############################################################################### +# Note that in the above implementation we added centering :math:`t \mapsto t - n_y/2` and +# :math:`r \mapsto r - n_y/2` so that the origin of the integration lines is exactly in the +# center of the image (centering for :math:`x` is not needed because we will use +# ``centeredh=True`` in the constructor of ``Radon2D``). + + x = np.load("../testdata/optimization/shepp_logan_phantom.npy").T x = x / x.max() nx, ny = x.shape @@ -81,14 +94,36 @@ def radoncurve(x, r, theta): axs[1].set_title("Data") axs[1].axis("tight") axs[2].imshow(xrec.T, cmap="gray") -axs[2].set_title("Adjoint model") +axs[2].set_title("Adjoint model in PyLops") axs[2].axis("tight") fig.tight_layout() +############################################################################### +# Note that our raw data ``y`` *does not represent exactly* classical sinograms +# in medical imaging. Integration along curves in the adjoint form of +# :func:`pylops.signalprocessing.Radon2D` is performed with respect to +# :math:`dx`, whereas canonically it is assumed to be with respect to the natural +# parametrization :math:`dl = \sqrt{(dx)^2 + (dt)^2}`. To retrieve back the +# classical sinogram we have to divide data by the jacobian +# :math:`j(x,l) = \left\vert dx/dl \right\vert = |\sin(\theta)|`. + +sinogram = np.divide( + y.T, np.abs(np.sin(theta) + 1e-15) +) # small shift to avoid zero-division +fig, axs = plt.subplots(1, 2, figsize=(10, 4)) +axs[0].imshow(y.T, cmap="gray") +axs[0].set_title("Data") +axs[0].axis("tight") +axs[1].imshow(sinogram, cmap="gray") +axs[1].set_title("Sinogram in medical imaging") +axs[1].axis("tight") +fig.tight_layout() ############################################################################### -# 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. +# From now on, we will not pursue further working with the "true sinogram", instead +# we will reconstruct the original phantom directly from ``y``. For this 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 diff --git a/tutorials/linearoperator.py b/tutorials/linearoperator.py index e1b86669..6c85f9c4 100755 --- a/tutorials/linearoperator.py +++ b/tutorials/linearoperator.py @@ -5,11 +5,17 @@ library for both new users and developers. 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 +different ways to apply the forward and adjoint operations. We will then 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, +have been implemented for such a class and will allow summing, subtracting, chaining, etc. multiple operators in very easy and expressive way. + +Finally, we will consider a scenario where the input and output vectors of an +operator are naturally represented by nd-arrays. Since pylops V2, we can both +pass to the operator a flattened array (original behaviour) as well as the +nd-array with shape `dims` (new behaviour); according to input, the output +will be either a vector or a nd-array with shape `dimsd`. """ ############################################################################### @@ -112,10 +118,10 @@ cmd2 = "Dop.rmatvec(x)" # .H* (pre-computed H) -cmd3 = "DopH*x" +cmd3 = "DopH@x" # .H* -cmd4 = "Dop.H*x" +cmd4 = "Dop.H@x" # timing t1 = 1.0e3 * np.array(timeit.repeat(cmd1, setup=cmd_setup, number=500, repeat=5)) @@ -171,7 +177,7 @@ print(Dop**3) # * and / -y = Dop * x +y = Dop @ x print(Dop / y) # eigs @@ -194,14 +200,14 @@ x = np.ones(n) Dop = pylops.Diagonal(d) -print(f"y = Dx = {Dop * x}") -print(f"y = conj(D)x = {Dop.conj() * x}") +print(f"y = Dx = {Dop @ x}") +print(f"y = conj(D)x = {Dop.conj() @ x}") ############################################################################### -# At this point, the concept of linear operator may sound abstract. -# The convinience method :func:`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 +# The concept of linear operator may sound abstract. The convenience method +# :func:`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)) @@ -212,14 +218,14 @@ plt.tight_layout() ############################################################################### -# At this point it is worth reiterating that if two linear operators are +# 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 :py:class:`pylops.LinearOperator` operator. This means # that we can still apply any of the methods implemented in our class # like ``*`` or ``/``. Dop1 = Dop - Dop.conj() -y = Dop1 * x +y = Dop1 @ x print(f"x = (Dop - conj(Dop))/y = {Dop1 / y}") D1 = Dop1.todense() @@ -232,7 +238,7 @@ plt.tight_layout() ############################################################################### -# Finally, another important feature of PyLops linear operators is that we can +# 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 @@ -252,6 +258,27 @@ print(f"Forward evaluations: {Dop.matvec_count}") print(f"Adjoint evaluations: {Dop.rmatvec_count}") + +######################################################################## +# Finally, let's consider a case where the :py:class:`pylops.basicoperators.Diagonal` +# operator acts on one dimension of a 2d-array. We will see how in this case the forward +# and adjoint passes can be applied on the flattened array as well as on the 2d-array directly + +# +m, n = 10, 5 +d = np.arange(n) + 1.0 +x = np.ones((m, n)) +Dop = pylops.Diagonal(d, dims=(m, n), axis=-1) + +yflat = Dop @ x.ravel() +y = Dop @ x + +xadjflat = Dop.H @ yflat +xadj = Dop.H @ y + +print(f"yflat shape = {yflat.shape}, y shape = {y.shape}") +print(f"xadjflat shape = {xadjflat.shape}, xadj shape = {xadj.shape}") + ############################################################################### # This first tutorial is completed. You have seen the basic operations that # can be performed using :py:class:`pylops.LinearOperator` and you diff --git a/tutorials/solvers.py b/tutorials/solvers.py index 018b5f85..c0bfe72c 100755 --- a/tutorials/solvers.py +++ b/tutorials/solvers.py @@ -178,8 +178,8 @@ # derivative operator # # .. math:: -# \mathbf{x} = (\mathbf{R^TR}+\epsilon_\nabla^2\nabla^T\nabla)^{-1} -# \mathbf{R^Ty} +# \mathbf{x} = (\mathbf{R}^T\mathbf{R}+\epsilon_\nabla^2\nabla^T\nabla)^{-1} +# \mathbf{R}^T\mathbf{y} # Create regularization operator D2op = pylops.SecondDerivative(N, dtype="float64")