Merge pull request #747 from mrava87/fix-bilinearcast

Fix: dtype consistency in operators

Author: mrava87

pylops/avo/avo.py

@@ -677,7 +677,9 @@ def __init__(
             dimsd += self.spatdims
         super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dimsd, name=name)
 
-        self.G = self.ncp.concatenate([gs.T[:, self.ncp.newaxis] for gs in Gs], axis=1)
+        self.G = self.ncp.concatenate(
+            [gs.astype(dtype).T[:, self.ncp.newaxis] for gs in Gs], axis=1
+        )
         # add dimensions to G to account for horizonal axes
         for _ in range(len(self.spatdims)):
             self.G = self.G[..., np.newaxis]

pylops/basicoperators/laplacian.py

@@ -1,5 +1,6 @@
 __all__ = ["Laplacian"]
 
+import numpy as np
 
 from pylops import LinearOperator
 from pylops.basicoperators import SecondDerivative
@@ -124,13 +125,17 @@ def _calc_l2op(
         kind: Tderivkind,
         dtype: DTypeLike,
     ):
+        weights = np.array(weights, dtype=dtype)
+        sampling = np.array(sampling, dtype=dtype)
         l2op = SecondDerivative(
             dims, axis=axes[0], sampling=sampling[0], edge=edge, kind=kind, dtype=dtype
         )
         dims = l2op.dims
         l2op *= weights[0]
         for ax, samp, weight in zip(axes[1:], sampling[1:], weights[1:], strict=True):
