Merge pull request #8 from mrava87/master

Added prescaled option to MDC and Marchenko

Author: mrava87

pylops_distributed/waveeqprocessing/marchenko.py

@@ -38,10 +38,16 @@ class Marchenko():
     nsmooth : :obj:`int`, optional
         Number of samples of smoothing operator to apply to window
     saveRt : :obj:`bool`, optional
-        Save ``R`` and ``R^H`` to speed up the computation of adjoint of
+        Save ``R`` and ``R^H`` to speed up the computation of the adjoint of
         :class:`pylops_distributed.signalprocessing.Fredholm1` (``True``) or
         create ``R^H`` on-the-fly (``False``) Note that ``saveRt=True`` will be
         faster but double the amount of required memory
+    prescaled : :obj:`bool`, optional
+        Apply scaling to ``R`` (``False``) or not (``False``)
+        when performing spatial and temporal summations within the
+        :class:`pylops.waveeqprocessing.MDC` operator. In case
+        ``prescaled=True``, the ``R`` is assumed to have been pre-scaled by
+        the user.
     dtype : :obj:`bool`, optional
         Type of elements in input array.
 
@@ -73,7 +79,7 @@ class Marchenko():
 
     """
     def __init__(self, R, nt, dt=0.004, dr=1., wav=None, toff=0.0,
-                 nsmooth=10, saveRt=True, dtype='float64'):
+                 nsmooth=10, saveRt=True, prescaled=False, dtype='float64'):
         # Save inputs into class
         self.nt = nt
         self.dt = dt
@@ -82,6 +88,7 @@ def __init__(self, R, nt, dt=0.004, dr=1., wav=None, toff=0.0,
         self.toff = toff
         self.nsmooth = nsmooth
         self.saveRt = saveRt
+        self.prescaled = prescaled
         self.dtype = dtype
         self.explicit = False
 
@@ -163,9 +170,11 @@ def apply_onepoint(self, trav, dist=None, G0=None, nfft=None,
 
         # Create operators
         Rop = MDC(self.Rtwosided_fft, self.nt2, nv=1, dt=self.dt, dr=self.dr,
-                  twosided=True, conj=False, saveGt=self.saveRt)
+                  twosided=True, conj=False, saveGt=self.saveRt,
+                  prescaled=self.prescaled)
         R1op = MDC(self.Rtwosided_fft, self.nt2, nv=1, dt=self.dt, dr=self.dr,
-                   twosided=True, conj=True, saveGt=self.saveRt)
+                   twosided=True, conj=True, saveGt=self.saveRt,
+                   prescaled=self.prescaled)
         Rollop = Roll(self.nt2 * self.ns,
                       dims=(self.nt2, self.ns),
                       dir=0, shift=-1, dtype=self.dtype)
@@ -327,9 +336,11 @@ def apply_multiplepoints(self, trav, dist=None, G0=None, nfft=None,
 
         # Create operators
         Rop = MDC(self.Rtwosided_fft, self.nt2, nv=nvs, dt=self.dt,
-                  dr=self.dr, twosided=True, conj=False, saveGt=self.saveRt)
+                  dr=self.dr, twosided=True, conj=False, saveGt=self.saveRt,
+                  prescaled=self.prescaled)
         R1op = MDC(self.Rtwosided_fft, self.nt2, nv=nvs, dt=self.dt,
-                   dr=self.dr, twosided=True, conj=True, saveGt=self.saveRt)
+                   dr=self.dr, twosided=True, conj=True, saveGt=self.saveRt,
+                   prescaled=self.prescaled)
         Rollop = Roll(self.ns * nvs * self.nt2,
                       dims=(self.nt2, self.ns, nvs),
                       dir=0, shift=-1, dtype=self.dtype)

pylops_distributed/waveeqprocessing/mdd.py

@@ -8,8 +8,8 @@
 
 
 def MDC(G, nt, nv, dt=1., dr=1., twosided=True,
-        saveGt=True, conj=False, compute=(False, False),
-        todask=(False, False)):
+        saveGt=True, conj=False, prescaled=False,
+        compute=(False, False), todask=(False, False)):
     r"""Multi-dimensional convolution.
 
     Apply multi-dimensional convolution between two datasets.
