From 32a40c20ecf450e11937b869a6a3ae3768d3a09e Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Sun, 1 Oct 2023 22:08:10 +0300
Subject: [PATCH 01/24] feature: restyling of Kirchhoff operator

---
 pylops/waveeqprocessing/kirchhoff.py | 607 +++++++--------------------
 pytests/test_kirchhoff.py            |  18 +-
 2 files changed, 149 insertions(+), 476 deletions(-)

diff --git a/pylops/waveeqprocessing/kirchhoff.py b/pylops/waveeqprocessing/kirchhoff.py
index a743ec6d..3ea6b6f4 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]_:
 
     .. 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}
+        \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 a filtered version of the wavelet :math:`w(t)`
+    as explained below (or the wavelet itself when ``wavfilter=False``).
+    Moreover, an angle scaling is included in the modelling operator,
+    where :math:`\theta=(\theta_s-\theta_r)/2` is half of the opening angle,
+    where :math:`\theta_s` and :math:`\theta_r` are the angles between the source-side
+    and receiver-side rays and the vertical at the image point.
 
-    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:
+    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{dist(\mathbf{x}, \mathbf{y})}}`
+    * ``3D``: :math:`a(\mathbf{x}, \mathbf{y})=\frac{1}{dist(\mathbf{x}, \mathbf{y})}`
 
-    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.
+    approximating the geometrical spreading of the wavefront. For ``mode=analytic``,
+    :math:`dist(\mathbf{x}, \mathbf{y})=\|\mathbf{x} - \mathbf{y}\|`, whilst for
+    ``mode=eikonal``, this is instead computed internally by the eikonal solver.
+
+    The wavelet filtering is applied as follows:
+
+    * ``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,15 +200,10 @@ 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
@@ -252,7 +239,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 +274,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 +295,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-4 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 +381,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,18 +627,19 @@ 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")
         elif dimsrc == 3 and dimrec == 3 and dimv == 3:
             # 3D
-            Wfilt = W * (-1j * 2 * np.pi * f)
+            Wfilt = W * (1j * 2 * np.pi * f)
         wavfilt = np.fft.irfft(Wfilt, n=len(wav))
         return wavfilt
 
@@ -750,225 +743,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 +755,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 +767,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 +825,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 +863,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 +875,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 +890,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 +936,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 +974,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 +987,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 +1009,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 +1021,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 +1034,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 +1058,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 +1070,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 +1083,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/pytests/test_kirchhoff.py b/pytests/test_kirchhoff.py
index 317821c3..b210a95b 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 reccommended 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,

From 7ed1d7f657ec37d602eb65d0df8d0fcaa7d51fbb Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Sun, 1 Oct 2023 22:10:17 +0300
Subject: [PATCH 02/24] minor: coverage action to use only python 3.8

---
 .github/workflows/codacy-coverage-reporter.yaml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/.github/workflows/codacy-coverage-reporter.yaml b/.github/workflows/codacy-coverage-reporter.yaml
index 5e3ed9fe..b16411b4 100644
--- a/.github/workflows/codacy-coverage-reporter.yaml
+++ b/.github/workflows/codacy-coverage-reporter.yaml
@@ -8,8 +8,8 @@ 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:

From ce165368c550e54db6102bca3085ca955f80e36c Mon Sep 17 00:00:00 2001
From: cako <c.dacostaf@gmail.com>
Date: Mon, 9 Oct 2023 15:59:41 -0700
Subject: [PATCH 03/24] update isort precommit hook

---
 .pre-commit-config.yaml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

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)

From e782fd3bfe56bd9294db4f23f1d4b23c3bc1b3f1 Mon Sep 17 00:00:00 2001
From: cako <c.dacostaf@gmail.com>
Date: Mon, 9 Oct 2023 15:59:56 -0700
Subject: [PATCH 04/24] improve language

---
 pylops/waveeqprocessing/kirchhoff.py | 22 +++++++++++-----------
 pytests/test_kirchhoff.py            |  2 +-
 2 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/pylops/waveeqprocessing/kirchhoff.py b/pylops/waveeqprocessing/kirchhoff.py
index 3ea6b6f4..5c7fef59 100644
--- a/pylops/waveeqprocessing/kirchhoff.py
+++ b/pylops/waveeqprocessing/kirchhoff.py
@@ -137,18 +137,18 @@ class Kirchhoff(LinearOperator):
 
     .. math::
         d(\mathbf{x_r}, \mathbf{x_s}, t) =
-        \widetilde{w}(t) * \int_V \frac{2 cos(\theta)} {v(\mathbf{x})}
+        \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 a filtered version of the wavelet :math:`w(t)`
-    as explained below (or the wavelet itself when ``wavfilter=False``).
+    :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 :math:`\theta=(\theta_s-\theta_r)/2` is half of the opening angle,
-    where :math:`\theta_s` and :math:`\theta_r` are the angles between the source-side
-    and receiver-side rays and the vertical at the image point.
+    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:
@@ -168,14 +168,14 @@ class Kirchhoff(LinearOperator):
 
     On the  other hand, when ``dynamic=True``, the amplitude scaling is defined as
 
-    * ``2D``: :math:`a(\mathbf{x}, \mathbf{y})=\frac{1}{\sqrt{dist(\mathbf{x}, \mathbf{y})}}`
-    * ``3D``: :math:`a(\mathbf{x}, \mathbf{y})=\frac{1}{dist(\mathbf{x}, \mathbf{y})}`
+    * ``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:`dist(\mathbf{x}, \mathbf{y})=\|\mathbf{x} - \mathbf{y}\|`, whilst for
-    ``mode=eikonal``, this is instead computed internally by the eikonal solver.
+    :math:`\text{dist}(\mathbf{x}, \mathbf{y})=\|\mathbf{x} - \mathbf{y}\|`, whilst for
+    ``mode=eikonal``, this is computed internally by the Eikonal solver.
 
-    The wavelet filtering is applied as follows:
+    The wavelet filtering is applied as follows [3]_:
 
     * ``2D``: :math:`\tilde{W}(f)=\sqrt{j\omega} \cdot W(f)`
     * ``3D``: :math:`\tilde{W}(f)=j\omega \cdot W(f)`
diff --git a/pytests/test_kirchhoff.py b/pytests/test_kirchhoff.py
index b210a95b..64139544 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 reccommended behaviour)"""
+    or trav_src and trav_rec (new recommended behaviour)"""
 
     # new behaviour
     Dop = Kirchhoff(

From 81466584955fb08d59086f84f8ab569592a7ba6a Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Tue, 10 Oct 2023 15:59:25 +0300
Subject: [PATCH 05/24] minor: fixed epsdist and sign of 3d wavelet

---
 pylops/waveeqprocessing/kirchhoff.py | 16 ++++++++++------
 pytests/test_kirchhoff.py            |  2 +-
 2 files changed, 11 insertions(+), 7 deletions(-)

diff --git a/pylops/waveeqprocessing/kirchhoff.py b/pylops/waveeqprocessing/kirchhoff.py
index 5c7fef59..7cd4a7da 100644
--- a/pylops/waveeqprocessing/kirchhoff.py
+++ b/pylops/waveeqprocessing/kirchhoff.py
@@ -133,7 +133,7 @@ class Kirchhoff(LinearOperator):
     -----
     The Kirchhoff demigration operator synthesizes seismic data given a
     propagation velocity model :math:`v(\mathbf{x})` and a
-    reflectivity model :math:`m(\mathbf{x})`. In forward mode [1]_, [2]_:
+    reflectivity model :math:`m(\mathbf{x})`. In forward mode [1]_, [2]_, [3]_:
 
     .. math::
         d(\mathbf{x_r}, \mathbf{x_s}, t) =
@@ -175,7 +175,7 @@ class Kirchhoff(LinearOperator):
     :math:`\text{dist}(\mathbf{x}, \mathbf{y})=\|\mathbf{x} - \mathbf{y}\|`, whilst for
     ``mode=eikonal``, this is computed internally by the Eikonal solver.
 
-    The wavelet filtering is applied as follows [3]_:
+    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)`
@@ -209,14 +209,18 @@ class Kirchhoff(LinearOperator):
     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.
 
     """
@@ -302,7 +306,7 @@ def __init__(
                 ) = 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 1e-4 of max distance to avoid having
+                    # 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
@@ -639,7 +643,7 @@ def _wavelet_reshaping(
             raise NotImplementedError("2.D wavelet currently not available")
         elif dimsrc == 3 and dimrec == 3 and dimv == 3:
             # 3D
-            Wfilt = W * (1j * 2 * np.pi * f)
+            Wfilt = W * (-1j * 2 * np.pi * f)
         wavfilt = np.fft.irfft(Wfilt, n=len(wav))
         return wavfilt
 
diff --git a/pytests/test_kirchhoff.py b/pytests/test_kirchhoff.py
index 64139544..d601d91b 100755
--- a/pytests/test_kirchhoff.py
+++ b/pytests/test_kirchhoff.py
@@ -347,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(

From 00f7034def7a0d26f3e17d02ba234d9d43cbef52 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Sat, 14 Oct 2023 22:54:36 +0300
Subject: [PATCH 06/24] minor: fix doc to reflect wavelet sign

---
 pylops/waveeqprocessing/kirchhoff.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pylops/waveeqprocessing/kirchhoff.py b/pylops/waveeqprocessing/kirchhoff.py
index 7cd4a7da..06c29251 100644
--- a/pylops/waveeqprocessing/kirchhoff.py
+++ b/pylops/waveeqprocessing/kirchhoff.py
@@ -178,7 +178,7 @@ class Kirchhoff(LinearOperator):
     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)`
+    * ``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:

From 35aeae521d7a0d26aae12ad911572b032d006dc2 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Thu, 26 Oct 2023 21:56:56 +0300
Subject: [PATCH 07/24] doc: added pylops-mpi to extensions page

---
 docs/source/extensions.rst | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

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 <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.
\ No newline at end of file
+* `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!).

From 47b672d807ebb290d7da54162be65572a8a898ff Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Sun, 29 Oct 2023 11:37:20 +0300
Subject: [PATCH 08/24] minor: improved signal definition in sliding example

---
 examples/plot_sliding.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

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)
 

From c8f54fc844ae13ecb811565f73acd87fbde32248 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Fri, 10 Nov 2023 22:01:34 +0300
Subject: [PATCH 09/24] minor: prepare for v2.2.0

---
 CHANGELOG.md              | 15 +++++++++++++++
 docs/source/changelog.rst | 18 ++++++++++++++++++
 pyproject.toml            |  2 +-
 3 files changed, 34 insertions(+), 1 deletion(-)

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/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/pyproject.toml b/pyproject.toml
index 4deee5a6..d2f6c854 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",

From a6bca011bb279df962f2efb9ec36099de48163be Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Sat, 11 Nov 2023 21:29:51 +0300
Subject: [PATCH 10/24] build: added python3.11 to github actions

---
 .github/workflows/build.yaml                    | 2 +-
 .github/workflows/codacy-coverage-reporter.yaml | 4 ++++
 2 files changed, 5 insertions(+), 1 deletion(-)

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 b16411b4..2c81e290 100644
--- a/.github/workflows/codacy-coverage-reporter.yaml
+++ b/.github/workflows/codacy-coverage-reporter.yaml
@@ -14,6 +14,10 @@ jobs:
     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:

From 1c7c9e0e86283910425cfda2e36ad0340d91eee0 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Tue, 14 Nov 2023 21:22:22 +0300
Subject: [PATCH 11/24] fix: change BlendingContinuous to work with pylops
 solvers BlendingContinuous matvec/rmatvec used * and .H * for some internal
 calls to operators, which does not allow running this operator with pylops
 solvers as they disable ndarray_multiplication. Changing this to
 matvec/rmatvec calls followed by reshape enables the use of this operator
 with pylops solvers.

---
 pylops/waveeqprocessing/blending.py | 16 ++++++++++------
 1 file changed, 10 insertions(+), 6 deletions(-)

diff --git a/pylops/waveeqprocessing/blending.py b/pylops/waveeqprocessing/blending.py
index ca18ccb3..828c8ef5 100644
--- a/pylops/waveeqprocessing/blending.py
+++ b/pylops/waveeqprocessing/blending.py
@@ -160,7 +160,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 +176,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
 

From 9523ae13952ca050cde7df7cbf1dfc71876c66c8 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Wed, 22 Nov 2023 13:45:09 +0300
Subject: [PATCH 12/24] fix: ensure sliding ops work with fp32

---
 pylops/signalprocessing/sliding1d.py    |   6 +-
 pylops/signalprocessing/sliding1dNEW.py | 237 ++++++++++++++++++++++++
 pylops/signalprocessing/sliding2d.py    |   6 +-
 pylops/signalprocessing/sliding3d.py    |   7 +-
 4 files changed, 250 insertions(+), 6 deletions(-)
 create mode 100644 pylops/signalprocessing/sliding1dNEW.py

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/sliding1dNEW.py b/pylops/signalprocessing/sliding1dNEW.py
new file mode 100644
index 00000000..826e93a7
--- /dev/null
+++ b/pylops/signalprocessing/sliding1dNEW.py
@@ -0,0 +1,237 @@
+__all__ = [
+    "sliding1d_design",
+    "Sliding1D",
+]
+
+import logging
+from typing import Tuple, Union
+
+import numpy as np
+from numpy.lib.stride_tricks import sliding_window_view
+
+from pylops import LinearOperator
+from pylops.signalprocessing.sliding2d import _slidingsteps
+from pylops.utils._internal import _value_or_sized_to_tuple
+from pylops.utils.backend import (
+    get_array_module,
+    get_sliding_window_view,
+    to_cupy_conditional,
+)
+from pylops.utils.decorators import reshaped
+from pylops.utils.tapers import taper
+from pylops.utils.typing import InputDimsLike, NDArray
+
+logging.basicConfig(format="%(levelname)s: %(message)s", level=logging.WARNING)
+
+
+def sliding1d_design(
+    dimd: int,
+    nwin: int,
+    nover: int,
+    nop: int,
+) -> Tuple[int, int, Tuple[NDArray, NDArray], Tuple[NDArray, NDArray]]:
+    """Design Sliding1D operator
+
+    This routine can be used prior to creating the :class:`pylops.signalprocessing.Sliding1D`
+    operator to identify the correct number of windows to be used based on the dimension of the data (``dimsd``),
+    dimension of the window (``nwin``), overlap (``nover``),a and dimension of the operator acting in the model
+    space.
+
+    Parameters
+    ----------
+    dimsd : :obj:`tuple`
+        Shape of 2-dimensional data.
+    nwin : :obj:`tuple`
+        Number of samples of window.
+    nover : :obj:`tuple`
+        Number of samples of overlapping part of window.
+    nop : :obj:`tuple`
+        Size of model in the transformed domain.
+
+    Returns
+    -------
+    nwins : :obj:`int`
+        Number of windows.
+    dim : :obj:`int`
+        Shape of 2-dimensional model.
+    mwins_inends : :obj:`tuple`
+        Start and end indices for model patches.
+    dwins_inends : :obj:`tuple`
+        Start and end indices for data patches.
+
+    """
+    # data windows
+    dwin_ins, dwin_ends = _slidingsteps(dimd, nwin, nover)
+    dwins_inends = (dwin_ins, dwin_ends)
+    nwins = len(dwin_ins)
+
+    # model windows
+    dim = nwins * nop
+    mwin_ins, mwin_ends = _slidingsteps(dim, nop, 0)
+    mwins_inends = (mwin_ins, mwin_ends)
+
+    # print information about patching
+    logging.warning("%d windows required...", nwins)
+    logging.warning(
+        "data wins - start:%s, end:%s",
+        dwin_ins,
+        dwin_ends,
+    )
+    logging.warning(
+        "model wins - start:%s, end:%s",
+        mwin_ins,
+        mwin_ends,
+    )
+    return nwins, dim, mwins_inends, dwins_inends
+
+
+class Sliding1D(LinearOperator):
+    r"""1D Sliding transform operator.
+
+    Apply a transform operator ``Op`` repeatedly to slices of the model
+    vector in forward mode and slices of the data vector in adjoint mode.
+    More specifically, in forward mode the model vector is divided into
+    slices, each slice is transformed, and slices are then recombined in a
+    sliding window fashion.
+
+    This operator can be used to perform local, overlapping transforms (e.g.,
+    :obj:`pylops.signalprocessing.FFT`) on 1-dimensional arrays.
+
+    .. note:: The shape of the model has to be consistent with
+       the number of windows for this operator not to return an error. As the
+       number of windows depends directly on the choice of ``nwin`` and
+       ``nover``, it is recommended to first run ``sliding1d_design`` to obtain
+       the corresponding ``dims`` and number of windows.
+
+    .. warning:: Depending on the choice of `nwin` and `nover` as well as the
+       size of the data, sliding windows may not cover the entire data.
+       The start and end indices of each window will be displayed and returned
+       with running ``sliding1d_design``.
+
+    Parameters
+    ----------
+    Op : :obj:`pylops.LinearOperator`
+        Transform operator
+    dim : :obj:`tuple`
+        Shape of 1-dimensional model.
+    dimd : :obj:`tuple`
+        Shape of 1-dimensional data
+    nwin : :obj:`int`
+        Number of samples of window
+    nover : :obj:`int`
+        Number of samples of overlapping part of window
+    tapertype : :obj:`str`, optional
+        Type of taper (``hanning``, ``cosine``, ``cosinesquare`` or ``None``)
+    name : :obj:`str`, optional
+        .. versionadded:: 2.0.0
+
+        Name of operator (to be used by :func:`pylops.utils.describe.describe`)
+
+    Raises
+    ------
+    ValueError
+        Identified number of windows is not consistent with provided model
+        shape (``dims``).
+
+    """
+
+    def __init__(
+        self,
+        Op: LinearOperator,
+        dim: Union[int, InputDimsLike],
+        dimd: Union[int, InputDimsLike],
+        nwin: int,
+        nover: int,
+        tapertype: str = "hanning",
+        name: str = "S",
+    ) -> None:
+
+        dim: Tuple[int, ...] = _value_or_sized_to_tuple(dim)
+        dimd: Tuple[int, ...] = _value_or_sized_to_tuple(dimd)
+
+        # data windows
+        dwin_ins, dwin_ends = _slidingsteps(dimd[0], nwin, nover)
+        self.dwin_inends = (dwin_ins, dwin_ends)
+        nwins = len(dwin_ins)
+        self.nwin = nwin
+        self.nover = nover
+
+        # check windows
+        if nwins * Op.shape[1] != dim[0] and Op.shape[1] != dim[0]:
+            raise ValueError(
+                f"Model shape (dim={dim}) is not consistent with chosen "
+                f"number of windows. Run sliding1d_design to identify the "
+                f"correct number of windows for the current "
+                "model size..."
+            )
+
+        # create tapers
+        self.tapertype = tapertype
+        if self.tapertype is not None:
+            tap = taper(nwin, nover, tapertype=self.tapertype)
+            tapin = tap.copy()
+            tapin[:nover] = 1
+            tapend = tap.copy()
+            tapend[-nover:] = 1
+            self.taps = [
+                tapin,
+            ]
+            for i in range(1, nwins - 1):
+                self.taps.append(tap)
+            self.taps.append(tapend)
+            self.taps = np.vstack(self.taps)
+
+        # check if operator is applied to all windows simultaneously
+        self.simOp = False
+        if Op.shape[1] == dim[0]:
+            self.simOp = True
+        self.Op = Op
+
+        # create temporary shape and strides for cpy
+        self.shape_wins = None
+        self.strides_wins = None
+
+        super().__init__(
+            dtype=Op.dtype,
+            dims=(nwins, int(dim[0] // nwins)),
+            dimsd=dimd,
+            clinear=False,
+            name=name,
+        )
+
+    @reshaped
+    def _matvec(self, x: NDArray) -> NDArray:
+        ncp = get_array_module(x)
+        if self.tapertype is not None:
+            self.taps = to_cupy_conditional(x, self.taps)
+        y = ncp.zeros(self.dimsd, dtype=self.dtype)
+        if self.simOp:
+            x = self.Op @ x
+        for iwin0 in range(self.dims[0]):
+            if self.simOp:
+                xx = x[iwin0]
+            else:
+                xx = self.Op.matvec(x[iwin0])
+            if self.tapertype is not None:
+                xxwin = self.taps[iwin0] * xx
+            else:
+                xxwin = xx
+            y[self.dwin_inends[0][iwin0] : self.dwin_inends[1][iwin0]] += xxwin
+        return y
+
+    @reshaped
+    def _rmatvec(self, x: NDArray) -> NDArray:
+        ncp = get_array_module(x)
+        ncp_sliding_window_view = get_sliding_window_view(x)
+        if self.tapertype is not None:
+            self.taps = to_cupy_conditional(x, self.taps)
+        ywins = ncp_sliding_window_view(x, self.nwin)[:: self.nwin - self.nover]
+        if self.tapertype is not None:
+            ywins = ywins * self.taps
+        if self.simOp:
+            y = self.Op.H @ ywins
+        else:
+            y = ncp.zeros(self.dims, dtype=self.dtype)
+            for iwin0 in range(self.dims[0]):
+                y[iwin0] = self.Op.rmatvec(ywins[iwin0])
+        return y
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(
         [

From 01568dad194b9f9737f9e6d318975f9eaefebd07 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Wed, 22 Nov 2023 13:47:32 +0300
Subject: [PATCH 13/24] minor: removed sliding1dNEW

---
 pylops/signalprocessing/sliding1dNEW.py | 237 ------------------------
 1 file changed, 237 deletions(-)
 delete mode 100644 pylops/signalprocessing/sliding1dNEW.py

diff --git a/pylops/signalprocessing/sliding1dNEW.py b/pylops/signalprocessing/sliding1dNEW.py
deleted file mode 100644
index 826e93a7..00000000
--- a/pylops/signalprocessing/sliding1dNEW.py
+++ /dev/null
@@ -1,237 +0,0 @@
-__all__ = [
-    "sliding1d_design",
-    "Sliding1D",
-]
-
-import logging
-from typing import Tuple, Union
-
-import numpy as np
-from numpy.lib.stride_tricks import sliding_window_view
-
-from pylops import LinearOperator
-from pylops.signalprocessing.sliding2d import _slidingsteps
-from pylops.utils._internal import _value_or_sized_to_tuple
-from pylops.utils.backend import (
-    get_array_module,
-    get_sliding_window_view,
-    to_cupy_conditional,
-)
-from pylops.utils.decorators import reshaped
-from pylops.utils.tapers import taper
-from pylops.utils.typing import InputDimsLike, NDArray
-
-logging.basicConfig(format="%(levelname)s: %(message)s", level=logging.WARNING)
-
-
-def sliding1d_design(
-    dimd: int,
-    nwin: int,
-    nover: int,
-    nop: int,
-) -> Tuple[int, int, Tuple[NDArray, NDArray], Tuple[NDArray, NDArray]]:
-    """Design Sliding1D operator
-
-    This routine can be used prior to creating the :class:`pylops.signalprocessing.Sliding1D`
-    operator to identify the correct number of windows to be used based on the dimension of the data (``dimsd``),
-    dimension of the window (``nwin``), overlap (``nover``),a and dimension of the operator acting in the model
-    space.
-
-    Parameters
-    ----------
-    dimsd : :obj:`tuple`
-        Shape of 2-dimensional data.
-    nwin : :obj:`tuple`
-        Number of samples of window.
-    nover : :obj:`tuple`
-        Number of samples of overlapping part of window.
-    nop : :obj:`tuple`
-        Size of model in the transformed domain.
-
-    Returns
-    -------
-    nwins : :obj:`int`
-        Number of windows.
-    dim : :obj:`int`
-        Shape of 2-dimensional model.
-    mwins_inends : :obj:`tuple`
-        Start and end indices for model patches.
-    dwins_inends : :obj:`tuple`
-        Start and end indices for data patches.
-
-    """
-    # data windows
-    dwin_ins, dwin_ends = _slidingsteps(dimd, nwin, nover)
-    dwins_inends = (dwin_ins, dwin_ends)
-    nwins = len(dwin_ins)
-
-    # model windows
-    dim = nwins * nop
-    mwin_ins, mwin_ends = _slidingsteps(dim, nop, 0)
-    mwins_inends = (mwin_ins, mwin_ends)
-
-    # print information about patching
-    logging.warning("%d windows required...", nwins)
-    logging.warning(
-        "data wins - start:%s, end:%s",
-        dwin_ins,
-        dwin_ends,
-    )
-    logging.warning(
-        "model wins - start:%s, end:%s",
-        mwin_ins,
-        mwin_ends,
-    )
-    return nwins, dim, mwins_inends, dwins_inends
-
-
-class Sliding1D(LinearOperator):
-    r"""1D Sliding transform operator.
-
-    Apply a transform operator ``Op`` repeatedly to slices of the model
-    vector in forward mode and slices of the data vector in adjoint mode.
-    More specifically, in forward mode the model vector is divided into
-    slices, each slice is transformed, and slices are then recombined in a
-    sliding window fashion.
-
-    This operator can be used to perform local, overlapping transforms (e.g.,
-    :obj:`pylops.signalprocessing.FFT`) on 1-dimensional arrays.
-
-    .. note:: The shape of the model has to be consistent with
-       the number of windows for this operator not to return an error. As the
-       number of windows depends directly on the choice of ``nwin`` and
-       ``nover``, it is recommended to first run ``sliding1d_design`` to obtain
-       the corresponding ``dims`` and number of windows.
-
-    .. warning:: Depending on the choice of `nwin` and `nover` as well as the
-       size of the data, sliding windows may not cover the entire data.
-       The start and end indices of each window will be displayed and returned
-       with running ``sliding1d_design``.
-
-    Parameters
-    ----------
-    Op : :obj:`pylops.LinearOperator`
-        Transform operator
-    dim : :obj:`tuple`
-        Shape of 1-dimensional model.
-    dimd : :obj:`tuple`
-        Shape of 1-dimensional data
-    nwin : :obj:`int`
-        Number of samples of window
-    nover : :obj:`int`
-        Number of samples of overlapping part of window
-    tapertype : :obj:`str`, optional
-        Type of taper (``hanning``, ``cosine``, ``cosinesquare`` or ``None``)
-    name : :obj:`str`, optional
-        .. versionadded:: 2.0.0
-
-        Name of operator (to be used by :func:`pylops.utils.describe.describe`)
-
-    Raises
-    ------
-    ValueError
-        Identified number of windows is not consistent with provided model
-        shape (``dims``).
-
-    """
-
-    def __init__(
-        self,
-        Op: LinearOperator,
-        dim: Union[int, InputDimsLike],
-        dimd: Union[int, InputDimsLike],
-        nwin: int,
-        nover: int,
-        tapertype: str = "hanning",
-        name: str = "S",
-    ) -> None:
-
-        dim: Tuple[int, ...] = _value_or_sized_to_tuple(dim)
-        dimd: Tuple[int, ...] = _value_or_sized_to_tuple(dimd)
-
-        # data windows
-        dwin_ins, dwin_ends = _slidingsteps(dimd[0], nwin, nover)
-        self.dwin_inends = (dwin_ins, dwin_ends)
-        nwins = len(dwin_ins)
-        self.nwin = nwin
-        self.nover = nover
-
-        # check windows
-        if nwins * Op.shape[1] != dim[0] and Op.shape[1] != dim[0]:
-            raise ValueError(
-                f"Model shape (dim={dim}) is not consistent with chosen "
-                f"number of windows. Run sliding1d_design to identify the "
-                f"correct number of windows for the current "
-                "model size..."
-            )
-
-        # create tapers
-        self.tapertype = tapertype
-        if self.tapertype is not None:
-            tap = taper(nwin, nover, tapertype=self.tapertype)
-            tapin = tap.copy()
-            tapin[:nover] = 1
-            tapend = tap.copy()
-            tapend[-nover:] = 1
-            self.taps = [
-                tapin,
-            ]
-            for i in range(1, nwins - 1):
-                self.taps.append(tap)
-            self.taps.append(tapend)
-            self.taps = np.vstack(self.taps)
-
-        # check if operator is applied to all windows simultaneously
-        self.simOp = False
-        if Op.shape[1] == dim[0]:
-            self.simOp = True
-        self.Op = Op
-
-        # create temporary shape and strides for cpy
-        self.shape_wins = None
-        self.strides_wins = None
-
-        super().__init__(
-            dtype=Op.dtype,
-            dims=(nwins, int(dim[0] // nwins)),
-            dimsd=dimd,
-            clinear=False,
-            name=name,
-        )
-
-    @reshaped
-    def _matvec(self, x: NDArray) -> NDArray:
-        ncp = get_array_module(x)
-        if self.tapertype is not None:
-            self.taps = to_cupy_conditional(x, self.taps)
-        y = ncp.zeros(self.dimsd, dtype=self.dtype)
-        if self.simOp:
-            x = self.Op @ x
-        for iwin0 in range(self.dims[0]):
-            if self.simOp:
-                xx = x[iwin0]
-            else:
-                xx = self.Op.matvec(x[iwin0])
-            if self.tapertype is not None:
-                xxwin = self.taps[iwin0] * xx
-            else:
-                xxwin = xx
-            y[self.dwin_inends[0][iwin0] : self.dwin_inends[1][iwin0]] += xxwin
-        return y
-
-    @reshaped
-    def _rmatvec(self, x: NDArray) -> NDArray:
-        ncp = get_array_module(x)
-        ncp_sliding_window_view = get_sliding_window_view(x)
-        if self.tapertype is not None:
-            self.taps = to_cupy_conditional(x, self.taps)
-        ywins = ncp_sliding_window_view(x, self.nwin)[:: self.nwin - self.nover]
-        if self.tapertype is not None:
-            ywins = ywins * self.taps
-        if self.simOp:
-            y = self.Op.H @ ywins
-        else:
-            y = ncp.zeros(self.dims, dtype=self.dtype)
-            for iwin0 in range(self.dims[0]):
-                y[iwin0] = self.Op.rmatvec(ywins[iwin0])
-        return y

From 55c42ac053c3456d596cf400891076e8b1972d9a Mon Sep 17 00:00:00 2001
From: alex-rakowski <alexrakowski90@gmail.com>
Date: Sun, 26 Nov 2023 12:49:26 -0800
Subject: [PATCH 14/24] changing to 128/256 bit to c/longdouble (#552)

---
 pytests/test_ffts.py | 64 +++++++++++++++++++++++++++++++++-----------
 1 file changed, 48 insertions(+), 16 deletions(-)

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

From 294f9517587618d4ac779529e34fdde705062618 Mon Sep 17 00:00:00 2001
From: alex-rakowski <alexrakowski90@gmail.com>
Date: Tue, 28 Nov 2023 09:18:44 -0800
Subject: [PATCH 15/24] adding make dev-install_conda_arm (#553)

Adding dev-install_conda_arm in Makefile to create environment on M-series MacOs (without icc_rt)
---
 Makefile                     |  3 +++
 docs/source/installation.rst |  3 ++-
 environment-dev-arm.yml      | 37 ++++++++++++++++++++++++++++++++++++
 3 files changed, 42 insertions(+), 1 deletion(-)
 create mode 100755 environment-dev-arm.yml

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/docs/source/installation.rst b/docs/source/installation.rst
index 3b237f23..55cfec7b 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.
 
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

From a227115035465a9196bcc7f716ff1a7013c548cf Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Tue, 28 Nov 2023 21:16:09 +0300
Subject: [PATCH 16/24] feature: handles
 https://github.com/PyLops/pylops/issues/554

---
 pylops/waveeqprocessing/oneway.py | 42 ++++++++++++++++++++++--------
 pytests/test_oneway.py            | 43 ++++++++++++++++++-------------
 2 files changed, 56 insertions(+), 29 deletions(-)

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/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,

From f130a967d3a3dbd1a0a1417548199cb1271ef085 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Wed, 29 Nov 2023 14:12:34 +0300
Subject: [PATCH 17/24] fix: force dtype for shift operator inputs

---
 pylops/signalprocessing/shift.py    | 4 ++--
 pylops/waveeqprocessing/blending.py | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

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/waveeqprocessing/blending.py b/pylops/waveeqprocessing/blending.py
index 828c8ef5..db69db5e 100644
--- a/pylops/waveeqprocessing/blending.py
+++ b/pylops/waveeqprocessing/blending.py
@@ -111,7 +111,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,

From 32d60d5ba94eb8bb351984e6ba439b7044088de1 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Sat, 16 Dec 2023 11:46:45 +0300
Subject: [PATCH 18/24] doc: adding info about decorator in new operator doc

---
 docs/source/adding.rst | 45 +++++++++++++++++++++++++++++-------------
 1 file changed, 31 insertions(+), 14 deletions(-)

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

From 77eca0f6693396b7bc1813d543b62440beaab2a5 Mon Sep 17 00:00:00 2001
From: alex-rakowski <alexrakowski90@gmail.com>
Date: Sun, 17 Dec 2023 10:00:27 -0800
Subject: [PATCH 19/24] Changing how optional dependencies are import state are
 checked  (#551)

* adding initial check

* adding check_module_enabled func

* silly typo "module" instead of module

* switch all optionals to check_module_enabled

* changing to c/longdouble for tests

* cleaning up old code and adding noqa

* Revert "changing to c/longdouble for tests"

This reverts commit 06534b1154ad323b91ccd39e2331f45438614da6.

* numpydoc style docstring

* adding cupy/cusignal import tests and how enables are set

* adding message to be optional

* fixing logic in cupy/cusignal_enabled

* cleaning up old code snippets and typos

* small format changes

* more formatting changes

* changing import to import_module for pylint

* minor: restyling deps.py

* minor: added alex-rakowski to contributors

---------

Co-authored-by: mrava87 <matteoravasi@gmail.com>
---
 README.md               |   1 +
 docs/source/credits.rst |   3 +-
 pylops/utils/deps.py    | 135 +++++++++++++++++++++++++++++++---------
 3 files changed, 107 insertions(+), 32 deletions(-)

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/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 <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
\ No newline at end of file
+*  `Fedor Goncharov <https://github.com/fedor-goncharov>`_, fedor-goncharov
+*  `Alex Rakowski <https://github.com/alex-rakowski>`_, alex-rakowski
diff --git a/pylops/utils/deps.py b/pylops/utils/deps.py
index 7fad2838..3b84488c 100644
--- a/pylops/utils/deps.py
+++ b/pylops/utils/deps.py
@@ -12,30 +12,78 @@
 ]
 
 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
-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 cusignal_import(message: Optional[str] = None) -> str:
+    cusignal_test = (
+        util.find_spec("cusignal") is not None
+        and int(os.getenv("CUSIGNAL_PYLOPS", 1)) == 1
+    )
+    if cusignal_test:
+        try:
+            import_module("cusignal")  # noqa: F401
+
+            cusignal_message = None
+        except (ImportError, ModuleNotFoundError) as e:
+            cusignal_message = (
+                f"Failed to import cusignal. Falling back to CPU (error: {e}) . "
+                "Please ensure your CUDA environment is set up correctly; "
+                "for more details visit 'https://github.com/rapidsai/cusignal#installation'"
+            )
+            print(UserWarning(cusignal_message))
+    else:
+        cusignal_message = (
+            "Cusignal not installed or os.getenv('CUSIGNAL_PYLOPS') == 0. "
+            f"In order to be able to use {message} "
+            "ensure 'os.getenv('CUSIGNAL_PYLOPS') == 1' and run "
+            "'conda install cusignal'; "
+            "for more details visit ''https://github.com/rapidsai/cusignal#installation''"
+        )
+
+    return cusignal_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:
@@ -49,10 +97,10 @@ def devito_import(message):
     return devito_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:
@@ -68,10 +116,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:
@@ -87,10 +135,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:
@@ -106,10 +154,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:
@@ -124,10 +172,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:
@@ -141,10 +189,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:
@@ -156,3 +204,28 @@ def sympy_import(message):
             f'"pip install sympy".'
         )
     return sympy_message
+
+
+# Set package availability booleans
+# cupy and cusignal: the package is imported to check everything is working correctly,
+# if not the package is disabled. We do this here as both libraries are 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
+)
+cusignal_enabled: bool = (
+    True
+    if (cusignal_import() is None and int(os.getenv("CUSIGNAL_PYLOPS", 1)) == 1)
+    else False
+)
+devito_enabled = util.find_spec("devito") 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 b67845e3c06a562a187f72cd243c3b99ba89dd45 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Sun, 17 Dec 2023 21:13:23 +0300
Subject: [PATCH 20/24] minor: added typing to metrics

---
 pylops/utils/metrics.py | 16 ++++++++++++----
 1 file changed, 12 insertions(+), 4 deletions(-)

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

From 68ee8b926a06857dcb02ecd3c1f22b9127547777 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Tue, 9 Jan 2024 22:37:01 +0300
Subject: [PATCH 21/24] fix: better handling of nttot in BlendingContinuous

---
 pylops/waveeqprocessing/blending.py | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/pylops/waveeqprocessing/blending.py b/pylops/waveeqprocessing/blending.py
index db69db5e..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)

From 666803818c8d7ffa86f19b9c1b5318a4dbaa768a Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Tue, 9 Jan 2024 22:51:16 +0300
Subject: [PATCH 22/24] feature: remove cusignal dependency

---
 docs/source/faq.rst                   |  4 +--
 docs/source/gpu.rst                   | 11 ++++----
 docs/source/installation.rst          | 14 +++-------
 pylops/signalprocessing/convolve1d.py |  8 +++---
 pylops/utils/backend.py               | 36 +++++++------------------
 pylops/utils/deps.py                  | 39 ++-------------------------
 6 files changed, 26 insertions(+), 86 deletions(-)

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 <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.
 
 
 
diff --git a/docs/source/installation.rst b/docs/source/installation.rst
index 55cfec7b..563b62aa 100755
--- a/docs/source/installation.rst
+++ b/docs/source/installation.rst
@@ -514,16 +514,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>`__.
\ No newline at end of file
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/utils/backend.py b/pylops/utils/backend.py
index ee2649f9..b727e2ef 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,9 @@ 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)
+        # import cusignal
+        # return cusignal.convolution.convolve
+        return cp_convolve
 
 
 def get_fftconvolve(x: npt.ArrayLike) -> Callable:
@@ -164,10 +158,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 +181,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 +204,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 3b84488c..cbd5b3e3 100644
--- a/pylops/utils/deps.py
+++ b/pylops/utils/deps.py
@@ -1,6 +1,5 @@
 __all__ = [
     "cupy_enabled",
-    "cusignal_enabled",
     "devito_enabled",
     "numba_enabled",
     "pyfftw_enabled",
@@ -51,35 +50,6 @@ def cupy_import(message: Optional[str] = None) -> str:
     return cupy_message
 
 
-def cusignal_import(message: Optional[str] = None) -> str:
-    cusignal_test = (
-        util.find_spec("cusignal") is not None
-        and int(os.getenv("CUSIGNAL_PYLOPS", 1)) == 1
-    )
-    if cusignal_test:
-        try:
-            import_module("cusignal")  # noqa: F401
-
-            cusignal_message = None
-        except (ImportError, ModuleNotFoundError) as e:
-            cusignal_message = (
-                f"Failed to import cusignal. Falling back to CPU (error: {e}) . "
-                "Please ensure your CUDA environment is set up correctly; "
-                "for more details visit 'https://github.com/rapidsai/cusignal#installation'"
-            )
-            print(UserWarning(cusignal_message))
-    else:
-        cusignal_message = (
-            "Cusignal not installed or os.getenv('CUSIGNAL_PYLOPS') == 0. "
-            f"In order to be able to use {message} "
-            "ensure 'os.getenv('CUSIGNAL_PYLOPS') == 1' and run "
-            "'conda install cusignal'; "
-            "for more details visit ''https://github.com/rapidsai/cusignal#installation''"
-        )
-
-    return cusignal_message
-
-
 def devito_import(message: Optional[str] = None) -> str:
     if devito_enabled:
         try:
@@ -207,8 +177,8 @@ def sympy_import(message: Optional[str] = None) -> str:
 
 
 # Set package availability booleans
-# cupy and cusignal: the package is imported to check everything is working correctly,
-# if not the package is disabled. We do this here as both libraries are used as drop-in
+# 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
@@ -216,11 +186,6 @@ def sympy_import(message: Optional[str] = None) -> str:
 cupy_enabled: bool = (
     True if (cupy_import() is None and int(os.getenv("CUPY_PYLOPS", 1)) == 1) else False
 )
-cusignal_enabled: bool = (
-    True
-    if (cusignal_import() is None and int(os.getenv("CUSIGNAL_PYLOPS", 1)) == 1)
-    else False
-)
 devito_enabled = util.find_spec("devito") is not None
 numba_enabled = util.find_spec("numba") is not None
 pyfftw_enabled = util.find_spec("pyfftw") is not None

From 05df62b280c1fadb64587420c9c391637e8f2988 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Tue, 9 Jan 2024 22:55:22 +0300
Subject: [PATCH 23/24] minor: remove commented code from backend.py

---
 pylops/utils/backend.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/pylops/utils/backend.py b/pylops/utils/backend.py
index b727e2ef..1cef8bdb 100644
--- a/pylops/utils/backend.py
+++ b/pylops/utils/backend.py
@@ -133,8 +133,6 @@ def get_convolve(x: npt.ArrayLike) -> Callable:
     if cp.get_array_module(x) == np:
         return convolve
     else:
-        # import cusignal
-        # return cusignal.convolution.convolve
         return cp_convolve
 
 

From 38fe5eeae3d926cb503e78ee6a5cade909fb19f1 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Wed, 24 Jan 2024 22:02:06 +0300
Subject: [PATCH 24/24] feat: added exponent to cosinetaper

---
 pylops/utils/tapers.py | 26 +++++++++++++++++++-------
 pytests/test_tapers.py | 20 ++++++++++++++++++--
 2 files changed, 37 insertions(+), 9 deletions(-)

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/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"])