-            l2op += weight * SecondDerivative(
+            tmpop = SecondDerivative(
                 dims, axis=ax, sampling=samp, edge=edge, kind=kind, dtype=dtype
             )
+            tmpop *= weight
+            l2op += tmpop
         return l2op

pylops/basicoperators/smoothing1d.py

@@ -87,7 +87,7 @@ def __init__(
     ):
         if nsmooth % 2 == 0:
             nsmooth += 1
-        h = np.ones(nsmooth) / float(nsmooth)
+        h = np.ones(nsmooth, dtype=dtype) / float(nsmooth)
         offset = (nsmooth - 1) // 2
         super().__init__(
             dims=dims, h=h, axis=axis, offset=offset, dtype=dtype, name=name

pylops/basicoperators/smoothing2d.py

@@ -80,7 +80,9 @@ def __init__(
             nsmooth[0] += 1
         if nsmooth[1] % 2 == 0:
             nsmooth[1] += 1
-        h = np.ones((nsmooth[0], nsmooth[1])) / float(nsmooth[0] * nsmooth[1])
+        h = np.ones((nsmooth[0], nsmooth[1]), dtype=dtype) / float(
+            nsmooth[0] * nsmooth[1]
+        )
         offset = [(nsmooth[0] - 1) // 2, (nsmooth[1] - 1) // 2]
         super().__init__(
             dims=dims, h=h, offset=offset, axes=axes, dtype=dtype, name=name

pylops/basicoperators/smoothingnd.py

@@ -78,7 +78,7 @@ def __init__(
         for i in range(len(nsmooth)):
             if nsmooth[i] % 2 == 0:
                 nsmooth[i] += 1
-        h = np.ones(nsmooth) / float(np.prod(nsmooth))
+        h = np.ones(nsmooth, dtype=dtype) / float(np.prod(nsmooth))
         offset = [(nsm - 1) // 2 for nsm in nsmooth]
         super().__init__(
             dims=dims, h=h, offset=offset, axes=axes, dtype=dtype, name=name

pylops/signalprocessing/bilinear.py

@@ -132,6 +132,7 @@ def __init__(
         )
 
         ncp = get_array_module(iava)
+
         # check non-unique pairs (works only with numpy arrays)
         _ensure_iava_is_unique(
             iava=to_numpy(iava),
@@ -141,10 +142,10 @@ def __init__(
         # find indices and weights
         self.iava_t = ncp.floor(iava[0]).astype(int)
         self.iava_b = self.iava_t + 1
-        self.weights_tb = iava[0] - self.iava_t
+        self.weights_tb = (iava[0] - self.iava_t).astype(self.dtype)
         self.iava_l = ncp.floor(iava[1]).astype(int)
         self.iava_r = self.iava_l + 1
-        self.weights_lr = iava[1] - self.iava_l
+        self.weights_lr = (iava[1] - self.iava_l).astype(self.dtype)
 
         # expand dims to weights for nd-arrays
         if ndims > 2:

pylops/signalprocessing/interp.py

@@ -55,7 +55,7 @@ def _linearinterp(
     # find indices and weights
     iva_l = ncp.floor(iava).astype(int)
     iva_r = iva_l + 1
-    weights = iava - iva_l
+    weights = (iava - iva_l).astype(dtype)
 
     # create operators
     Op = Diagonal(1 - weights, dims=dimsd, axis=axis, dtype=dtype) * Restriction(
@@ -81,7 +81,7 @@ def _sincinterp(
     nreg = dims[axis]
     ireg = ncp.arange(nreg)
     sinc = ncp.tile(iava[:, np.newaxis], (1, nreg)) - ncp.tile(ireg, (len(iava), 1))
-    sinc = ncp.sinc(sinc)
+    sinc = ncp.sinc(sinc).astype(dtype)
 
     # identify additional dimensions and create MatrixMult operator
     otherdims = np.array(dims)

pylops/signalprocessing/nonstatconvolve2d.py

@@ -484,7 +484,7 @@ def __rmatvec(
         # the same loop (see https://numba.pydata.org/numba-doc/latest/user/parallel.html?highlight=njit).
         # Until atomic operations are provided we create a temporary filter array and store intermediate
         # results from each ix and reduce them at the end.
-        hstmp = np.zeros((xdims[0], *hs.shape))
+        hstmp = np.zeros((xdims[0], *hs.shape), dtype=x.dtype)
         for ix in prange(xdims[0]):
             for iz in range(xdims[1]):
                 # find extremes of model where to apply h (in case h is going out of model)

pylops/signalprocessing/seislet.py

@@ -451,7 +451,7 @@ def __init__(
         super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dimsd, name=name)
 
         pad = [(0, ndimpow2 - self.dims[0])] + [(0, 0)] * (len(self.dims) - 1)
-        self.pad = Pad(self.dims, pad)
+        self.pad = Pad(self.dims, pad, dtype=self.dtype)
         self.nx, self.nt = self.dimsd
 
         # define levels
@@ -473,7 +473,9 @@ def __init__(
     def _matvec(self, x: NDArray) -> NDArray:
         x = self.pad.matvec(x)
         x = np.reshape(x, self.dimsd)
-        y = np.zeros((np.sum(self.levels_size) + self.levels_size[-1], self.nt))
+        y = np.zeros(
+            (np.sum(self.levels_size) + self.levels_size[-1], self.nt), dtype=self.dtype
+        )
         for ilevel in range(self.level):
             odd = x[1::2]
             even = x[::2]
@@ -519,7 +521,7 @@ def _rmatvec(self, x: NDArray) -> NDArray:
                     backward=False,
                     adj=True,
                 )
-                y = np.zeros((2 * even.shape[0], self.nt))
+                y = np.zeros((2 * even.shape[0], self.nt), dtype=self.dtype)
                 y[1::2] = odd
                 y[::2] = even
             y = self.pad.rmatvec(y.ravel())
@@ -542,7 +544,7 @@ def inverse(self, x: NDArray) -> NDArray:
             odd = res + self.predict(
                 even, self.dt, self.dx, self.slopes, repeat=ilevel - 1, backward=False
             )
-            y = np.zeros((2 * even.shape[0], self.nt))
+            y = np.zeros((2 * even.shape[0], self.nt), dtype=self.dtype)
             y[1::2] = odd
             y[::2] = even
         y = self.pad.rmatvec(y.ravel())

pylops/signalprocessing/sliding1d.py

@@ -199,7 +199,7 @@ def __init__(
         self.tapertype = tapertype
         self.savetaper = savetaper
         if self.tapertype is not None:
-            tap = taper(nwin, nover, tapertype=self.tapertype)
+            tap = taper(nwin, nover, tapertype=self.tapertype).astype(Op.dtype)
             tapin = tap.copy()
             tapin[:nover] = 1
             tapend = tap.copy()
@@ -247,7 +247,7 @@ def _matvec_savetaper(self, x: NDArray) -> NDArray:
             self.taps = to_cupy_conditional(x, self.taps)
         y = ncp.zeros(self.dimsd, dtype=self.dtype)
         if self.simOp:
-            x = self.Op @ x
+            x = self.Op.matvec(x.ravel()).reshape(self.Op.dimsd)
             if self.tapertype is not None:
                 x = self.taps * x
         for iwin0 in range(self.dims[0]):
@@ -270,7 +270,7 @@ def _rmatvec_savetaper(self, x: NDArray) -> NDArray:
         if self.tapertype is not None:
             ywins = ywins * self.taps
         if self.simOp:
-            y = self.Op.H @ ywins
+            y = self.Op.rmatvec(ywins.ravel()).reshape(self.dims)
         else:
             y = ncp.zeros(self.dims, dtype=self.dtype)
             for iwin0 in range(self.dims[0]):
@@ -284,7 +284,7 @@ def _matvec_nosavetaper(self, x: NDArray) -> NDArray:
             self.taps = to_cupy_conditional(x, self.taps)
         y = ncp.zeros(self.dimsd, dtype=self.dtype)
         if self.simOp:
-            x = self.Op @ x
+            x = self.Op.matvec(x.ravel()).reshape(self.Op.dimsd)
         for iwin0 in range(self.dims[0]):
             if self.simOp:
                 xxwin = x[iwin0]
@@ -311,7 +311,7 @@ def _rmatvec_nosavetaper(self, x: NDArray) -> NDArray:
             if self.tapertype is not None:
                 for iwin0 in range(self.dims[0]):
                     ywins = self._apply_taper(ywins, iwin0)
-            y = self.Op.H @ ywins
+            y = self.Op.rmatvec(ywins.ravel())
         else:
             y = ncp.zeros(self.dims, dtype=self.dtype)
             for iwin0 in range(self.dims[0]):