NonStationaryFilters1D and NonStationaryFilters2D operators (#478)

* feat: Added NonStationaryFilters1D and NonStationaryFilters2D
First implemenation of non-stationary filter estimation operators
(point 4 in #466).

* minor: added tests for nonstatfilters operators

* minor: added checks that filter locations are within dims

* minor: reintroduced parallel in rmatvec

docs/source/api/index.rst

@@ -91,6 +91,8 @@ Signal processing
     ConvolveND
     NonStationaryConvolve1D
     NonStationaryConvolve2D
+    NonStationaryFilters1D
+    NonStationaryFilters2D
     Interp
     Bilinear
     FFT

examples/plot_nonstatfilter.py

@@ -0,0 +1,135 @@
+"""
+Non-stationary Filter Estimation
+================================
+This example shows how to use the :py:class:`pylops.signalprocessing.NonStationaryFilters1D`
+and :py:class:`pylops.signalprocessing.NonStationaryFilters2D` operators to perform non-stationary
+filter estimation that when convolved with an input signal produces a desired output signal.
+
+"""
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.signal.windows import gaussian
+
+import pylops
+from pylops.utils.wavelets import ricker
+
+plt.close("all")
+
+
+###############################################################################
+# We will start by creating a zero signal of length `nt` and we will
+# place a comb of unitary spikes. We also create a non-stationary filter defined by
+# 5 equally spaced `Ricker wavelets <http://subsurfwiki.org/wiki/Ricker_wavelet>`_
+# with dominant frequencies of :math:`f = 10, 15, 20, 25` and :math:`30` Hz.
+nt = 601
+dt = 0.004
+t = np.arange(nt) * dt
+tw = ricker(t[:51], f0=5)[1]
+
+fs = [10, 15, 20, 25, 30]
+wavs = np.stack([ricker(t[:51], f0=f)[0] for f in fs])
+
+x = np.zeros(nt)
+x[64 : nt - 64 : 64] = 1.0
+
+Cop = pylops.signalprocessing.NonStationaryFilters1D(
+    x, hsize=101, ih=(101, 201, 301, 401, 501)
+)
+
+y = Cop @ wavs
+
+wavsinv = Cop.div(y, niter=20)
+wavsinv = wavsinv.reshape(wavs.shape)
+
+fig, axs = plt.subplots(1, len(fs), sharey=True, figsize=(14, 3))
+fig.suptitle("1D Non-stationary filter estimation")
+for i, ax in enumerate(axs):
+    ax.plot(wavs[i], "k", lw=4, label="True")
+    ax.plot(wavsinv[i], "r", lw=2, label="Estimate")
+    ax.set_xlabel("Time [sec]")
+axs[0].legend()
+plt.tight_layout()
+
+###############################################################################
+# Finally, we repeat the same exercise with a 2-dimensional non-stationary
+# filter
+nx, nz = 201, 101
+
+wav1a, _, _ = ricker(t[:9], f0=12)
+wav1b, _, _ = ricker(t[:9], f0=30)
+wav2a = gaussian(15, 2.0)
+wav2b = gaussian(15, 4.0)
+
+wav11 = np.outer(wav1a, wav2a[np.newaxis]).T
+wav12 = np.outer(wav1b, wav2a[np.newaxis]).T
+wav21 = np.outer(wav1b, wav2b[np.newaxis]).T
+wav22 = np.outer(wav1a, wav2b[np.newaxis]).T
+wavsize = wav11.shape
+
+hs = np.zeros((2, 2, *wavsize))
+hs[0, 0] = wav11
+hs[0, 1] = wav12
+hs[1, 0] = wav21
+hs[1, 1] = wav22
+
+x = np.zeros((nx, nz))
+x[:, 21] = 1.0
+x[:, 41] = -1.0
+
+Cop = pylops.signalprocessing.NonStationaryFilters2D(
+    inp=x, hshape=wavsize, ihx=(21, 41), ihz=(21, 41), engine="numba"
+)
+
+y = Cop @ hs
+
+hsinv = Cop.div(y.ravel(), niter=50)
+hsinv = hsinv.reshape(hs.shape)
+
+fig, axs = plt.subplots(2, 2, figsize=(10, 5))
+fig.suptitle("True filters")
+axs[0, 0].imshow(hs[0, 0], cmap="gray", vmin=-1, vmax=1)
+axs[0, 0].axis("tight")
+axs[0, 0].set_title(r"$H_{1,1}$")
+axs[0, 1].imshow(hs[0, 1], cmap="gray", vmin=-1, vmax=1)
+axs[0, 1].axis("tight")
+axs[0, 1].set_title(r"$H_{1,2}$")
+axs[1, 0].imshow(hs[1, 0], cmap="gray", vmin=-1, vmax=1)
+axs[1, 0].axis("tight")
+axs[1, 0].set_title(r"$H_{2,1}$")
+axs[1, 1].imshow(hs[1, 1], cmap="gray", vmin=-1, vmax=1)
+axs[1, 1].axis("tight")
+axs[1, 1].set_title(r"$H_{2,2}$")
+plt.tight_layout()
+
+fig, axs = plt.subplots(2, 2, figsize=(10, 5))
+fig.suptitle("Estimated filters")
+axs[0, 0].imshow(hsinv[0, 0], cmap="gray", vmin=-1, vmax=1)
+axs[0, 0].axis("tight")
+axs[0, 0].set_title(r"$H_{1,1}$")
+axs[0, 1].imshow(hsinv[0, 1], cmap="gray", vmin=-1, vmax=1)
+axs[0, 1].axis("tight")
+axs[0, 1].set_title(r"$H_{1,2}$")
+axs[1, 0].imshow(hsinv[1, 0], cmap="gray", vmin=-1, vmax=1)
+axs[1, 0].axis("tight")
+axs[1, 0].set_title(r"$H_{2,1}$")
+axs[1, 1].imshow(hsinv[1, 1], cmap="gray", vmin=-1, vmax=1)
+axs[1, 1].axis("tight")
+axs[1, 1].set_title(r"$H_{2,2}$")
+plt.tight_layout()
+
+
+fig, axs = plt.subplots(2, 2, figsize=(10, 5))
+fig.suptitle("Estimation error")
+axs[0, 0].imshow(hs[0, 0] - hsinv[0, 0], cmap="gray", vmin=-1, vmax=1)
+axs[0, 0].axis("tight")
+axs[0, 0].set_title(r"$H_{1,1}$")
+axs[0, 1].imshow(hs[0, 1] - hsinv[0, 1], cmap="gray", vmin=-1, vmax=1)
+axs[0, 1].axis("tight")
+axs[0, 1].set_title(r"$H_{1,2}$")
+axs[1, 0].imshow(hs[1, 0] - hsinv[1, 0], cmap="gray", vmin=-1, vmax=1)
+axs[1, 0].axis("tight")
+axs[1, 0].set_title(r"$H_{2,1}$")
+axs[1, 1].imshow(hs[1, 1] - hsinv[1, 1], cmap="gray", vmin=-1, vmax=1)
+axs[1, 1].axis("tight")
+axs[1, 1].set_title(r"$H_{2,2}$")
+plt.tight_layout()

pylops/signalprocessing/__init__.py

@@ -12,6 +12,8 @@
     ConvolveND                      ND convolution operator.
     NonStationaryConvolve1D         1D nonstationary convolution operator.
     NonStationaryConvolve2D         2D nonstationary convolution operator.
+    NonStationaryFilters1D          1D nonstationary filter estimation operator.
+    NonStationaryFilters2D          2D nonstationary filter estimation operator.
     Interp                          Interpolation operator.
     Bilinear                        Bilinear interpolation operator.
     FFT                             One dimensional Fast-Fourier Transform.
@@ -67,6 +69,8 @@
     "Convolve2D",
     "NonStationaryConvolve1D",
     "NonStationaryConvolve2D",
+    "NonStationaryFilters1D",
+    "NonStationaryFilters2D",
     "Interp",
     "Bilinear",
     "Radon2D",

pylops/signalprocessing/nonstatconvolve1d.py

@@ -1,4 +1,7 @@
-__all__ = ["NonStationaryConvolve1D"]
+__all__ = [
+    "NonStationaryConvolve1D",
+    "NonStationaryFilters1D",
+]
 
 from typing import Union
 
@@ -50,13 +53,15 @@ class NonStationaryConvolve1D(LinearOperator):
         If filters ``hs`` have even size
     ValueError
         If ``ih`` is not regularly sampled
+    ValueError
+        If ``ih`` is outside the bounds defined by ``dims[axis]``
 
     Notes
     -----
     The NonStationaryConvolve1D operator applies non-stationary
     one-dimensional convolution between the input signal :math:`d(t)`
     and a bank of compact filter kernels :math:`h(t; t_i)`. Assuming
-    an input signal composed of :math:`N=4` samples, and filters at locations
+    an input signal composed of :math:`N=5` samples, and filters at locations
     :math:`t_1` and :math:`t_3`, the forward operator can be represented as follows:
 
     .. math::
@@ -100,12 +105,16 @@ def __init__(
             raise ValueError(
                 "the indices of filters 'ih' are must be regularly sampled"
             )
+        dims = _value_or_sized_to_tuple(dims)
+        if min(ih) < 0 or max(ih) >= dims[axis]:
+            raise ValueError(
+                "the indices of filters 'ih' must be larger than 0 and smaller than `dims`"
+            )
         self.hs = hs
         self.hsize = hs.shape[1]
         self.oh, self.dh, self.nh, self.eh = ih[0], ih[1] - ih[0], len(ih), ih[-1]
         self.axis = axis
 
-        dims = _value_or_sized_to_tuple(dims)
         super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dims, name=name)
 
     @property
@@ -183,3 +192,173 @@ def todense(self):
         )
         H = H[:, int(self.hsize // 2) : -int(self.hsize // 2) - 1]
         return H
+
+
+class NonStationaryFilters1D(LinearOperator):
+    r"""1D non-stationary filter estimation operator.
+
+    Estimate a non-stationary one-dimensional filter by non-stationary convolution.
+    In forward mode, a varying compact filter on a coarser grid is on-the-fly linearly interpolated prior
+    to being convolved with a fixed input signal. In adjoint mode, the output signal is first weighted by the
+    fixed input signal and then spread across multiple filters (i.e., adjoint of linear interpolation).
+
+    Parameters
+    ----------
+    inp : :obj:`numpy.ndarray`
+        Fixed input signal of size :math:`n_x`.
+    hsize : :obj:`int`
+        Size of the 1d compact filters (filters must have odd number of samples and are assumed to be
+        centered in the middle of the filter support).
+    ih : :obj:`tuple`
+        Indices of the locations of the filters ``hs`` in the model (and data). Note
+        that the filters must be regularly sampled, i.e. :math:`dh=\text{diff}(ih)=\text{const.}`
+    dtype : :obj:`str`, optional
+        Type of elements in input array.
+    name : :obj:`str`, optional
+        Name of operator (to be used by :func:`pylops.utils.describe.describe`)
+
+    Attributes
+    ----------
+    shape : :obj:`tuple`
+        Operator shape
+    explicit : :obj:`bool`
+        Operator contains a matrix that can be solved explicitly (``True``) or
+        not (``False``)
+
+    Raises
+    ------
+    ValueError
+        If filters ``hsize`` is a even number
+    ValueError
+        If ``ih`` is not regularly sampled
+    ValueError
+        If ``ih`` is outside the bounds defined by ``dims[axis]``
+
+    See Also
+    --------
+    NonStationaryConvolve1D : 1D non-stationary convolution operator.
+    NonStationaryFilters2D : 2D non-stationary filter estimation operator.
+
+    Notes
+    -----
+    The NonStationaryFilters1D operates in a similar fashion to the
+    :class:`pylops.signalprocessing.NonStationaryConvolve1D` operator. In practical applications,
+    this operator shall be used when interested to estimate a 1-dimensional non-stationary filter
+    given an input and output signal.
+
+    In forward mode, this operator uses the same implementation of the
+    :class:`pylops.signalprocessing.NonStationaryConvolve1D`, with the main difference that
+    the role of the filters and the input signal is swapped. Nevertheless, to understand how
+    to implement adjoint, mathematically we arrange the forward operator in a slightly different way.
+    Assuming once again an input signal composed of :math:`N=5` samples, and filters at locations
+    :math:`t_1` and :math:`t_3`, the forward operator can be represented as follows:
+
+    .. math::
+        \mathbf{y} =
+        \begin{bmatrix}
+           \mathbf{X}_0 & \mathbf{X}_1 & \vdots & \mathbf{X}_4
+        \end{bmatrix} \mathbf{L}
+        \begin{bmatrix}
+            \mathbf{h}_1 \\ \mathbf{h}_3
+        \end{bmatrix}
+
+    where :math:`\mathbf{L}` is an operator that linearly interpolates the filters from the available locations to
+    the entire input grid -- i.e., :math:`[\hat{\mathbf{h}}_0 \quad \mathbf{h}_1 \quad \hat{\mathbf{h}}_2 \quad
+    \mathbf{h}_3 \quad \hat{\mathbf{h}}_4]^T = \mathbf{L} [ \mathbf{h}_1 \quad \mathbf{h}_3]`. Finally,
+    :math:`\mathbf{X}_i` is a diagonal matrix containing the value :math:`x_i` along the
+    main diagonal. Note that in practice the filter may be shorter than the input and output signals and
+    the :math:`x_i` values are placed only at the effective positions of the filter along the diagonal matrices
+    :math:`\mathbf{X}_i`.
+
+    In adjoint mode, the output signal is first weighted by the fixed input signal and then spread across
+    multiple filters (i.e., adjoint of linear interpolation) as follows
+
+    .. math::
+        \mathbf{h} =
+        \mathbf{L}^H
+        \begin{bmatrix}
+           \mathbf{X}_0 \\ \mathbf{X}_1 \\ \vdots \\ \mathbf{X}_4
+        \end{bmatrix}
+        \mathbf{y}
+
+    """
+
+    def __init__(
+        self,
+        inp: NDArray,
+        hsize: int,
+        ih: InputDimsLike,
+        dtype: DTypeLike = "float64",
+        name: str = "C",
+    ) -> None:
+        if hsize % 2 == 0:
+            raise ValueError("filters hs must have odd length")
+        if len(np.unique(np.diff(ih))) > 1:
+            raise ValueError(
+                "the indices of filters 'ih' are must be regularly sampled"
+            )
+        if min(ih) < 0 or max(ih) >= inp.size:
+            raise ValueError(
+                "the indices of filters 'ih' must be larger than 0 and smaller than `dims`"
+            )
+        self.inp = inp
+        self.hsize = hsize
+        self.oh, self.dh, self.nh, self.eh = ih[0], ih[1] - ih[0], len(ih), ih[-1]
+
+        super().__init__(
+            dtype=np.dtype(dtype), dims=(len(ih), hsize), dimsd=inp.shape, name=name
+        )
+
+    # use same interpolation method as inNonStationaryConvolve1D
+    _interpolate_h = staticmethod(NonStationaryConvolve1D._interpolate_h)
+
+    @staticmethod
+    def _interpolate_hadj(htmp, hs, hextremes, ix, oh, dh, nh):
+        """find closest filters and spread weighted psf"""
+        ih_closest = int(np.floor((ix - oh) / dh))
+        if ih_closest < 0:
+            hs[0, hextremes[0] : hextremes[1]] += htmp
+        elif ih_closest >= nh - 1:
+            hs[nh - 1, hextremes[0] : hextremes[1]] += htmp
+        else:
+            dh_closest = (ix - oh) / dh - ih_closest
+            hs[ih_closest, hextremes[0] : hextremes[1]] += (1 - dh_closest) * htmp
+            hs[ih_closest + 1, hextremes[0] : hextremes[1]] += dh_closest * htmp
+        return hs
+
+    @reshaped
+    def _matvec(self, x: NDArray) -> NDArray:
+        y = np.zeros(self.dimsd, dtype=self.dtype)
+        for ix in range(self.dimsd[0]):
+            h = self._interpolate_h(x, ix, self.oh, self.dh, self.nh)
+            xextremes = (
+                max(0, ix - self.hsize // 2),
+                min(ix + self.hsize // 2 + 1, self.dimsd[0]),
+            )
+            hextremes = (
+                max(0, -ix + self.hsize // 2),
+                min(self.hsize, self.hsize // 2 + (self.dimsd[0] - ix)),
+            )
+            y[..., xextremes[0] : xextremes[1]] += (
+                self.inp[..., ix : ix + 1] * h[hextremes[0] : hextremes[1]]
+            )
+        return y
+
+    @reshaped
+    def _rmatvec(self, x: NDArray) -> NDArray:
+        hs = np.zeros(self.dims, dtype=self.dtype)
+        for ix in range(self.dimsd[0]):
+            xextremes = (
+                max(0, ix - self.hsize // 2),
+                min(ix + self.hsize // 2 + 1, self.dimsd[0]),
+            )
+            hextremes = (
+                max(0, -ix + self.hsize // 2),
+                min(self.hsize, self.hsize // 2 + (self.dimsd[0] - ix)),
+            )
+
+            htmp = self.inp[ix] * x[..., xextremes[0] : xextremes[1]]
+            hs = self._interpolate_hadj(
+                htmp, hs, hextremes, ix, self.oh, self.dh, self.nh
+            )
+        return hs