diff --git a/.github/workflows/build.yaml b/.github/workflows/build.yaml index abf7ff7e..09795a41 100644 --- a/.github/workflows/build.yaml +++ b/.github/workflows/build.yaml @@ -7,7 +7,7 @@ jobs: strategy: matrix: platform: [ ubuntu-latest, macos-latest ] - python-version: ["3.8", "3.9", "3.10"] + python-version: ["3.8", "3.9", "3.10", "3.11"] runs-on: ${{ matrix.platform }} steps: diff --git a/.github/workflows/codacy-coverage-reporter.yaml b/.github/workflows/codacy-coverage-reporter.yaml index 5e3ed9fe..2c81e290 100644 --- a/.github/workflows/codacy-coverage-reporter.yaml +++ b/.github/workflows/codacy-coverage-reporter.yaml @@ -8,12 +8,16 @@ jobs: build: strategy: matrix: - platform: [ ubuntu-latest, macos-latest ] - python-version: ["3.8", "3.9", "3.10"] + platform: [ ubuntu-latest, ] + python-version: ["3.8", ] 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: 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/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/Makefile b/Makefile index 2db298d1..341c53be 100755 --- a/Makefile +++ b/Makefile @@ -29,6 +29,9 @@ install_conda: dev-install_conda: conda env create -f environment-dev.yml && conda activate pylops && pip install -e . +dev-install_conda_arm: + conda env create -f environment-dev-arm.yml && conda activate pylops && pip install -e . + tests: make pythoncheck pytest diff --git a/README.md b/README.md index 3d1fc89f..17512736 100644 --- a/README.md +++ b/README.md @@ -149,3 +149,4 @@ A list of video tutorials to learn more about PyLops: * Rohan Babbar, rohanbabbar04 * Wei Zhang, ZhangWeiGeo * Fedor Goncharov, fedor-goncharov +* Alex Rakowski, alex-rakowski diff --git a/docs/source/adding.rst b/docs/source/adding.rst index 87d68bea..8582cf46 100755 --- a/docs/source/adding.rst +++ b/docs/source/adding.rst @@ -4,9 +4,9 @@ Implementing new operators ========================== Users are welcome to create new operators and add them to the PyLops library. -In this tutorial, we will go through the key steps in the definition of an operator, using the -:py:class:`pylops.Diagonal` as an example. This is a very simple operator that applies a diagonal matrix to the model -in forward mode and to the data in adjoint mode. +In this tutorial, we will go through the key steps in the definition of an operator, using a simplified version of the +:py:class:`pylops.Diagonal` operator as an example. This is a very simple operator that applies a diagonal matrix +to the model in forward mode and to the data in adjoint mode. Creating the operator @@ -45,14 +45,17 @@ Initialization (``__init__``) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ We then need to create the ``__init__`` where the input parameters are passed and saved as members of our class. -While the input parameters change from operator to operator, it is always required to create three members, the first -called ``shape`` with a tuple containing the dimensions of the operator in the data and model space, the second -called ``dtype`` with the data type object (:obj:`np.dtype`) of the model and data, and the third -called ``explicit`` with a boolean (``True`` or ``False``) identifying if the operator can be inverted by a direct -solver or requires an iterative solver. This member is ``True`` if the operator has also a member ``A`` that contains -the matrix to be inverted like for example in the :py:class:`pylops.MatrixMult` operator, and it will be ``False`` otherwise. -In this case we have another member called ``d`` which is equal to the input vector containing the diagonal elements -of the matrix we want to multiply to the model and data. +While the input parameters change from operator to operator, it is always required to create three members: + +- ``dtype``: data type object (of type :obj:`str` or :obj:`np.dtype`) of the model and data; +- ``shape``: a tuple containing the dimensions of the operator in the data and model space; +- ``explicit``: a boolean (``True`` or ``False``) identifying if the operator can be inverted by a direct solver or + requires an iterative solver. This member is ``True`` if the operator has also a member ``A`` that contains + the matrix to be inverted like for example in the :py:class:`pylops.MatrixMult` operator, and it will be + ``False`` otherwise. + +In this specific case, we have another member called ``d`` which is equal to the input vector containing the diagonal +elements of the matrix we want to multiply to the model and data. .. code-block:: python @@ -62,7 +65,7 @@ of the matrix we want to multiply to the model and data. self.dtype = np.dtype(dtype) self.explicit = False -Alternatively, since version 2.0.0, the recommended way of initializing operators derived from the base +Alternatively, since version ``v2.0.0``, the recommended way of initializing operators derived from the base :py:class:`pylops.LinearOperator` class is to invoke ``super`` to assign the required attributes: .. code-block:: python @@ -72,8 +75,14 @@ Alternatively, since version 2.0.0, the recommended way of initializing operator super().__init__(dtype=np.dtype(dtype), shape=(len(self.d), len(self.d))) In this case, there is no need to declare ``explicit`` as it already defaults to ``False``. -Since version 2.0.0, every :py:class:`pylops.LinearOperator` class is imbued with ``dims``, -``dimsd``, ``clinear`` and ``explicit``, in addition to the required ``dtype`` and ``shape``. + +Moreover, since version ``v2.0.0``, every :py:class:`pylops.LinearOperator` class is imbued with ``dims``, +``dimsd``, and ``clinear`` in addition to the required ``dtype``, ``shape``, and ``explicit``. Note that +``dims`` and ``dimsd`` can be defined in spite of ``shape``, which will be automatically assigned within the +``super`` method: the main difference between ``dims``/``dimsd`` and ``shape`` is the the former variables can be +used the define the n-dimensional nature of the input of an operator, whilst the latter variable refers to their overall +shape when the input is flattened. + See the docs of :py:class:`pylops.LinearOperator` for more information about what these attributes mean. @@ -91,6 +100,10 @@ We will finally need to ``return`` the result of this operation: def _matvec(self, x): return self.d * x +Note that since version ``v2.0.0``, this method can be decorated by the decorator ``@reshaped``. As discussed in +more details in the decorator documentation, by adding such decorator the input ``x`` is initially reshaped into +a nd-array of shape ``dims``, fed to the actual code in ``_matvec`` and then flattened. + Adjoint mode (``_rmatvec``) ^^^^^^^^^^^^^^^^^^^^^^^^^^^ Finally we need to implement the *adjoint mode* in the method ``_rmatvec``. In other words, we will need to write @@ -106,6 +119,10 @@ different from operator to operator): And that's it, we have implemented our first linear operator! +Similar to ``_matvec``, since version ``v2.0.0``, this method can also be decorated by the decorator ``@reshaped``. +When doing so, the input ``x`` is initially reshaped into +a nd-array of shape ``dimsd``, fed to the actual code in ``_rmatvec`` and then flattened. + Testing the operator -------------------- Being able to write an operator is not yet a guarantee of the fact that the operator is correct, or in other words 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 2a5daf61..6310549d 100755 --- a/docs/source/credits.rst +++ b/docs/source/credits.rst @@ -20,4 +20,5 @@ Contributors * `Aniket Singh Rawat `_, dikwickley * `Rohan Babbar `_, rohanbabbar04 * `Wei Zhang `_, ZhangWeiGeo -* `Fedor Goncharov `_, fedor-goncharov \ No newline at end of file +* `Fedor Goncharov `_, fedor-goncharov +* `Alex Rakowski `_, alex-rakowski 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/faq.rst b/docs/source/faq.rst index 08f6a59c..ab057392 100755 --- a/docs/source/faq.rst +++ b/docs/source/faq.rst @@ -14,9 +14,9 @@ working with linear operators is indeed that you don't really need to access the of an operator. -**2. Can I have an older version of** ``cupy`` **or** ``cusignal`` **installed in my system (** ``cupy-cudaXX<8.1.0`` **or** ``cusignal>=0.16.0`` **)?** +**2. Can I have an older version of** ``cupy`` **installed in my system (** ``cupy-cudaXX<10.6.0``)?** Yes. Nevertheless you need to tell PyLops that you don't want to use its ``cupy`` -backend by setting the environment variable ``CUPY_PYLOPS=0`` or ``CUPY_SIGNAL=0``. +backend by setting the environment variable ``CUPY_PYLOPS=0``. Failing to do so will lead to an error when you import ``pylops`` because some of the ``cupyx`` routines that we use are not available in earlier version of ``cupy``. \ No newline at end of file diff --git a/docs/source/gpu.rst b/docs/source/gpu.rst index 7bc361d9..864451c4 100755 --- a/docs/source/gpu.rst +++ b/docs/source/gpu.rst @@ -5,15 +5,14 @@ GPU Support Overview -------- -From ``v1.12.0``, PyLops supports computations on GPUs powered by -`CuPy `_ (``cupy-cudaXX>=8.1.0``) and `cuSignal `_ (``cusignal>=0.16.0``). -They must be installed *before* PyLops is installed. +PyLops supports computations on GPUs powered by `CuPy `_ (``cupy-cudaXX>=10.6.0``). +This library must be installed *before* PyLops is installed. .. note:: - Set environment variables ``CUPY_PYLOPS=0`` and/or ``CUSIGNAL_PYLOPS=0`` to force PyLops to ignore - ``cupy`` and ``cusignal`` backends. - This can be also used if a previous version of ``cupy`` or ``cusignal`` is installed in your system, otherwise you will get an error when importing PyLops. + Set environment variable ``CUPY_PYLOPS=0`` to force PyLops to ignore the ``cupy`` backend. + This can be also used if a previous (or faulty) version of ``cupy`` is installed in your system, + otherwise you will get an error when importing PyLops. diff --git a/docs/source/installation.rst b/docs/source/installation.rst index 973b3840..a9c2d52a 100755 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -99,7 +99,8 @@ For a ``conda`` environment, run .. code-block:: bash - >> make dev-install_conda + >> make dev-install_conda # for x86 (Intel or AMD CPUs) + >> make dev-install_conda_arm # for arm (M-series Mac) This will create and activate an environment called ``pylops``, with all required and optional dependencies. @@ -524,16 +525,8 @@ disable this option. For more details of GPU-accelerated PyLops read :ref:`gpu`. CuPy ---- -`CuPy `_ is a library used as a drop-in replacement to NumPy -for GPU-accelerated -computations. Since many different versions of CuPy exist (based on the +`CuPy `_ is a library used as a drop-in replacement to NumPy and some parts of SciPy +for GPU-accelerated computations. Since many different versions of CuPy exist (based on the CUDA drivers of the GPU), users must install CuPy prior to installing PyLops. To do so, follow their -`installation instructions `__. - -cuSignal --------- -`cuSignal `_ is a library is used as a drop-in replacement to `SciPy Signal `_ for -GPU-accelerated computations. Similar to CuPy, users must install -cuSignal prior to installing PyLops. To do so, follow their -`installation instructions `__. +`installation instructions `__. \ No newline at end of file diff --git a/environment-dev-arm.yml b/environment-dev-arm.yml new file mode 100755 index 00000000..e04c5af7 --- /dev/null +++ b/environment-dev-arm.yml @@ -0,0 +1,37 @@ +name: pylops +channels: + - defaults + - conda-forge + - numba + - pytorch +dependencies: + - python>=3.6.4 + - pip + - numpy>=1.21.0 + - scipy>=1.4.0 + - pytorch>=1.2.0 + - pyfftw + - pywavelets + - sympy + - matplotlib + - ipython + - pytest + - Sphinx + - numpydoc + - numba + - pre-commit + - autopep8 + - isort + - black + - pip: + - devito + - scikit-fmm + - spgl1 + - pytest-runner + - setuptools_scm + - pydata-sphinx-theme + - sphinx-gallery + - nbsphinx + - image + - flake8 + - mypy 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/signalprocessing/convolve1d.py b/pylops/signalprocessing/convolve1d.py index 74057aa8..fd82eb91 100644 --- a/pylops/signalprocessing/convolve1d.py +++ b/pylops/signalprocessing/convolve1d.py @@ -100,7 +100,7 @@ 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 + self.h, self.method, self.dims, self.axis ) return self.convfunc(x, self.h, mode="same") @@ -109,7 +109,7 @@ 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 + self.hstar, self.method, self.dims, self.axis ) return self.convfunc(x, self.hstar, mode="same") @@ -165,7 +165,7 @@ 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 + self.h, self.method, self.dims, self.axis ) x = np.pad(x, self.pad) y = self.convfunc(self.h, x, mode="same") @@ -177,7 +177,7 @@ def _rmatvec(self, x: NDArray) -> NDArray: 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 + self.hstar, self.method, self.dims, self.axis ) x = np.pad(x, self.padd) y = self.convfunc(self.hstar, x) diff --git a/pylops/signalprocessing/shift.py b/pylops/signalprocessing/shift.py index 43a9a500..b7d17ed7 100644 --- a/pylops/signalprocessing/shift.py +++ b/pylops/signalprocessing/shift.py @@ -109,7 +109,7 @@ def Shift( shift = _value_or_sized_to_array(shift) if shift.size == 1: - shift = np.exp(-1j * 2 * np.pi * Fop.f * shift) + shift = np.exp(-1j * 2 * np.pi * Fop.f * shift).astype(Fop.cdtype) Sop = Diagonal(shift, dims=dimsdiag, axis=axis, dtype=Fop.cdtype) else: # add dimensions to shift to match dimensions of model and data @@ -120,7 +120,7 @@ def Shift( sdims = np.ones(shift.ndim + 1, dtype=int) sdims[:axis] = shift.shape[:axis] sdims[axis + 1 :] = shift.shape[axis:] - shift = np.exp(-1j * 2 * np.pi * f * shift.reshape(sdims)) + shift = np.exp(-1j * 2 * np.pi * f * shift.reshape(sdims)).astype(Fop.cdtype) Sop = Diagonal(shift, dtype=Fop.cdtype) Op = Fop.H * Sop * Fop Op.dims = Op.dimsd = Fop.dims diff --git a/pylops/signalprocessing/sliding1d.py b/pylops/signalprocessing/sliding1d.py index 5a9232f9..5cb8c213 100644 --- a/pylops/signalprocessing/sliding1d.py +++ b/pylops/signalprocessing/sliding1d.py @@ -157,7 +157,7 @@ def Sliding1D( # create tapers if tapertype is not None: - tap = taper(nwin, nover, tapertype=tapertype) + tap = taper(nwin, nover, tapertype=tapertype).astype(Op.dtype) tapin = tap.copy() tapin[:nover] = 1 tapend = tap.copy() @@ -172,7 +172,9 @@ def Sliding1D( if tapertype is None: OOp = BlockDiag([Op for _ in range(nwins)]) else: - OOp = BlockDiag([Diagonal(taps[itap].ravel()) * Op for itap in range(nwins)]) + OOp = BlockDiag( + [Diagonal(taps[itap].ravel(), dtype=Op.dtype) * Op for itap in range(nwins)] + ) combining = HStack( [ diff --git a/pylops/signalprocessing/sliding2d.py b/pylops/signalprocessing/sliding2d.py index 5af92b0b..c97802c6 100644 --- a/pylops/signalprocessing/sliding2d.py +++ b/pylops/signalprocessing/sliding2d.py @@ -191,7 +191,7 @@ def Sliding2D( # create tapers if tapertype is not None: - tap = taper2d(dimsd[1], nwin, nover, tapertype=tapertype) + tap = taper2d(dimsd[1], nwin, nover, tapertype=tapertype).astype(Op.dtype) tapin = tap.copy() tapin[:nover] = 1 tapend = tap.copy() @@ -206,7 +206,9 @@ def Sliding2D( if tapertype is None: OOp = BlockDiag([Op for _ in range(nwins)]) else: - OOp = BlockDiag([Diagonal(taps[itap].ravel()) * Op for itap in range(nwins)]) + OOp = BlockDiag( + [Diagonal(taps[itap].ravel(), dtype=Op.dtype) * Op for itap in range(nwins)] + ) combining = HStack( [ diff --git a/pylops/signalprocessing/sliding3d.py b/pylops/signalprocessing/sliding3d.py index f0a2b64a..7d21018c 100644 --- a/pylops/signalprocessing/sliding3d.py +++ b/pylops/signalprocessing/sliding3d.py @@ -183,13 +183,16 @@ def Sliding3D( # create tapers if tapertype is not None: - tap = taper3d(dimsd[2], nwin, nover, tapertype=tapertype) + tap = taper3d(dimsd[2], nwin, nover, tapertype=tapertype).astype(Op.dtype) # transform to apply if tapertype is None: OOp = BlockDiag([Op for _ in range(nwins)], nproc=nproc) else: - OOp = BlockDiag([Diagonal(tap.ravel()) * Op for _ in range(nwins)], nproc=nproc) + OOp = BlockDiag( + [Diagonal(tap.ravel(), dtype=Op.dtype) * Op for _ in range(nwins)], + nproc=nproc, + ) hstack = HStack( [ diff --git a/pylops/utils/backend.py b/pylops/utils/backend.py index ee2649f9..1cef8bdb 100644 --- a/pylops/utils/backend.py +++ b/pylops/utils/backend.py @@ -37,18 +37,13 @@ import cupyx.scipy.fft as cp_fft from cupyx.scipy.linalg import block_diag as cp_block_diag from cupyx.scipy.linalg import toeplitz as cp_toeplitz + from cupyx.scipy.signal import convolve as cp_convolve + from cupyx.scipy.signal import correlate as cp_correlate + from cupyx.scipy.signal import fftconvolve as cp_fftconvolve + from cupyx.scipy.signal import oaconvolve as cp_oaconvolve from cupyx.scipy.sparse import csc_matrix as cp_csc_matrix from cupyx.scipy.sparse import eye as cp_eye -if deps.cusignal_enabled: - import cusignal - -cu_message = "cupy package not installed. Use numpy arrays of " "install cupy." - -cusignal_message = ( - "cusignal package not installed. Use numpy arrays of" "install cusignal." -) - def get_module(backend: str = "numpy") -> ModuleType: """Returns correct numerical module based on backend string @@ -138,10 +133,7 @@ def get_convolve(x: npt.ArrayLike) -> Callable: if cp.get_array_module(x) == np: return convolve else: - if deps.cusignal_enabled: - return cusignal.convolution.convolve - else: - raise ModuleNotFoundError(cusignal_message) + return cp_convolve def get_fftconvolve(x: npt.ArrayLike) -> Callable: @@ -164,10 +156,7 @@ def get_fftconvolve(x: npt.ArrayLike) -> Callable: if cp.get_array_module(x) == np: return fftconvolve else: - if deps.cusignal_enabled: - return cusignal.convolution.fftconvolve - else: - raise ModuleNotFoundError(cusignal_message) + return cp_fftconvolve def get_oaconvolve(x: npt.ArrayLike) -> Callable: @@ -190,11 +179,7 @@ def get_oaconvolve(x: npt.ArrayLike) -> Callable: if cp.get_array_module(x) == np: return oaconvolve else: - raise NotImplementedError( - "oaconvolve not implemented in " - "cupy/cusignal. Consider using a different" - "option..." - ) + return cp_oaconvolve def get_correlate(x: npt.ArrayLike) -> Callable: @@ -217,10 +202,7 @@ def get_correlate(x: npt.ArrayLike) -> Callable: if cp.get_array_module(x) == np: return correlate else: - if deps.cusignal_enabled: - return cusignal.convolution.correlate - else: - raise ModuleNotFoundError(cusignal_message) + return cp_correlate def get_add_at(x: npt.ArrayLike) -> Callable: diff --git a/pylops/utils/deps.py b/pylops/utils/deps.py index c805a9a3..4b2d21e7 100644 --- a/pylops/utils/deps.py +++ b/pylops/utils/deps.py @@ -1,6 +1,5 @@ __all__ = [ "cupy_enabled", - "cusignal_enabled", "devito_enabled", "dtcwt_enabled", "numba_enabled", @@ -13,31 +12,49 @@ ] import os -from importlib import util - -# check package availability -cupy_enabled = ( - util.find_spec("cupy") is not None and int(os.getenv("CUPY_PYLOPS", 1)) == 1 -) -cusignal_enabled = ( - util.find_spec("cusignal") is not None and int(os.getenv("CUSIGNAL_PYLOPS", 1)) == 1 -) -devito_enabled = util.find_spec("devito") is not None -dtcwt_enabled = util.find_spec("dtcwt") is not None -numba_enabled = util.find_spec("numba") is not None -pyfftw_enabled = util.find_spec("pyfftw") is not None -pywt_enabled = util.find_spec("pywt") is not None -skfmm_enabled = util.find_spec("skfmm") is not None -spgl1_enabled = util.find_spec("spgl1") is not None -sympy_enabled = util.find_spec("sympy") is not None -torch_enabled = util.find_spec("torch") is not None +from importlib import import_module, util +from typing import Optional # error message at import of available package -def devito_import(message): +def cupy_import(message: Optional[str] = None) -> str: + # detect if cupy is available and the user is expecting to be used + cupy_test = ( + util.find_spec("cupy") is not None and int(os.getenv("CUPY_PYLOPS", 1)) == 1 + ) + # if cupy should be importable + if cupy_test: + # try importing it + try: + import_module("cupy") # noqa: F401 + + # if successful set the message to None. + cupy_message = None + # if unable to import but the package is installed + except (ImportError, ModuleNotFoundError) as e: + cupy_message = ( + f"Failed to import cupy, Falling back to CPU (error: {e}). " + "Please ensure your CUDA environment is set up correctly " + "for more details visit 'https://docs.cupy.dev/en/stable/install.html'" + ) + print(UserWarning(cupy_message)) + # if cupy_test is False, it means not installed or environment variable set to 0 + else: + cupy_message = ( + "Cupy package not installed or os.getenv('CUPY_PYLOPS') == 0. " + f"In order to be able to use {message} " + "ensure 'os.getenv('CUPY_PYLOPS') == 1' and run " + "'pip install cupy'; " + "for more details visit 'https://docs.cupy.dev/en/stable/install.html'" + ) + + return cupy_message + + +def devito_import(message: Optional[str] = None) -> str: if devito_enabled: try: - import devito # noqa: F401 + import_module("devito") # noqa: F401 devito_message = None except Exception as e: @@ -51,7 +68,7 @@ def devito_import(message): return devito_message -def dtcwt_import(message): +def dtcwt_import(message: Optional[str] = None) -> str: if dtcwt_enabled: try: import dtcwt # noqa: F401 @@ -68,10 +85,10 @@ def dtcwt_import(message): return dtcwt_message -def numba_import(message): +def numba_import(message: Optional[str] = None) -> str: if numba_enabled: try: - import numba # noqa: F401 + import_module("numba") # noqa: F401 numba_message = None except Exception as e: @@ -87,10 +104,10 @@ def numba_import(message): return numba_message -def pyfftw_import(message): +def pyfftw_import(message: Optional[str] = None) -> str: if pyfftw_enabled: try: - import pyfftw # noqa: F401 + import_module("pyfftw") # noqa: F401 pyfftw_message = None except Exception as e: @@ -106,10 +123,10 @@ def pyfftw_import(message): return pyfftw_message -def pywt_import(message): +def pywt_import(message: Optional[str] = None) -> str: if pywt_enabled: try: - import pywt # noqa: F401 + import_module("pywt") # noqa: F401 pywt_message = None except Exception as e: @@ -125,10 +142,10 @@ def pywt_import(message): return pywt_message -def skfmm_import(message): +def skfmm_import(message: Optional[str] = None) -> str: if skfmm_enabled: try: - import skfmm # noqa: F401 + import_module("skfmm") # noqa: F401 skfmm_message = None except Exception as e: @@ -143,10 +160,10 @@ def skfmm_import(message): return skfmm_message -def spgl1_import(message): +def spgl1_import(message: Optional[str] = None) -> str: if spgl1_enabled: try: - import spgl1 # noqa: F401 + import_module("spgl1") # noqa: F401 spgl1_message = None except Exception as e: @@ -160,10 +177,10 @@ def spgl1_import(message): return spgl1_message -def sympy_import(message): +def sympy_import(message: Optional[str] = None) -> str: if sympy_enabled: try: - import sympy # noqa: F401 + import_module("sympy") # noqa: F401 sympy_message = None except Exception as e: @@ -175,3 +192,24 @@ def sympy_import(message): f'"pip install sympy".' ) return sympy_message + + +# Set package availability booleans +# cupy: the package is imported to check everything is working correctly, +# if not the package is disabled. We do this here as this library is used as drop-in +# replacement for many numpy and scipy routines when cupy arrays are provided. +# all other libraries: we simply check if the package is available and postpone its import +# to check everything is working correctly when a user tries to create an operator that requires +# such a package +cupy_enabled: bool = ( + True if (cupy_import() is None and int(os.getenv("CUPY_PYLOPS", 1)) == 1) else False +) +devito_enabled = util.find_spec("devito") is not None +dtcwt_enabled = util.find_spec("dtcwt") is not None +numba_enabled = util.find_spec("numba") is not None +pyfftw_enabled = util.find_spec("pyfftw") is not None +pywt_enabled = util.find_spec("pywt") is not None +skfmm_enabled = util.find_spec("skfmm") is not None +spgl1_enabled = util.find_spec("spgl1") is not None +sympy_enabled = util.find_spec("sympy") is not None +torch_enabled = util.find_spec("torch") is not None diff --git a/pylops/utils/metrics.py b/pylops/utils/metrics.py index e1e7a55c..6393dec2 100644 --- a/pylops/utils/metrics.py +++ b/pylops/utils/metrics.py @@ -5,10 +5,13 @@ "psnr", ] +from typing import Optional + import numpy as np +import numpy.typing as npt -def mae(xref, xcmp): +def mae(xref: npt.ArrayLike, xcmp: npt.ArrayLike) -> float: """Mean Absolute Error (MAE) Compute Mean Absolute Error between two vectors @@ -30,7 +33,7 @@ def mae(xref, xcmp): return mae -def mse(xref, xcmp): +def mse(xref: npt.ArrayLike, xcmp: npt.ArrayLike) -> float: """Mean Square Error (MSE) Compute Mean Square Error between two vectors @@ -52,7 +55,7 @@ def mse(xref, xcmp): return mse -def snr(xref, xcmp): +def snr(xref: npt.ArrayLike, xcmp: npt.ArrayLike) -> float: """Signal to Noise Ratio (SNR) Compute Signal to Noise Ratio between two vectors @@ -75,7 +78,12 @@ def snr(xref, xcmp): return snr -def psnr(xref, xcmp, xmax=None, xmin=0.0): +def psnr( + xref: npt.ArrayLike, + xcmp: npt.ArrayLike, + xmax: Optional[float] = None, + xmin: Optional[float] = 0.0, +) -> float: """Peak Signal to Noise Ratio (PSNR) Compute Peak Signal to Noise Ratio between two vectors diff --git a/pylops/utils/tapers.py b/pylops/utils/tapers.py index ea109a00..52c95e8e 100644 --- a/pylops/utils/tapers.py +++ b/pylops/utils/tapers.py @@ -7,7 +7,7 @@ "tapernd", ] -from typing import Tuple, Union +from typing import Optional, Tuple, Union import numpy as np import numpy.typing as npt @@ -59,6 +59,7 @@ def cosinetaper( nmask: int, ntap: int, square: bool = False, + exponent: Optional[float] = None, ) -> npt.ArrayLike: r"""1D Cosine or Cosine square taper @@ -71,8 +72,10 @@ def cosinetaper( Number of samples of mask ntap : :obj:`int` Number of samples of hanning tapering at edges - square : :obj:`bool` - Cosine square taper (``True``)or Cosine taper (``False``) + square : :obj:`bool`, optional + Cosine square taper (``True``) or Cosine taper (``False``) + exponent : :obj:`float`, optional + Exponent to apply to Cosine taper. If provided, takes precedence over ``square`` Returns ------- @@ -81,7 +84,8 @@ def cosinetaper( """ ntap = 0 if ntap == 1 else ntap - exponent = 1 if not square else 2 + if exponent is None: + exponent = 1 if not square else 2 cos_win = ( 0.5 * ( @@ -123,7 +127,8 @@ def taper( ntap : :obj:`int` Number of samples of hanning tapering at edges tapertype : :obj:`str`, optional - Type of taper (``hanning``, ``cosine``, ``cosinesquare`` or ``None``) + Type of taper (``hanning``, ``cosine``, + ``cosinesquare``, ``cosinesqrt`` or ``None``) Returns ------- @@ -137,6 +142,8 @@ def taper( tpr_1d = cosinetaper(nmask, ntap, False) elif tapertype == "cosinesquare": tpr_1d = cosinetaper(nmask, ntap, True) + elif tapertype == "cosinesqrt": + tpr_1d = cosinetaper(nmask, ntap, False, 0.5) else: tpr_1d = np.ones(nmask) return tpr_1d @@ -214,7 +221,7 @@ def taper3d( Number of samples of tapering at edges of first and second dimensions tapertype : :obj:`int` Type of taper (``hanning``, ``cosine``, - ``cosinesquare`` or ``None``) + ``cosinesquare``, ``cosinesqrt`` or ``None``) Returns ------- @@ -236,6 +243,9 @@ def taper3d( elif tapertype == "cosinesquare": tpr_y = cosinetaper(nmasky, ntapy, True) tpr_x = cosinetaper(nmaskx, ntapx, True) + elif tapertype == "cosinesqrt": + tpr_y = cosinetaper(nmasky, ntapy, False, 0.5) + tpr_x = cosinetaper(nmaskx, ntapx, False, 0.5) else: tpr_y = np.ones(nmasky) tpr_x = np.ones(nmaskx) @@ -266,7 +276,7 @@ def tapernd( Number of samples of tapering at edges of every dimension tapertype : :obj:`int` Type of taper (``hanning``, ``cosine``, - ``cosinesquare`` or ``None``) + ``cosinesquare``, ``cosinesqrt`` or ``None``) Returns ------- @@ -282,6 +292,8 @@ def tapernd( tpr = [cosinetaper(nm, nt, False) for nm, nt in zip(nmask, ntap)] elif tapertype == "cosinesquare": tpr = [cosinetaper(nm, nt, True) for nm, nt in zip(nmask, ntap)] + elif tapertype == "cosinesqrt": + tpr = [cosinetaper(nm, nt, False, 0.5) for nm, nt in zip(nmask, ntap)] else: tpr = [np.ones(nm) for nm in nmask] diff --git a/pylops/waveeqprocessing/blending.py b/pylops/waveeqprocessing/blending.py index ca18ccb3..adca6d93 100644 --- a/pylops/waveeqprocessing/blending.py +++ b/pylops/waveeqprocessing/blending.py @@ -76,7 +76,12 @@ def __init__( self.dt = dt self.times = times self.shiftall = shiftall - self.nttot = int(np.max(self.times) / self.dt + self.nt + 1) + if np.max(self.times) // dt == np.max(self.times) / dt: + # do not add extra sample as no shift will be applied + self.nttot = int(np.max(self.times) / self.dt + self.nt) + else: + # add 1 extra sample at the end + self.nttot = int(np.max(self.times) / self.dt + self.nt + 1) 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) @@ -111,7 +116,7 @@ def __init__( # 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) + diff = np.repeat(diff[:, np.newaxis], self.nr, axis=1).astype(self.dtype) self.ShiftOp = Shift( (self.ns, self.nr, self.nt + 1), diff, @@ -160,7 +165,11 @@ def _matvec_largerecs(self, x: NDArray) -> NDArray: if self.ShiftOps[i] is None: blended_data[:, shift_int : shift_int + self.nt] += x[i, :, :] else: - shifted_data = self.ShiftOps[i] * self.PadOp * x[i, :, :] + shifted_data = ( + self.ShiftOps[i] + .matvec(self.PadOp.matvec(x[i, :, :].ravel())) + .reshape(self.ShiftOps[i].dimsd) + ) blended_data[:, shift_int : shift_int + self.nt + 1] += shifted_data return blended_data @@ -172,11 +181,11 @@ def _rmatvec_largerecs(self, x: NDArray) -> NDArray: if self.ShiftOps[i] is None: deblended_data[i, :, :] = x[:, shift_int : shift_int + self.nt] else: - shifted_data = ( - self.PadOp.H - * self.ShiftOps[i].H - * x[:, shift_int : shift_int + self.nt + 1] - ) + shifted_data = self.PadOp.rmatvec( + self.ShiftOps[i].rmatvec( + x[:, shift_int : shift_int + self.nt + 1].ravel() + ) + ).reshape(self.PadOp.dims) deblended_data[i, :, :] = shifted_data return deblended_data diff --git a/pylops/waveeqprocessing/kirchhoff.py b/pylops/waveeqprocessing/kirchhoff.py index a743ec6d..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 @@ -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 @@ -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) @@ -1245,9 +991,6 @@ def _register_multiplications(self, engine: str) -> None: 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/pylops/waveeqprocessing/oneway.py b/pylops/waveeqprocessing/oneway.py index ee4d9474..b41779db 100644 --- a/pylops/waveeqprocessing/oneway.py +++ b/pylops/waveeqprocessing/oneway.py @@ -193,6 +193,7 @@ def Deghosting( dr: Sequence[float], vel: float, zrec: float, + kind: Optional[str] = "p", pd: Optional[NDArray] = None, win: Optional[NDArray] = None, npad: Union[Tuple[int], Tuple[int, int]] = (11, 11), @@ -206,13 +207,15 @@ def Deghosting( ) -> Tuple[NDArray, NDArray]: r"""Wavefield deghosting. - Apply seismic wavefield decomposition from single-component (pressure) - data. This process is also generally referred to as model-based deghosting. + Apply seismic wavefield decomposition from single-component (pressure or + vertical velocity) data. This process is also generally referred to as + model-based deghosting. Parameters ---------- p : :obj:`np.ndarray` - Pressure data of of size :math:`\lbrack n_{r_x}\,(\times n_{r_y}) + Pressure (or vertical velocity) data of of size + :math:`\lbrack n_{r_x}\,(\times n_{r_y}) \times n_t \rbrack` (or :math:`\lbrack n_{r_{x,\text{sub}}}\, (\times n_{r_{y,\text{sub}}}) \times n_t \rbrack` in case a ``restriction`` operator is provided. Note that @@ -231,6 +234,8 @@ def Deghosting( Velocity along the receiver array (must be constant) zrec : :obj:`float` Depth of receiver array + kind : :obj:`str`, optional + Type of data (``p`` or ``vz``) pd : :obj:`np.ndarray`, optional Direct arrival to be subtracted from ``p`` win : :obj:`np.ndarray`, optional @@ -260,14 +265,19 @@ def Deghosting( Returns ------- pup : :obj:`np.ndarray` - Up-going wavefield + Up-going pressure (or particle velocity) wavefield pdown : :obj:`np.ndarray` - Down-going wavefield + Down-going (or particle velocity) wavefield + + Raises + ------ + ValueError + If ``kind`` is not "p" or "vz". Notes ----- - Up- and down-going components of seismic data :math:`p^-(x, t)` - and :math:`p^+(x, t)` can be estimated from single-component data + The up- and down-going components of a seismic data (:math:`p^-(x, t)` + and :math:`p^+(x, t)`) can be estimated from single-component data :math:`p(x, t)` using a ghost model. The basic idea [1]_ is that of using a one-way propagator in the f-k domain @@ -284,16 +294,22 @@ def Deghosting( In a matrix form we can thus write the total wavefield as: .. math:: - \mathbf{p} - \mathbf{p_d} = (\mathbf{I} + \Phi) \mathbf{p}^- + \mathbf{p} - \mathbf{p_d} = (\mathbf{I} \pm \Phi) \mathbf{p}^- where :math:`\Phi` is one-way propagator implemented via the - :class:`pylops.waveeqprocessing.PhaseShift` operator. + :class:`pylops.waveeqprocessing.PhaseShift` operator. Note that :math:`+` is + used for the pressure data, whilst :math:`-` is used for the vertical velocity + data. .. [1] Amundsen, L., 1993, Wavenumber-based filtering of marine point-source data: GEOPHYSICS, 58, 1335–1348. - """ + # Check kind + if kind not in ["p", "vz"]: + raise ValueError("kind must be p or vz") + + # Identify dimensions ndims = p.ndim if ndims == 2: dims = (nt, nr) @@ -328,7 +344,11 @@ def Deghosting( ) # Decomposition operator - Dupop = Identity(nt * nrs, dtype=p.dtype) + Pop + if kind == "p": + Dupop = Identity(nt * nrs, dtype=p.dtype) + Pop + else: + Dupop = Identity(nt * nrs, dtype=p.dtype) - Pop + if dottest: Dottest(Dupop, nt * nrs, nt * nrs, verb=True) diff --git a/pyproject.toml b/pyproject.toml index cd6759fa..9c338a48 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,7 +27,7 @@ classifiers = [ "Programming Language :: Python :: 3.8", "Programming Language :: Python :: 3.9", "Programming Language :: Python :: 3.10", - "Topic :: Scientific/Engineering : Mathematics", + "Topic :: Scientific/Engineering :: Mathematics", ] dependencies = [ "numpy >= 1.21.0", diff --git a/pytests/test_ffts.py b/pytests/test_ffts.py index 6d6c6c72..f912291c 100755 --- a/pytests/test_ffts.py +++ b/pytests/test_ffts.py @@ -162,7 +162,7 @@ def test_unknown_engine(par): (np.float16, 1), (np.float32, 4), (np.float64, 11), - (np.float128, 11), + (np.longdouble, 11), ], norm=["ortho", "none", "1/n"], ifftshift_before=[False, True], @@ -234,7 +234,7 @@ def test_FFT_small_real(par): (np.float16, 1), (np.float32, 3), (np.float64, 11), - (np.float128, 11), + (np.longdouble, 11), ], ifftshift_before=[False, True], engine=["numpy", "fftw", "scipy"], @@ -280,7 +280,7 @@ def test_FFT_random_real(par): par_lists_fft_small_cpx = dict( - dtype_precision=[(np.complex64, 4), (np.complex128, 11), (np.complex256, 11)], + dtype_precision=[(np.complex64, 4), (np.complex128, 11), (np.clongdouble, 11)], norm=["ortho", "none", "1/n"], ifftshift_before=[False, True], fftshift_after=[False, True], @@ -344,10 +344,10 @@ def test_FFT_small_complex(par): (np.float16, 1), (np.float32, 3), (np.float64, 11), - (np.float128, 11), + (np.longdouble, 11), (np.complex64, 3), (np.complex128, 11), - (np.complex256, 11), + (np.clongdouble, 11), ], ifftshift_before=[False, True], fftshift_after=[False, True], @@ -426,7 +426,7 @@ def test_FFT_random_complex(par): (np.float16, 1), (np.float32, 3), (np.float64, 11), - (np.float128, 11), + (np.longdouble, 11), ], ifftshift_before=[False, True], engine=["numpy", "scipy"], @@ -484,10 +484,10 @@ def test_FFT2D_random_real(par): (np.float16, 1), (np.float32, 3), (np.float64, 11), - (np.float128, 11), + (np.longdouble, 11), (np.complex64, 3), (np.complex128, 11), - (np.complex256, 11), + (np.clongdouble, 11), ], ifftshift_before=itertools.product([False, True], [False, True]), fftshift_after=itertools.product([False, True], [False, True]), @@ -566,7 +566,7 @@ def test_FFT2D_random_complex(par): (np.float16, 1), (np.float32, 3), (np.float64, 11), - (np.float128, 11), + (np.longdouble, 11), ], engine=["numpy", "scipy"], ) @@ -625,10 +625,10 @@ def test_FFTND_random_real(par): (np.float16, 1), (np.float32, 3), (np.float64, 11), - (np.float128, 11), + (np.longdouble, 11), (np.complex64, 3), (np.complex128, 11), - (np.complex256, 11), + (np.clongdouble, 11), ], engine=["numpy", "scipy"], ) @@ -700,7 +700,7 @@ def test_FFTND_random_complex(par): par_lists_fft2dnd_small_cpx = dict( - dtype_precision=[(np.complex64, 5), (np.complex128, 11), (np.complex256, 11)], + dtype_precision=[(np.complex64, 5), (np.complex128, 11), (np.clongdouble, 11)], norm=["ortho", "none", "1/n"], engine=["numpy", "scipy"], ) @@ -874,7 +874,15 @@ def test_FFT_1dsignal(par): assert_array_almost_equal(y_fftshift, np.fft.fftshift(y)) xadj = FFTop_fftshift.H * y_fftshift # adjoint is same as inverse for fft - xinv = lsqr(FFTop_fftshift, y_fftshift, damp=1e-10, iter_lim=10, atol=1e-8, btol=1e-8, show=0)[0] + xinv = lsqr( + FFTop_fftshift, + y_fftshift, + damp=1e-10, + iter_lim=10, + atol=1e-8, + btol=1e-8, + show=0, + )[0] assert_array_almost_equal(x[:imax], xadj[:imax], decimal=decimal) assert_array_almost_equal(x[:imax], xinv[:imax], decimal=decimal) @@ -958,7 +966,15 @@ def test_FFT_2dsignal(par): assert_array_almost_equal(D_fftshift, D2) dadj = FFTop_fftshift.H * D_fftshift # adjoint is same as inverse for fft - dinv = lsqr(FFTop_fftshift, D_fftshift, damp=1e-10, iter_lim=10, atol=1e-8, btol=1e-8, show=0)[0] + dinv = lsqr( + FFTop_fftshift, + D_fftshift, + damp=1e-10, + iter_lim=10, + atol=1e-8, + btol=1e-8, + show=0, + )[0] dadj = np.real(dadj.reshape(nt, nx)) dinv = np.real(dinv.reshape(nt, nx)) @@ -1016,7 +1032,15 @@ def test_FFT_2dsignal(par): assert_array_almost_equal(D_fftshift, D2) dadj = FFTop_fftshift.H * D_fftshift # adjoint is same as inverse for fft - dinv = lsqr(FFTop_fftshift, D_fftshift, damp=1e-10, iter_lim=10, atol=1e-8, btol=1e-8, show=0)[0] + dinv = lsqr( + FFTop_fftshift, + D_fftshift, + damp=1e-10, + iter_lim=10, + atol=1e-8, + btol=1e-8, + show=0, + )[0] dadj = np.real(dadj.reshape(nt, nx)) dinv = np.real(dinv.reshape(nt, nx)) @@ -1193,7 +1217,15 @@ def test_FFT_3dsignal(par): assert_array_almost_equal(D_fftshift, D2) dadj = FFTop_fftshift.H * D_fftshift # adjoint is same as inverse for fft - dinv = lsqr(FFTop_fftshift, D_fftshift, damp=1e-10, iter_lim=10, atol=1e-8, btol=1e-8, show=0)[0] + dinv = lsqr( + FFTop_fftshift, + D_fftshift, + damp=1e-10, + iter_lim=10, + atol=1e-8, + btol=1e-8, + show=0, + )[0] dadj = np.real(dadj.reshape(nt, nx, ny)) dinv = np.real(dinv.reshape(nt, nx, ny)) diff --git a/pytests/test_kirchhoff.py b/pytests/test_kirchhoff.py index 317821c3..d601d91b 100755 --- a/pytests/test_kirchhoff.py +++ b/pytests/test_kirchhoff.py @@ -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_oneway.py b/pytests/test_oneway.py index 2e8162d4..48f73a9e 100755 --- a/pytests/test_oneway.py +++ b/pytests/test_oneway.py @@ -22,8 +22,10 @@ "f0": 40, } -par1 = {"ny": 8, "nx": 10, "nt": 20, "dtype": "float32"} # even -par2 = {"ny": 9, "nx": 11, "nt": 21, "dtype": "float32"} # odd +par1 = {"ny": 8, "nx": 10, "nt": 20, "kind": "p", "dtype": "float32"} # even, p +par2 = {"ny": 9, "nx": 11, "nt": 21, "kind": "p", "dtype": "float32"} # odd, p +par1v = {"ny": 8, "nx": 10, "nt": 20, "kind": "vz", "dtype": "float32"} # even, vz +par2v = {"ny": 9, "nx": 11, "nt": 21, "kind": "vz", "dtype": "float32"} # odd, vz # deghosting params vel_sep = 1000.0 # velocity at separation level @@ -34,27 +36,31 @@ wav = ricker(t[:41], f0=parmod["f0"])[0] -@pytest.fixture(scope="module") +@pytest.fixture def create_data2D(): """Create 2d dataset""" - t0_plus = np.array([0.02, 0.08]) - t0_minus = t0_plus + 0.04 - vrms = np.array([1400.0, 1800.0]) - amp = np.array([1.0, -0.6]) - p2d_minus = hyperbolic2d(x, t, t0_minus, vrms, amp, wav)[1].T + def core(datakind): + t0_plus = np.array([0.02, 0.08]) + t0_minus = t0_plus + 0.04 + vrms = np.array([1400.0, 1800.0]) + amp = np.array([1.0, -0.6]) - kx = np.fft.ifftshift(np.fft.fftfreq(parmod["nx"], parmod["dx"])) - freq = np.fft.rfftfreq(parmod["nt"], parmod["dt"]) + p2d_minus = hyperbolic2d(x, t, t0_minus, vrms, amp, wav)[1].T - Pop = -PhaseShift(vel_sep, 2 * zrec, parmod["nt"], freq, kx) + kx = np.fft.ifftshift(np.fft.fftfreq(parmod["nx"], parmod["dx"])) + freq = np.fft.rfftfreq(parmod["nt"], parmod["dt"]) - # Decomposition operator - Dupop = Identity(parmod["nt"] * parmod["nx"]) + Pop + Pop = -PhaseShift(vel_sep, 2 * zrec, parmod["nt"], freq, kx) - p2d = Dupop * p2d_minus.ravel() - p2d = p2d.reshape(parmod["nt"], parmod["nx"]) - return p2d, p2d_minus + # Decomposition operator + Dupop = Identity(parmod["nt"] * parmod["nx"]) + datakind * Pop + + p2d = Dupop * p2d_minus.ravel() + p2d = p2d.reshape(parmod["nt"], parmod["nx"]) + return p2d, p2d_minus + + return core @pytest.mark.parametrize("par", [(par1), (par2)]) @@ -87,10 +93,10 @@ def test_PhaseShift_3dsignal(par): ) -@pytest.mark.parametrize("par", [(par1), (par2)]) +@pytest.mark.parametrize("par", [(par1), (par2), (par1v), (par2v)]) def test_Deghosting_2dsignal(par, create_data2D): """Deghosting of 2d data""" - p2d, p2d_minus = create_data2D + p2d, p2d_minus = create_data2D(1 if par["kind"] is "p" else -1) p2d_minus_inv, p2d_plus_inv = Deghosting( p2d, @@ -100,6 +106,7 @@ def test_Deghosting_2dsignal(par, create_data2D): parmod["dx"], vel_sep, zrec, + kind=par["kind"], win=np.ones_like(p2d), npad=0, ntaper=0, diff --git a/pytests/test_tapers.py b/pytests/test_tapers.py index 823097a4..320af932 100755 --- a/pytests/test_tapers.py +++ b/pytests/test_tapers.py @@ -40,9 +40,23 @@ "ntap": (4, 6), "tapertype": "cosinesquare", } # cosinesquare, even samples and taper +par7 = { + "nt": 21, + "nspat": (11, 13), + "ntap": (3, 5), + "tapertype": "cosinesqrt", +} # cosinesqrt, odd samples and taper +par8 = { + "nt": 20, + "nspat": (12, 16), + "ntap": (4, 6), + "tapertype": "cosinesqrt", +} # cosinesqrt, even samples and taper -@pytest.mark.parametrize("par", [(par1), (par2), (par3), (par4), (par5), (par6)]) +@pytest.mark.parametrize( + "par", [(par1), (par2), (par3), (par4), (par5), (par6), (par7), (par8)] +) def test_taper2d(par): """Create taper wavelet and check size and values""" tap = taper2d(par["nt"], par["nspat"][0], par["ntap"][0], par["tapertype"]) @@ -54,7 +68,9 @@ def test_taper2d(par): assert_array_equal(tap[par["nspat"][0] // 2], np.ones(par["nt"])) -@pytest.mark.parametrize("par", [(par1), (par2), (par3), (par4), (par5), (par6)]) +@pytest.mark.parametrize( + "par", [(par1), (par2), (par3), (par4), (par5), (par6), (par7), (par8)] +) def test_taper3d(par): """Create taper wavelet and check size and values""" tap = taper3d(par["nt"], par["nspat"], par["ntap"], par["tapertype"])