Improved efficiency of Fredholm1 when saveGt=False

Author: mrava87

pylops_distributed/signalprocessing/Fredholm1.py

@@ -61,10 +61,13 @@ def __init__(self, G, nz=1, saveGt=True, compute=(False, False),
         self.nz = nz
         self.nsl, self.nx, self.ny = G.shape
         self.G = G
-        if saveGt:
-            self.GT = G.transpose((0, 2, 1)).conj()
         self.shape = (self.nsl * self.nx * self.nz,
                       self.nsl * self.ny * self.nz)
+        if saveGt:
+            self.GT = G.transpose((0, 2, 1)).conj().persist()
+        else:
+            # choose what to transpose if Gt is not saved
+            self.transposeG = True if self.G.size < self.shape[0] else False
         self.dtype = np.dtype(dtype)
         self.compute = compute
         self.chunks = chunks
@@ -90,5 +93,9 @@ def _rmatvec(self, x):
         if hasattr(self, 'GT'):
             y = da.matmul(self.GT, x)
         else:
-            y = da.matmul(self.G.transpose((0, 2, 1)).conj(), x)
+            if self.transposeG:
+                y = da.matmul(self.G.transpose((0, 2, 1)).conj(), x)
+            else:
+                y = da.matmul(x.transpose(0, 2, 1).conj(),
+                              self.G).transpose(0, 2, 1).conj()
         return y.ravel()

pylops_distributed/waveeqprocessing/marchenko.py

@@ -153,8 +153,6 @@ def apply_onepoint(self, trav, dist=None, G0=None, nfft=None,
         # Create window
         trav_off = trav - self.toff
         trav_off = np.round(trav_off / self.dt).astype(np.int)
-
-        # window
         w = np.zeros((self.nr, self.nt))
         for ir in range(self.nr):
             w[ir, :trav_off[ir]] = 1
@@ -392,7 +390,6 @@ def apply_multiplepoints(self, trav, dist=None, G0=None, nfft=None,
         f1_inv_tot = \
             f1_inv + da.concatenate((np.zeros((self.nt2, self.nr, nvs)),
                                      fd_plus))
-
         if greens:
             # Create Green's functions
             g_inv = Gop * f1_inv_tot.flatten()

pylops_distributed/waveeqprocessing/mdd.py

@@ -1,9 +1,8 @@
 import logging
 from pylops.waveeqprocessing.mdd import _MDC
 
-from pylops_distributed import Diagonal, Identity, Transpose
+from pylops_distributed import Identity, Transpose
 from pylops_distributed.signalprocessing import FFT, Fredholm1
-from pylops_distributed.utils import dottest as Dottest
 
 logging.basicConfig(format='%(levelname)s: %(message)s', level=logging.WARNING)
 

pytests/test_fredholm.py

@@ -13,8 +13,8 @@
         'saveGt': False, 'imag': 0, 'dtype': 'float32'}  # real, unsaved Gt
 par3 = {'nsl': 4, 'ny': 6, 'nx': 4, 'nz': 5,
         'saveGt':True, 'imag': 1j, 'dtype': 'complex64'}  # complex, saved Gt
-par4 = {'nsl': 4, 'ny': 6, 'nx': 4, 'nz': 5, 'saveGt': False,
-        'imag': 1j, 'dtype': 'complex64'}  # complex, unsaved Gt
+par4 = {'nsl': 4, 'ny': 6, 'nx': 4, 'nz': 5,
+        'saveGt': False, 'imag': 1j, 'dtype': 'complex64'}  # complex, unsaved Gt
 par5 = {'nsl': 4, 'ny': 6, 'nx': 4, 'nz': 1,
         'saveGt': True, 'imag': 0, 'dtype': 'float32'}  # real, saved Gt, nz=1
 par6 = {'nsl': 4, 'ny': 6, 'nx': 4, 'nz': 1,
@@ -54,3 +54,4 @@ def test_Fredholm1(par):
     dy = dFop * x.ravel()
     y = Fop * x.ravel().compute()
     assert_array_almost_equal(dy, y, decimal=5)
+    

pytests/test_marchenko.py

@@ -75,12 +75,17 @@
 dRtwosided_fft = da.from_zarr(inputzarr)
 
 
-def test_Marchenko():
+par1 = {'saveRt': True}  # square real
+par2 = {'saveRt': False}  # overdetermined real
+
+
+@pytest.mark.parametrize("par", [(par1), (par2)])
+def test_Marchenko(par):
     """Dot-test and comparison with pylops for Marchenko.apply_onepoint
     """
     dMarchenkoWM = dMarchenko(dRtwosided_fft, nt=nt, dt=dt, dr=dr,
-                              wav=wav, toff=toff, nsmooth=nsmooth)
-
+                              wav=wav, toff=toff, nsmooth=nsmooth,
+                              saveRt=par['saveRt'])
 
     MarchenkoWM = Marchenko(Rtwosided_fft, nt=nt, dt=dt, dr=dr,
                             wav=wav, toff=toff, nsmooth=nsmooth)
@@ -102,11 +107,13 @@ def test_Marchenko():
            np.linalg.norm(gsub_norm) < 1e-1
 
 
-def test_Marchenko__multi():
+@pytest.mark.parametrize("par", [(par1), (par2)])
+def test_Marchenko__multi(par):
     """Dot-test and comparison with pylops for Marchenko.apply_multiplepoints
     """
     dMarchenkoWM = dMarchenko(dRtwosided_fft, nt=nt, dt=dt, dr=dr,
-                              wav=wav, toff=toff, nsmooth=nsmooth)
+                              wav=wav, toff=toff, nsmooth=nsmooth,
+                              saveRt=par['saveRt'])
 
     MarchenkoWM = Marchenko(Rtwosided_fft, nt=nt, dt=dt, dr=dr,
                             wav=wav, toff=toff, nsmooth=nsmooth)