@@ -43,6 +43,11 @@ def MDC(G, nt, nv, dt=1., dr=1., twosided=True,
         faster but double the amount of required memory
     conj : :obj:`str`, optional
         Perform Fredholm integral computation with complex conjugate of ``G``
+    prescaled : :obj:`bool`, optional
+        Apply scaling to kernel (``False``) or not (``False``) when performing
+        spatial and temporal summations. In case ``prescaled=True``, the
+        kernel is assumed to have been pre-scaled when passed to the MDC
+        routine.
     compute : :obj:`tuple`, optional
         Compute the outcome of forward and adjoint or simply define the graph
         and return a :obj:`dask.array`
@@ -57,7 +62,7 @@ def MDC(G, nt, nv, dt=1., dr=1., twosided=True,
 
     """
     return _MDC(G, nt, nv, dt=dt, dr=dr, twosided=twosided,
-                transpose=False, saveGt=saveGt, conj=conj,
+                transpose=False, saveGt=saveGt, conj=conj, prescaled=prescaled,
                 _Identity=Identity, _Transpose=Transpose,
                 _FFT=FFT, _Fredholm1=Fredholm1,
                 args_Fredholm1={'chunks': ((G.chunks[0], G.shape[2], nv),

pytests/test_marchenko.py

@@ -74,21 +74,30 @@
 # Load distributed Rs in frequency domain
 dRtwosided_fft = da.from_zarr(inputzarr)
 
+par1 = {'saveRt': True, 'prescaled':False}
+par2 = {'saveRt': False, 'prescaled':False}
+par3 = {'saveRt': True, 'prescaled':True}
+par4 = {'saveRt': False, 'prescaled':True}
 
-par1 = {'saveRt': True}  # square real
-par2 = {'saveRt': False}  # overdetermined real
 
-
-@pytest.mark.parametrize("par", [(par1), (par2)])
+@pytest.mark.parametrize("par", [(par1), (par2), (par3), (par4)])
 def test_Marchenko(par):
     """Dot-test and comparison with pylops for Marchenko.apply_onepoint
     """
-    dMarchenkoWM = dMarchenko(dRtwosided_fft, nt=nt, dt=dt, dr=dr,
+    if par['prescaled']:
+        dRtwosided_fft_sc = np.sqrt(2 * nt - 1) * dt * dr * dRtwosided_fft
+        Rtwosided_fft_sc = np.sqrt(2 * nt - 1) * dt * dr * Rtwosided_fft
+    else:
+        dRtwosided_fft_sc = dRtwosided_fft
+        Rtwosided_fft_sc = Rtwosided_fft
+
+    dMarchenkoWM = dMarchenko(dRtwosided_fft_sc, nt=nt, dt=dt, dr=dr,
                               wav=wav, toff=toff, nsmooth=nsmooth,
-                              saveRt=par['saveRt'])
+                              saveRt=par['saveRt'], prescaled=par['prescaled'])
 
-    MarchenkoWM = Marchenko(Rtwosided_fft, nt=nt, dt=dt, dr=dr,
-                            wav=wav, toff=toff, nsmooth=nsmooth)
+    MarchenkoWM = Marchenko(Rtwosided_fft_sc, nt=nt, dt=dt, dr=dr,
+                            wav=wav, toff=toff, nsmooth=nsmooth,
+                            prescaled=par['prescaled'])
 
     _, _, dp0_minus, dg_inv_minus, dg_inv_plus = \
         dMarchenkoWM.apply_onepoint(trav, nfft=2 ** 11, rtm=True, greens=True,
@@ -107,16 +116,24 @@ def test_Marchenko(par):
            np.linalg.norm(gsub_norm) < 1e-1
 
 
-@pytest.mark.parametrize("par", [(par1), (par2)])
-def test_Marchenko__multi(par):
+@pytest.mark.parametrize("par", [(par1), (par2), (par3), (par4)])
+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,
+    if par['prescaled']:
+        dRtwosided_fft_sc = np.sqrt(2 * nt - 1) * dt * dr * dRtwosided_fft
+        Rtwosided_fft_sc = np.sqrt(2 * nt - 1) * dt * dr * Rtwosided_fft
+    else:
+        dRtwosided_fft_sc = dRtwosided_fft
+        Rtwosided_fft_sc = Rtwosided_fft
+
+    dMarchenkoWM = dMarchenko(dRtwosided_fft_sc, nt=nt, dt=dt, dr=dr,
                               wav=wav, toff=toff, nsmooth=nsmooth,
-                              saveRt=par['saveRt'])
+                              saveRt=par['saveRt'], prescaled=par['prescaled'])
 
-    MarchenkoWM = Marchenko(Rtwosided_fft, nt=nt, dt=dt, dr=dr,
-                            wav=wav, toff=toff, nsmooth=nsmooth)
+    MarchenkoWM = Marchenko(Rtwosided_fft_sc, nt=nt, dt=dt, dr=dr,
+                            wav=wav, toff=toff, nsmooth=nsmooth,
+                            prescaled=par['prescaled'])
 
     _, _, dp0_minus, dg_inv_minus, dg_inv_plus = \
         dMarchenkoWM.apply_multiplepoints(trav_multi, nfft=2 ** 11, rtm=True,