Fix: dtype consistency in operators

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/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/waveeqprocessing/marchenko.py

@@ -301,13 +301,14 @@ def __init__(
 
         # Add negative time to reflection data and convert to frequency
         if not np.iscomplexobj(R):
+            dtypec = (np.empty(0, dtype=dtype) + 1j * np.empty(0, dtype=dtype)).dtype
             Rtwosided = np.concatenate(
                 (self.ncp.zeros((self.ns, self.nr, self.nt - 1), dtype=R.dtype), R),
                 axis=-1,
             )
             Rtwosided_fft = np.fft.rfft(Rtwosided, self.nt2, axis=-1) / np.sqrt(
                 self.nt2
-            )
+            ).astype(dtypec)
             self.Rtwosided_fft = Rtwosided_fft[..., :nfmax]
         else:
             self.Rtwosided_fft = R
@@ -427,13 +428,12 @@ def apply_onepoint(
             shift=-1,
             dtype=self.dtype,
         )
-        Wop = Diagonal(w.T.ravel())
-        Iop = Identity(self.nr * self.nt2)
+        Wop = Diagonal(w.T.ravel(), dtype=self.dtype)
+        Iop = Identity(self.nr * self.nt2, dtype=self.dtype)
         Mop = Block(
-            [[Iop, -1 * Wop * Rop], [-1 * Wop * Rollop * R1op, Iop]]
-        ) * BlockDiag([Wop, Wop])
-        Gop = Block([[Iop, -1 * Rop], [-1 * Rollop * R1op, Iop]])
-
+            [[Iop, -1 * Wop * Rop], [-1 * Wop * Rollop * R1op, Iop]], dtype=self.dtype
+        ) * BlockDiag([Wop, Wop], dtype=self.dtype)
+        Gop = Block([[Iop, -1 * Rop], [-1 * Rollop * R1op, Iop]], dtype=self.dtype)
         if dottest:
             Dottest(
                 Gop,
@@ -443,7 +443,6 @@ def apply_onepoint(
                 verb=True,
                 backend=get_module_name(self.ncp),
             )
-        if dottest:
             Dottest(
                 Mop,
                 2 * self.ns * self.nt2,
@@ -457,10 +456,14 @@ def apply_onepoint(
         if G0 is None:
             if self.wav is not None and nfft is not None:
                 G0 = (
-                    directwave(
-                        self.wav, trav, self.nt, self.dt, nfft=nfft, derivative=True
+                    (
+                        directwave(
+                            self.wav, trav, self.nt, self.dt, nfft=nfft, derivative=True
+                        )
                     )
-                ).T
+                    .astype(self.dtype)
+                    .T
+                )
                 G0 = to_cupy_conditional(self.Rtwosided_fft, G0)
             else:
                 msg = "wav and/or nfft are not provided. Provide either G0 or wav and nfft..."
@@ -634,12 +637,12 @@ def apply_multiplepoints(
             shift=-1,
             dtype=self.dtype,
         )
-        Wop = Diagonal(w.transpose(2, 0, 1).ravel())
+        Wop = Diagonal(w.transpose(2, 0, 1).ravel(), dtype=self.dtype)
         Iop = Identity(self.nr * nvs * self.nt2)
         Mop = Block(
-            [[Iop, -1 * Wop * Rop], [-1 * Wop * Rollop * R1op, Iop]]
-        ) * BlockDiag([Wop, Wop])
-        Gop = Block([[Iop, -1 * Rop], [-1 * Rollop * R1op, Iop]])
+            [[Iop, -1 * Wop * Rop], [-1 * Wop * Rollop * R1op, Iop]], dtype=self.dtype
+        ) * BlockDiag([Wop, Wop], dtype=self.dtype)
+        Gop = Block([[Iop, -1 * Rop], [-1 * Rollop * R1op, Iop]], dtype=self.dtype)
 
         if dottest:
             Dottest(
@@ -650,7 +653,6 @@ def apply_multiplepoints(
                 verb=True,
                 backend=get_module_name(self.ncp),
             )
-        if dottest:
             Dottest(
                 Mop,
                 2 * self.ns * nvs * self.nt2,
@@ -666,15 +668,19 @@ def apply_multiplepoints(
                 G0 = np.zeros((self.nr, nvs, self.nt), dtype=self.dtype)
                 for ivs in range(nvs):
                     G0[:, ivs] = (
-                        directwave(
-                            self.wav,
-                            trav[:, ivs],
-                            self.nt,
-                            self.dt,
-                            nfft=nfft,
-                            derivative=True,
+                        (
+                            directwave(
+                                self.wav,
+                                trav[:, ivs],
+                                self.nt,
+                                self.dt,
+                                nfft=nfft,
+                                derivative=True,
+                            )
                         )
-                    ).T
+                        .astype(self.dtype)
+                        .T
+                    )
                 G0 = to_cupy_conditional(self.Rtwosided_fft, G0)
             else:
                 msg = "wav and/or nfft are not provided. Provide either G0 or wav and nfft..."

pytests/test_avo.py

@@ -141,24 +141,39 @@ def test_zoeppritz_and_approx_multipleangles():
 
 
 @pytest.mark.parametrize("par", [(par1), (par2), (par3), (par4)])
-def test_AVOLinearModelling(par):
+@pytest.mark.parametrize("dtype", [np.float32, np.float64])
+def test_AVOLinearModelling(par, dtype):
     """Dot-test and inversion for AVOLinearModelling"""
     AVOop = AVOLinearModelling(
-        np.asarray(theta),
-        vsvp=par["vsvp"] if isinstance(par["vsvp"], float) else np.asarray(par["vsvp"]),
+        np.asarray(theta).astype(dtype),
+        vsvp=par["vsvp"]
+        if isinstance(par["vsvp"], float)
+        else np.asarray(par["vsvp"]).astype(dtype),
         nt0=nt0,
         linearization=par["linearization"],
+        dtype=dtype,
+    )
+    assert dottest(
+        AVOop,
+        ntheta * nt0,
+        3 * nt0,
+        rtol=1e-4 if dtype == np.float32 else 1e-6,
+        backend=backend,
     )
-    assert dottest(AVOop, ntheta * nt0, 3 * nt0, backend=backend)
 
+    d = AVOop * np.asarray(m).astype(dtype)
+    madj = AVOop.H * d
     minv = lsqr(
         AVOop,
-        AVOop * np.asarray(m),
-        x0=np.zeros_like(m),
+        d,
+        x0=np.zeros_like(m).astype(dtype),
         damp=1e-20,
         niter=1000,
         atol=1e-8,
         btol=1e-8,
         show=0,
     )[0]
+
+    assert d.dtype == dtype
+    assert madj.dtype == dtype
     assert_array_almost_equal(m, minv, decimal=3)