Skip to content

Added prescaled option to MDC and Marchenko #8

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Apr 13, 2020
Merged

Added prescaled option to MDC and Marchenko #8

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 17 additions & 6 deletions pylops_distributed/waveeqprocessing/marchenko.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand Down Expand Up @@ -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
Expand All @@ -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

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
11 changes: 8 additions & 3 deletions pylops_distributed/waveeqprocessing/mdd.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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`
Expand All @@ -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),
Expand Down
45 changes: 31 additions & 14 deletions pytests/test_marchenko.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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,
Expand Down