fix: fixed Seislet to work with correct dtype (and added tests)

Author: mrava87

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

pytests/test_seislet.py

@@ -15,23 +15,20 @@
     "dx": 10,
     "dt": 0.004,
     "level": None,
-    "dtype": "float32",
 }  # nx power of 2, max level
 par2 = {
     "nx": 16,
     "nt": 30,
     "dx": 10,
     "dt": 0.004,
     "level": 2,
-    "dtype": "float32",
 }  # nx power of 2, smaller level
 par3 = {
     "nx": 13,
     "nt": 30,
     "dx": 10,
     "dt": 0.004,
     "level": 2,
-    "dtype": "float32",
 }  # nx not power of 2, max level
 
 np.random.seed(10)
@@ -118,21 +115,30 @@ def _predict_reshape(
     int(os.environ.get("TEST_CUPY_PYLOPS", 0)) == 1, reason="Not CuPy enabled"
 )
 @pytest.mark.parametrize("par", [(par1), (par2), (par3)])
-def test_Seislet(par):
-    """Dot-test and forward-inverse for Seislet"""
-    slope = np.random.normal(0, 0.1, (par["nx"], par["nt"]))
+@pytest.mark.parametrize("dtype", [np.float32, np.float64])
+def test_Seislet(par, dtype):
+    """Dot-test and forward/adjoint/inverse for Seislet"""
+    slope = np.random.normal(0, 0.1, (par["nx"], par["nt"])).astype(dtype)
 
     for kind in ("haar", "linear"):
         Sop = Seislet(
             slope,
             sampling=(par["dx"], par["dt"]),
             level=par["level"],
             kind=kind,
-            dtype=par["dtype"],
+            dtype=dtype,
+        )
+        dottest(
+            Sop,
+            Sop.shape[0],
+            par["nx"] * par["nt"],
+            rtol=1e-3 if dtype == np.float32 else 1e-6,
         )
-        dottest(Sop, Sop.shape[0], par["nx"] * par["nt"])
 
-        x = np.random.normal(0, 0.1, par["nx"] * par["nt"])
+        x = np.random.normal(0, 0.1, par["nx"] * par["nt"]).astype(dtype)
         y = Sop * x
+        xadj = Sop.H * y
         xinv = Sop.inverse(y)
-        assert_array_almost_equal(x, xinv)
+        assert y.dtype == dtype
+        assert xadj.dtype == dtype
+        assert_array_almost_equal(x, xinv, decimal=3 if dtype == np.float32 else 6)