Merge pull request #336 from prisae/dottest-allclose

Check passed with np.isclose

Author: cako

MIGRATION_V1_V2.md

@@ -8,3 +8,5 @@ should be used as a checklist when converting a piece of code using PyLops from
 - XX
 - XX
 - XX
+
+- `utils.dottest`: The relative tolerance is new set via `rtol` (before `tol`), and absolute tolerance is new supported via the keyword `atol`. When calling it with purely positional arguments, note that after `rtol` comes now first `atol` before `complexflag`. When using `raiseerror=True` it now emits an `AttributeError` instead of a `ValueError`.

examples/plot_nmo.py

@@ -296,7 +296,7 @@ def _rmatvec(self, y):
 # and adjoint transforms truly are adjoints of each other.
 
 NMOOp = NMO(t, x, vel_t)
-dottest(NMOOp, *NMOOp.shape, tol=1e-4)
+dottest(NMOOp, *NMOOp.shape, rtol=1e-4)
 
 ###############################################################################
 # NMO using :py:class:`pylops.Spread`
@@ -370,7 +370,7 @@ def create_tables(taxis, haxis, vels_rms):
     dtable=nmo_dtable,  # Table of weights for linear interpolation
     engine="numba",  # numba or numpy
 ).H  # To perform NMO *correction*, we need the adjoint
-dottest(SpreadNMO, *SpreadNMO.shape, tol=1e-4)
+dottest(SpreadNMO, *SpreadNMO.shape, rtol=1e-4)
 
 ###############################################################################
 # We see it passes the dot test, but are the results right? Let's find out.

pylops/utils/dottest.py

@@ -7,7 +7,8 @@ def dottest(
     Op,
     nr=None,
     nc=None,
-    tol=1e-6,
+    rtol=1e-6,
+    atol=1e-21,
     complexflag=0,
     raiseerror=True,
     verb=False,
@@ -28,8 +29,11 @@ def dottest(
         Number of rows of operator (i.e., elements in data)
     nc : :obj:`int`
         Number of columns of operator (i.e., elements in model)
-    tol : :obj:`float`, optional
-        Dottest tolerance
+    rtol : :obj:`float`, optional
+        Relative dottest tolerance
+    atol : :obj:`float`, optional
+        Absolute dottest tolerance
+        .. versionadded:: 2.0.0
     complexflag : :obj:`bool`, optional
         Generate random vectors with
 
@@ -50,8 +54,8 @@ def dottest(
 
     Raises
     ------
-    ValueError
-        If dot-test is not verified within chosen tolerance.
+    AssertionError
+        If dot-test is not verified within chosen tolerances.
 
     Notes
     -----
@@ -74,25 +78,19 @@ def dottest(
     if nc is None:
         nc = Op.shape[1]
 
-    assert (nr, nc) == Op.shape, "Provided nr and nc do not match operator shape"
+    if (nr, nc) != Op.shape:
+        raise AssertionError("Provided nr and nc do not match operator shape")
 
     # make u and v vectors
-    if complexflag != 0:
-        rdtype = np.real(np.ones(1, Op.dtype)).dtype
+    rdtype = np.ones(1, Op.dtype).real.dtype
 
-    if complexflag in (0, 2):
-        u = ncp.random.randn(nc).astype(Op.dtype)
-    else:
-        u = ncp.random.randn(nc).astype(rdtype) + 1j * ncp.random.randn(nc).astype(
-            rdtype
-        )
+    u = ncp.random.randn(nc).astype(rdtype)
+    if complexflag not in (0, 2):
+        u = u + 1j * ncp.random.randn(nc).astype(rdtype)
 
-    if complexflag in (0, 1):
-        v = ncp.random.randn(nr).astype(Op.dtype)
-    else:
-        v = ncp.random.randn(nr).astype(rdtype) + 1j * ncp.random.randn(nr).astype(
-            rdtype
-        )
+    v = ncp.random.randn(nr).astype(rdtype)
+    if complexflag not in (0, 1):
+        v = v + 1j * ncp.random.randn(nr).astype(rdtype)
 
     y = Op.matvec(u)  # Op * u
     x = Op.rmatvec(v)  # Op'* v
@@ -110,17 +108,16 @@ def dottest(
     # complex numbers in subsequent prints also when using cupy arrays.
     xx, yy = np.array([to_numpy(xx)])[0], np.array([to_numpy(yy)])[0]
 
-    def passes(xx, yy, tol=tol):
-        """True if xx and yy are  the same within tolerance."""
-        return abs((yy - xx) / ((yy + xx + 1e-15) / 2)) < tol
-
-    # evaluate if dot test is passed (both real and imag parts)
-    passed = passes(np.real(xx), np.real(yy)) and passes(np.imag(xx), np.imag(yy))
+    # evaluate if dot test passed
+    passed = np.isclose(xx, yy, rtol, atol)
 
-    if not passed and raiseerror:
-        raise ValueError(f"Dot test failed, v^H(Opu)={yy} - u^H(Op^Hv)={xx}")
-    elif verb:
+    # verbosity or error raising
+    if (not passed and raiseerror) or verb:
         passed_status = "passed" if passed else "failed"
-        print(f"Dot test {passed_status}, v^H(Opu)={yy} - u^H(Op^Hv)={xx}")
+        msg = f"Dot test {passed_status}, v^H(Opu)={yy} - u^H(Op^Hv)={xx}"
+        if not passed and raiseerror:
+            raise AssertionError(msg)
+        else:
+            print(msg)
 
     return passed