minor: fix deps

.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:

.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:

.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)

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

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

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

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

docs/source/changelog.rst

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

docs/source/credits.rst

@@ -20,4 +20,5 @@ Contributors
 *  `Aniket Singh Rawat <https://github.com/dikwickley>`_, dikwickley
 *  `Rohan Babbar <https://github.com/rohanbabbar04>`_, rohanbabbar04
 *  `Wei Zhang <https://github.com/ZhangWeiGeo>`_, ZhangWeiGeo
-*  `Fedor Goncharov <https://github.com/fedor-goncharov>`_, fedor-goncharov
+*  `Fedor Goncharov <https://github.com/fedor-goncharov>`_, fedor-goncharov
+*  `Alex Rakowski <https://github.com/alex-rakowski>`_, alex-rakowski

docs/source/extensions.rst

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

docs/source/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``.

docs/source/gpu.rst

@@ -5,15 +5,14 @@ GPU Support
 
 Overview
 --------
-From ``v1.12.0``, PyLops supports computations on GPUs powered by
-`CuPy <https://cupy.dev/>`_ (``cupy-cudaXX>=8.1.0``) and `cuSignal <https://docs.rapids.ai/api/cusignal/stable/>`_ (``cusignal>=0.16.0``).
-They must be installed *before* PyLops is installed.
+PyLops supports computations on GPUs powered by `CuPy <https://cupy.dev/>`_ (``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.
 
 
 

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 <https://cupy.dev/>`_ 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 <https://cupy.dev/>`_ 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 <https://docs.cupy.dev/en/stable/install.html>`__.
-
-cuSignal
---------
-`cuSignal <https://docs.rapids.ai/api/cusignal/stable/>`_ is a library is used as a drop-in replacement to `SciPy Signal <https://docs.scipy.org/doc/scipy/reference/signal.html>`_ for
-GPU-accelerated computations. Similar to CuPy, users must install
-cuSignal prior to installing PyLops. To do so, follow their
-`installation instructions <https://github.com/rapidsai/cusignal#installation>`__.
+`installation instructions <https://docs.cupy.dev/en/stable/install.html>`__.

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

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)
 

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)

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

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(
         [