From 27f0849b68b6613bb3ddcb15be0d405d6ca093bb Mon Sep 17 00:00:00 2001 From: Anang Sahroni Date: Sat, 25 Jun 2022 13:26:31 +0700 Subject: [PATCH 01/12] add harmonics find azim static method decorated with numba jit --- rfpy/harmonics.py | 124 +++++++++++++++++++++++++++++++++ rfpy/scripts/rfpy_harmonics.py | 51 +++++++++++++- 2 files changed, 173 insertions(+), 2 deletions(-) diff --git a/rfpy/harmonics.py b/rfpy/harmonics.py index a50bf91..8d54949 100644 --- a/rfpy/harmonics.py +++ b/rfpy/harmonics.py @@ -28,6 +28,7 @@ import numpy as np from obspy.core import Stream, Trace import matplotlib.pyplot as plt +import numba class Harmonics(object): @@ -106,6 +107,129 @@ def __init__(self, radialRF, transvRF=None, azim=0, xmin=0., xmax=10.): self.xmin = xmin self.xmax = xmax + @staticmethod + @numba.jit(nopython=True) + def dcomp_find_azim_numba(xmin=None, xmax=None, + nbin=None, nz=None, delta=None, + baz_list=None, + dataR_list=None, + dataT_list=None, + + ): + """ + Method to decompose radial and transverse receiver function + streams into back-azimuth harmonics and determine the main + orientation ``azim``, obtained by minimizing the B1 component + between ``xmin`` and ``xmax`` (i.e., time or depth). + + Parameters + ---------- + xmin : float + Minimum x axis value over which to calculate ``azim`` + xmax : float + Maximum x axis value over which to calculate ``azim`` + + Attributes + ---------- + hstream : :class:`~obspy.core.Stream` + Stream containing the 5 harmonics, oriented in direction ``azim`` + azim : float + Direction (azimuth) along which the B1 component of the stream + is minimized (between ``xmin`` and ``xmax``) + var : :class:`~numpy.ndarray` + Variance of the 5 harmonics between ``xmin`` and ``xmax`` + + """ + + if not xmin: + xmin = xmin + if not xmax: + xmax = xmax + + print() + print('Decomposing receiver functions into baz harmonics') + + # Some integers + nbin = nbin + nz = nz + naz = 180 + daz = np.float(360/naz) + deg2rad = np.pi/180. + + # Define depth range over which to calculate azimuth + indmin = int(xmin/delta) + indmax = int(xmax/delta) + + # Initialize work arrays + C0 = np.zeros((nz, naz)) + C1 = np.zeros((nz, naz)) + C2 = np.zeros((nz, naz)) + C3 = np.zeros((nz, naz)) + C4 = np.zeros((nz, naz)) + + # Loop over each depth step + for iz in range(nz): + + # Build matrices OBS and H for each azimuth + for iaz in range(naz): + + # Initialize work arrays + OBS = np.zeros(2*nbin) + H = np.zeros((2*nbin, 5)) + + azim = iaz*daz + + # Radial component + for irow, dataR, baz in zip(range(len(dataR_list)),dataR_list,baz_list): + + baz = baz + OBS[irow] = dataR[iz] + H[irow, 0] = 1.0 + H[irow, 1] = np.cos(deg2rad*(baz-azim)) + H[irow, 2] = np.sin(deg2rad*(baz-azim)) + H[irow, 3] = np.cos(2.*deg2rad*(baz-azim)) + H[irow, 4] = np.sin(2.*deg2rad*(baz-azim)) + + shift = 90. + + # Transverse component + for irow, dataT, baz in zip(range(len(dataT_list)),dataT_list,baz_list): + + baz = baz + OBS[irow+nbin] = dataT[iz] + H[irow+nbin, 0] = 0.0 + H[irow+nbin, 1] = np.cos(deg2rad*(baz+shift-azim)) + H[irow+nbin, 2] = np.sin(deg2rad*(baz+shift-azim)) + H[irow+nbin, 3] = np.cos(2.*deg2rad*(baz+shift/2.0-azim)) + H[irow+nbin, 4] = np.sin(2.*deg2rad*(baz+shift/2.0-azim)) + + # Solve system of equations with truncated SVD + u, s, v = np.linalg.svd(H) + s[s < 0.001] = 0. + CC = np.linalg.solve(s.reshape(s.shape[0],1) * v, u.T.dot(OBS)[:5]) + + # Fill up arrays + C0[iz, iaz] = np.float(CC[0]) + C1[iz, iaz] = np.float(CC[1]) + C2[iz, iaz] = np.float(CC[2]) + C3[iz, iaz] = np.float(CC[3]) + C4[iz, iaz] = np.float(CC[4]) + + # Minimize variance of third component over specific depth range to + # find azim + C1var = np.zeros(naz) + for iaz in range(naz): + C1var[iaz] = np.sqrt(np.mean(np.square(C1[indmin:indmax, iaz]))) + indaz = np.argmin(C1var) + + C0var = np.sqrt(np.mean(np.square(C0[indmin:indmax, indaz]))) + C1var = np.sqrt(np.mean(np.square(C1[indmin:indmax, indaz]))) + C2var = np.sqrt(np.mean(np.square(C2[indmin:indmax, indaz]))) + C3var = np.sqrt(np.mean(np.square(C3[indmin:indmax, indaz]))) + C4var = np.sqrt(np.mean(np.square(C4[indmin:indmax, indaz]))) + + return C0,C1,C2,C3,C4,C0var,C1var,C2var,C3var,C4var,indaz,daz + def dcomp_find_azim(self, xmin=None, xmax=None): """ Method to decompose radial and transverse receiver function diff --git a/rfpy/scripts/rfpy_harmonics.py b/rfpy/scripts/rfpy_harmonics.py index 273cb4d..9a805ff 100644 --- a/rfpy/scripts/rfpy_harmonics.py +++ b/rfpy/scripts/rfpy_harmonics.py @@ -88,6 +88,14 @@ def get_harmonics_arguments(argv=None): "Default behaviour uses short key form (NET.STN) for the folder " + "names, regardless of the key type of the database." ) + parser.add_argument( + "-UN", "--use-numba", + action="store_true", + dest="unumba", + default=False, + help="Use numba jit on calculation [Default False]" + ) + # Event Selection Criteria TimeGroup = parser.add_argument_group( @@ -257,7 +265,7 @@ def get_harmonics_arguments(argv=None): default="png", help="Specify format of figure. Can be any one of the valid" + "matplotlib formats: 'png', 'jpg', 'eps', 'pdf'. [Default 'png']") - + args = parser.parse_args(argv) # Check inputs @@ -502,7 +510,46 @@ def main(): # Stack with or without dip if args.find_azim: - harmonics.dcomp_find_azim(xmin=args.trange[0], xmax=args.trange[1]) + import time + t0=time.time() + + if args.unumba: + from obspy.core import Trace + str_stats = harmonics.radialRF[0].stats + baz_list = np.array([trace0.stats.baz for trace0 in harmonics.radialRF]) + dataR_list = np.array([trace0.data for trace0 in harmonics.radialRF]) + dataT_list = np.array([trace0.data for trace0 in harmonics.transvRF]) + C0,C1,C2,C3,C4,C0var,C1var,C2var,C3var,C4var,indaz,daz = harmonics.dcomp_find_azim_numba( + xmin=args.trange[0], xmax=args.trange[1], + nbin=len(harmonics.radialRF), + nz=len(harmonics.radialRF[0].data), + delta=harmonics.radialRF[0].stats.delta, + baz_list=baz_list, + dataR_list=dataR_list, + dataT_list=dataT_list, + + ) + # Put back into traces + A = Trace(data=C0[:, indaz], header=str_stats) + B1 = Trace(data=C1[:, indaz], header=str_stats) + B2 = Trace(data=C2[:, indaz], header=str_stats) + C1 = Trace(data=C3[:, indaz], header=str_stats) + C2 = Trace(data=C4[:, indaz], header=str_stats) + + # Put all treaces into stream + hstream = Stream(traces=[A, B1, B2, C1, C2]) + azim = indaz*daz + var = [C0var, C1var, C2var, C3var, C4var] + + harmonics.hstream = hstream + harmonics.azim = azim + harmonics.var = var + else: + harmonics.dcomp_find_azim(xmin=args.trange[0], xmax=args.trange[1]) + + t1=time.time() + print("elapsed time ",t1-t0) + print("Optimal azimuth for trange between " + str(args.trange[0])+" and "+str(args.trange[1]) + " seconds is: "+str(harmonics.azim)) From 30b0192110aac1442c8509f1af548786b99314d4 Mon Sep 17 00:00:00 2001 From: Anang Sahroni Date: Sat, 25 Jun 2022 14:03:03 +0700 Subject: [PATCH 02/12] move parameter preparation to the main harmonics file and add wrapper function --- rfpy/harmonics.py | 54 ++++++++++++++++++++++++++++++---- rfpy/scripts/rfpy_harmonics.py | 41 ++------------------------ 2 files changed, 51 insertions(+), 44 deletions(-) diff --git a/rfpy/harmonics.py b/rfpy/harmonics.py index 8d54949..4581a29 100644 --- a/rfpy/harmonics.py +++ b/rfpy/harmonics.py @@ -109,7 +109,7 @@ def __init__(self, radialRF, transvRF=None, azim=0, xmin=0., xmax=10.): @staticmethod @numba.jit(nopython=True) - def dcomp_find_azim_numba(xmin=None, xmax=None, + def _dcomp_find_azim_numba(xmin=None, xmax=None, nbin=None, nz=None, delta=None, baz_list=None, dataR_list=None, @@ -141,11 +141,6 @@ def dcomp_find_azim_numba(xmin=None, xmax=None, """ - if not xmin: - xmin = xmin - if not xmax: - xmax = xmax - print() print('Decomposing receiver functions into baz harmonics') @@ -230,6 +225,48 @@ def dcomp_find_azim_numba(xmin=None, xmax=None, return C0,C1,C2,C3,C4,C0var,C1var,C2var,C3var,C4var,indaz,daz + def dcomp_find_azim_numba(self, xmin=None, xmax=None): + import time + t0 = time.time() + + if not xmin: + xmin = self.xmin + if not xmax: + xmax = self.xmax + + str_stats = self.radialRF[0].stats + baz_list = np.array([trace0.stats.baz for trace0 in self.radialRF]) + dataR_list = np.array([trace0.data for trace0 in self.radialRF]) + dataT_list = np.array([trace0.data for trace0 in self.transvRF]) + C0,C1,C2,C3,C4,C0var,C1var,C2var,C3var,C4var,indaz,daz = self._dcomp_find_azim_numba( + xmin=xmin, xmax=xmax, + nbin=len(self.radialRF), + nz=len(self.radialRF[0].data), + delta=self.radialRF[0].stats.delta, + baz_list=baz_list, + dataR_list=dataR_list, + dataT_list=dataT_list, + ) + + # Put back into traces + A = Trace(data=C0[:, indaz], header=str_stats) + B1 = Trace(data=C1[:, indaz], header=str_stats) + B2 = Trace(data=C2[:, indaz], header=str_stats) + C1 = Trace(data=C3[:, indaz], header=str_stats) + C2 = Trace(data=C4[:, indaz], header=str_stats) + + # Put all treaces into stream + hstream = Stream(traces=[A, B1, B2, C1, C2]) + azim = indaz*daz + var = [C0var, C1var, C2var, C3var, C4var] + + self.hstream = hstream + self.azim = azim + self.var = var + + t1 = time.time() + print("Elapsed time ", t1-t0) + def dcomp_find_azim(self, xmin=None, xmax=None): """ Method to decompose radial and transverse receiver function @@ -255,6 +292,8 @@ def dcomp_find_azim(self, xmin=None, xmax=None): Variance of the 5 harmonics between ``xmin`` and ``xmax`` """ + import time + t0 = time.time() if not xmin: xmin = self.xmin @@ -358,6 +397,9 @@ def dcomp_find_azim(self, xmin=None, xmax=None): self.azim = indaz*daz self.var = [C0var, C1var, C2var, C3var, C4var] + t1 = time.time() + print("Elapsed time ", t1-t0) + def dcomp_fix_azim(self, azim=None): """ Method to decompose radial and transverse receiver function diff --git a/rfpy/scripts/rfpy_harmonics.py b/rfpy/scripts/rfpy_harmonics.py index 9a805ff..a62cccb 100644 --- a/rfpy/scripts/rfpy_harmonics.py +++ b/rfpy/scripts/rfpy_harmonics.py @@ -91,7 +91,7 @@ def get_harmonics_arguments(argv=None): parser.add_argument( "-UN", "--use-numba", action="store_true", - dest="unumba", + dest="use_numba", default=False, help="Use numba jit on calculation [Default False]" ) @@ -510,46 +510,11 @@ def main(): # Stack with or without dip if args.find_azim: - import time - t0=time.time() - - if args.unumba: - from obspy.core import Trace - str_stats = harmonics.radialRF[0].stats - baz_list = np.array([trace0.stats.baz for trace0 in harmonics.radialRF]) - dataR_list = np.array([trace0.data for trace0 in harmonics.radialRF]) - dataT_list = np.array([trace0.data for trace0 in harmonics.transvRF]) - C0,C1,C2,C3,C4,C0var,C1var,C2var,C3var,C4var,indaz,daz = harmonics.dcomp_find_azim_numba( - xmin=args.trange[0], xmax=args.trange[1], - nbin=len(harmonics.radialRF), - nz=len(harmonics.radialRF[0].data), - delta=harmonics.radialRF[0].stats.delta, - baz_list=baz_list, - dataR_list=dataR_list, - dataT_list=dataT_list, - - ) - # Put back into traces - A = Trace(data=C0[:, indaz], header=str_stats) - B1 = Trace(data=C1[:, indaz], header=str_stats) - B2 = Trace(data=C2[:, indaz], header=str_stats) - C1 = Trace(data=C3[:, indaz], header=str_stats) - C2 = Trace(data=C4[:, indaz], header=str_stats) - - # Put all treaces into stream - hstream = Stream(traces=[A, B1, B2, C1, C2]) - azim = indaz*daz - var = [C0var, C1var, C2var, C3var, C4var] - - harmonics.hstream = hstream - harmonics.azim = azim - harmonics.var = var + if args.use_numba: + harmonics.dcomp_find_azim_numba(xmin=args.trange[0], xmax=args.trange[1]) else: harmonics.dcomp_find_azim(xmin=args.trange[0], xmax=args.trange[1]) - t1=time.time() - print("elapsed time ",t1-t0) - print("Optimal azimuth for trange between " + str(args.trange[0])+" and "+str(args.trange[1]) + " seconds is: "+str(harmonics.azim)) From ded2583318504c350de43c379571251045ae6370 Mon Sep 17 00:00:00 2001 From: Anang Sahroni Date: Sat, 25 Jun 2022 14:04:26 +0700 Subject: [PATCH 03/12] add numba to dependencies, not specified by version --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index c9f78dc..fc1e61d 100644 --- a/setup.py +++ b/setup.py @@ -32,7 +32,7 @@ def find_version(*paths): 'Programming Language :: Python :: 3.7', 'Programming Language :: Python :: 3.8', 'Programming Language :: Python :: 3.9'], - install_requires=['numpy', 'obspy', 'stdb>=0.2.0'], + install_requires=['numpy', 'numba', 'obspy', 'stdb>=0.2.0'], python_requires='>=3.6', packages=setuptools.find_packages(), entry_points={'console_scripts':[ From 09881e6fbb1b29a7a5a24f217a84afeac4d91c8f Mon Sep 17 00:00:00 2001 From: Anang Sahroni Date: Sat, 25 Jun 2022 14:26:50 +0700 Subject: [PATCH 04/12] add temporary progressbar --- rfpy/harmonics.py | 28 +++++++++++++++++----------- 1 file changed, 17 insertions(+), 11 deletions(-) diff --git a/rfpy/harmonics.py b/rfpy/harmonics.py index 4581a29..bec886e 100644 --- a/rfpy/harmonics.py +++ b/rfpy/harmonics.py @@ -29,6 +29,7 @@ from obspy.core import Stream, Trace import matplotlib.pyplot as plt import numba +from numba_progress import ProgressBar class Harmonics(object): @@ -109,12 +110,13 @@ def __init__(self, radialRF, transvRF=None, azim=0, xmin=0., xmax=10.): @staticmethod @numba.jit(nopython=True) - def _dcomp_find_azim_numba(xmin=None, xmax=None, + def _dcomp_find_azim_numba( + progress_proxy, + xmin=None, xmax=None, nbin=None, nz=None, delta=None, baz_list=None, dataR_list=None, dataT_list=None, - ): """ Method to decompose radial and transverse receiver function @@ -210,6 +212,8 @@ def _dcomp_find_azim_numba(xmin=None, xmax=None, C3[iz, iaz] = np.float(CC[3]) C4[iz, iaz] = np.float(CC[4]) + progress_proxy.update(1) + # Minimize variance of third component over specific depth range to # find azim C1var = np.zeros(naz) @@ -238,15 +242,17 @@ def dcomp_find_azim_numba(self, xmin=None, xmax=None): baz_list = np.array([trace0.stats.baz for trace0 in self.radialRF]) dataR_list = np.array([trace0.data for trace0 in self.radialRF]) dataT_list = np.array([trace0.data for trace0 in self.transvRF]) - C0,C1,C2,C3,C4,C0var,C1var,C2var,C3var,C4var,indaz,daz = self._dcomp_find_azim_numba( - xmin=xmin, xmax=xmax, - nbin=len(self.radialRF), - nz=len(self.radialRF[0].data), - delta=self.radialRF[0].stats.delta, - baz_list=baz_list, - dataR_list=dataR_list, - dataT_list=dataT_list, - ) + with ProgressBar(total=len(self.radialRF[0].data)) as progress: + C0,C1,C2,C3,C4,C0var,C1var,C2var,C3var,C4var,indaz,daz = self._dcomp_find_azim_numba( + xmin=xmin, xmax=xmax, + nbin=len(self.radialRF), + nz=len(self.radialRF[0].data), + delta=self.radialRF[0].stats.delta, + baz_list=baz_list, + dataR_list=dataR_list, + dataT_list=dataT_list, + progress_proxy=progress + ) # Put back into traces A = Trace(data=C0[:, indaz], header=str_stats) From 785deca2ecef2518fa20779b65bbc14ff55d6dd9 Mon Sep 17 00:00:00 2001 From: Anang Sahroni Date: Sat, 25 Jun 2022 15:45:06 +0700 Subject: [PATCH 05/12] add numba_progressbar --- rfpy/harmonics.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/rfpy/harmonics.py b/rfpy/harmonics.py index bec886e..ef84164 100644 --- a/rfpy/harmonics.py +++ b/rfpy/harmonics.py @@ -109,9 +109,9 @@ def __init__(self, radialRF, transvRF=None, azim=0, xmin=0., xmax=10.): self.xmax = xmax @staticmethod - @numba.jit(nopython=True) + @numba.jit(nopython=True, nogil=True) def _dcomp_find_azim_numba( - progress_proxy, + progress_proxy=None, xmin=None, xmax=None, nbin=None, nz=None, delta=None, baz_list=None, @@ -143,8 +143,7 @@ def _dcomp_find_azim_numba( """ - print() - print('Decomposing receiver functions into baz harmonics') + # Some integers nbin = nbin @@ -230,6 +229,10 @@ def _dcomp_find_azim_numba( return C0,C1,C2,C3,C4,C0var,C1var,C2var,C3var,C4var,indaz,daz def dcomp_find_azim_numba(self, xmin=None, xmax=None): + + print() + print('Decomposing receiver functions into baz harmonics') + import time t0 = time.time() @@ -242,7 +245,7 @@ def dcomp_find_azim_numba(self, xmin=None, xmax=None): baz_list = np.array([trace0.stats.baz for trace0 in self.radialRF]) dataR_list = np.array([trace0.data for trace0 in self.radialRF]) dataT_list = np.array([trace0.data for trace0 in self.transvRF]) - with ProgressBar(total=len(self.radialRF[0].data)) as progress: + with ProgressBar(total=len(self.radialRF[0].data), ascii=" #") as progress: C0,C1,C2,C3,C4,C0var,C1var,C2var,C3var,C4var,indaz,daz = self._dcomp_find_azim_numba( xmin=xmin, xmax=xmax, nbin=len(self.radialRF), From 94fbe472d016d69b512d57ab32f7a40643fcbdb0 Mon Sep 17 00:00:00 2001 From: Anang Sahroni Date: Sat, 25 Jun 2022 19:03:51 +0700 Subject: [PATCH 06/12] add numba_progress to dependency --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index fc1e61d..299beec 100644 --- a/setup.py +++ b/setup.py @@ -32,7 +32,7 @@ def find_version(*paths): 'Programming Language :: Python :: 3.7', 'Programming Language :: Python :: 3.8', 'Programming Language :: Python :: 3.9'], - install_requires=['numpy', 'numba', 'obspy', 'stdb>=0.2.0'], + install_requires=['numpy','obspy', 'stdb>=0.2.0', 'numba', 'numba_progress'], python_requires='>=3.6', packages=setuptools.find_packages(), entry_points={'console_scripts':[ From d9a9accefaf7fb35fd46a78b392be396016179f1 Mon Sep 17 00:00:00 2001 From: Anang Sahroni Date: Sat, 25 Jun 2022 19:04:28 +0700 Subject: [PATCH 07/12] move function to the main dcomp method and tidy up the code --- rfpy/harmonics.py | 349 ++++++++++++++++----------------- rfpy/scripts/rfpy_harmonics.py | 18 +- 2 files changed, 176 insertions(+), 191 deletions(-) diff --git a/rfpy/harmonics.py b/rfpy/harmonics.py index ef84164..9f3b72c 100644 --- a/rfpy/harmonics.py +++ b/rfpy/harmonics.py @@ -108,16 +108,7 @@ def __init__(self, radialRF, transvRF=None, azim=0, xmin=0., xmax=10.): self.xmin = xmin self.xmax = xmax - @staticmethod - @numba.jit(nopython=True, nogil=True) - def _dcomp_find_azim_numba( - progress_proxy=None, - xmin=None, xmax=None, - nbin=None, nz=None, delta=None, - baz_list=None, - dataR_list=None, - dataT_list=None, - ): + def dcomp_find_azim(self, xmin=None, xmax=None, use_numba=False): """ Method to decompose radial and transverse receiver function streams into back-azimuth harmonics and determine the main @@ -130,6 +121,8 @@ def _dcomp_find_azim_numba( Minimum x axis value over which to calculate ``azim`` xmax : float Maximum x axis value over which to calculate ``azim`` + use_numba : bool + Use ``numba`` just-in-time compiler to increase the processing speed, default to ``False`` Attributes ---------- @@ -142,120 +135,117 @@ def _dcomp_find_azim_numba( Variance of the 5 harmonics between ``xmin`` and ``xmax`` """ + import time + t0 = time.time() + + if not xmin: + xmin = self.xmin + if not xmax: + xmax = self.xmax - + print() + print('Decomposing receiver functions into baz harmonics') # Some integers - nbin = nbin - nz = nz + nbin = len(self.radialRF) + nz = len(self.radialRF[0].data) naz = 180 daz = np.float(360/naz) deg2rad = np.pi/180. # Define depth range over which to calculate azimuth - indmin = int(xmin/delta) - indmax = int(xmax/delta) + indmin = int(xmin/self.radialRF[0].stats.delta) + indmax = int(xmax/self.radialRF[0].stats.delta) - # Initialize work arrays - C0 = np.zeros((nz, naz)) - C1 = np.zeros((nz, naz)) - C2 = np.zeros((nz, naz)) - C3 = np.zeros((nz, naz)) - C4 = np.zeros((nz, naz)) + # Copy stream stats + str_stats = self.radialRF[0].stats - # Loop over each depth step - for iz in range(nz): + if use_numba == False: + # Initialize work arrays + C0 = np.zeros((nz, naz)) + C1 = np.zeros((nz, naz)) + C2 = np.zeros((nz, naz)) + C3 = np.zeros((nz, naz)) + C4 = np.zeros((nz, naz)) + # Loop over each depth step + for iz in range(nz): - # Build matrices OBS and H for each azimuth + # Build matrices OBS and H for each azimuth + for iaz in range(naz): + + # Initialize work arrays + OBS = np.zeros(2*nbin) + H = np.zeros((2*nbin, 5)) + + azim = iaz*daz + + # Radial component + for irow, trace in enumerate(self.radialRF): + + baz = trace.stats.baz + OBS[irow] = trace.data[iz] + H[irow, 0] = 1.0 + H[irow, 1] = np.cos(deg2rad*(baz-azim)) + H[irow, 2] = np.sin(deg2rad*(baz-azim)) + H[irow, 3] = np.cos(2.*deg2rad*(baz-azim)) + H[irow, 4] = np.sin(2.*deg2rad*(baz-azim)) + + shift = 90. + + # Transverse component + for irow, trace in enumerate(self.transvRF): + + baz = trace.stats.baz + OBS[irow+nbin] = trace.data[iz] + H[irow+nbin, 0] = 0.0 + H[irow+nbin, 1] = np.cos(deg2rad*(baz+shift-azim)) + H[irow+nbin, 2] = np.sin(deg2rad*(baz+shift-azim)) + H[irow+nbin, 3] = np.cos(2.*deg2rad*(baz+shift/2.0-azim)) + H[irow+nbin, 4] = np.sin(2.*deg2rad*(baz+shift/2.0-azim)) + + # Solve system of equations with truncated SVD + u, s, v = np.linalg.svd(H) + s[s < 0.001] = 0. + CC = np.linalg.solve(s[:, None] * v, u.T.dot(OBS)[:5]) + + # Fill up arrays + C0[iz, iaz] = np.float(CC[0]) + C1[iz, iaz] = np.float(CC[1]) + C2[iz, iaz] = np.float(CC[2]) + C3[iz, iaz] = np.float(CC[3]) + C4[iz, iaz] = np.float(CC[4]) + + # Minimize variance of third component over specific depth range to + # find azim + C1var = np.zeros(naz) for iaz in range(naz): - - # Initialize work arrays - OBS = np.zeros(2*nbin) - H = np.zeros((2*nbin, 5)) - - azim = iaz*daz - - # Radial component - for irow, dataR, baz in zip(range(len(dataR_list)),dataR_list,baz_list): - - baz = baz - OBS[irow] = dataR[iz] - H[irow, 0] = 1.0 - H[irow, 1] = np.cos(deg2rad*(baz-azim)) - H[irow, 2] = np.sin(deg2rad*(baz-azim)) - H[irow, 3] = np.cos(2.*deg2rad*(baz-azim)) - H[irow, 4] = np.sin(2.*deg2rad*(baz-azim)) - - shift = 90. - - # Transverse component - for irow, dataT, baz in zip(range(len(dataT_list)),dataT_list,baz_list): - - baz = baz - OBS[irow+nbin] = dataT[iz] - H[irow+nbin, 0] = 0.0 - H[irow+nbin, 1] = np.cos(deg2rad*(baz+shift-azim)) - H[irow+nbin, 2] = np.sin(deg2rad*(baz+shift-azim)) - H[irow+nbin, 3] = np.cos(2.*deg2rad*(baz+shift/2.0-azim)) - H[irow+nbin, 4] = np.sin(2.*deg2rad*(baz+shift/2.0-azim)) - - # Solve system of equations with truncated SVD - u, s, v = np.linalg.svd(H) - s[s < 0.001] = 0. - CC = np.linalg.solve(s.reshape(s.shape[0],1) * v, u.T.dot(OBS)[:5]) - - # Fill up arrays - C0[iz, iaz] = np.float(CC[0]) - C1[iz, iaz] = np.float(CC[1]) - C2[iz, iaz] = np.float(CC[2]) - C3[iz, iaz] = np.float(CC[3]) - C4[iz, iaz] = np.float(CC[4]) - - progress_proxy.update(1) - - # Minimize variance of third component over specific depth range to - # find azim - C1var = np.zeros(naz) - for iaz in range(naz): - C1var[iaz] = np.sqrt(np.mean(np.square(C1[indmin:indmax, iaz]))) - indaz = np.argmin(C1var) - - C0var = np.sqrt(np.mean(np.square(C0[indmin:indmax, indaz]))) - C1var = np.sqrt(np.mean(np.square(C1[indmin:indmax, indaz]))) - C2var = np.sqrt(np.mean(np.square(C2[indmin:indmax, indaz]))) - C3var = np.sqrt(np.mean(np.square(C3[indmin:indmax, indaz]))) - C4var = np.sqrt(np.mean(np.square(C4[indmin:indmax, indaz]))) - - return C0,C1,C2,C3,C4,C0var,C1var,C2var,C3var,C4var,indaz,daz - - def dcomp_find_azim_numba(self, xmin=None, xmax=None): - - print() - print('Decomposing receiver functions into baz harmonics') - - import time - t0 = time.time() - - if not xmin: - xmin = self.xmin - if not xmax: - xmax = self.xmax - - str_stats = self.radialRF[0].stats - baz_list = np.array([trace0.stats.baz for trace0 in self.radialRF]) - dataR_list = np.array([trace0.data for trace0 in self.radialRF]) - dataT_list = np.array([trace0.data for trace0 in self.transvRF]) - with ProgressBar(total=len(self.radialRF[0].data), ascii=" #") as progress: - C0,C1,C2,C3,C4,C0var,C1var,C2var,C3var,C4var,indaz,daz = self._dcomp_find_azim_numba( - xmin=xmin, xmax=xmax, - nbin=len(self.radialRF), - nz=len(self.radialRF[0].data), - delta=self.radialRF[0].stats.delta, - baz_list=baz_list, - dataR_list=dataR_list, - dataT_list=dataT_list, - progress_proxy=progress - ) + C1var[iaz] = np.sqrt(np.mean(np.square(C1[indmin:indmax, iaz]))) + indaz = np.argmin(C1var) + + C0var = np.sqrt(np.mean(np.square(C0[indmin:indmax, indaz]))) + C1var = np.sqrt(np.mean(np.square(C1[indmin:indmax, indaz]))) + C2var = np.sqrt(np.mean(np.square(C2[indmin:indmax, indaz]))) + C3var = np.sqrt(np.mean(np.square(C3[indmin:indmax, indaz]))) + C4var = np.sqrt(np.mean(np.square(C4[indmin:indmax, indaz]))) + + elif use_numba == True: + # convert stream objects to numpy array + baz_arr = np.array([trace0.stats.baz for trace0 in self.radialRF]) + radialRF_arr = np.array([trace0.data for trace0 in self.radialRF]) + transvRF_arr = np.array([trace0.data for trace0 in self.transvRF]) + # run calculation using numba + with ProgressBar(total=len(self.radialRF[0].data), ascii=" #") as progress: + C0,C1,C2,C3,C4,C0var,C1var,C2var,C3var,C4var,indaz = \ + self._dcomp_calculate_harmonics_isolated( + nbin=nbin, + nz=nz, + indmin=indmin, + indmax=indmax, + baz_arr=baz_arr, + radialRF_arr=radialRF_arr, + transvRF_arr=transvRF_arr, + progress_hook=progress + ) # Put back into traces A = Trace(data=C0[:, indaz], header=str_stats) @@ -265,67 +255,73 @@ def dcomp_find_azim_numba(self, xmin=None, xmax=None): C2 = Trace(data=C4[:, indaz], header=str_stats) # Put all treaces into stream - hstream = Stream(traces=[A, B1, B2, C1, C2]) - azim = indaz*daz - var = [C0var, C1var, C2var, C3var, C4var] - - self.hstream = hstream - self.azim = azim - self.var = var + self.hstream = Stream(traces=[A, B1, B2, C1, C2]) + self.azim = indaz*daz + self.var = [C0var, C1var, C2var, C3var, C4var] t1 = time.time() print("Elapsed time ", t1-t0) - def dcomp_find_azim(self, xmin=None, xmax=None): + @staticmethod + @numba.jit(nopython=True, nogil=True) + def _dcomp_calculate_harmonics_isolated(nbin=None, nz=None, baz_arr=None, \ + radialRF_arr=None, transvRF_arr=None, indmin=None, indmax=None, progress_hook=None): """ Method to decompose radial and transverse receiver function - streams into back-azimuth harmonics and determine the main - orientation ``azim``, obtained by minimizing the B1 component - between ``xmin`` and ``xmax`` (i.e., time or depth). + array into back-azimuth harmonics. This static method is isolated from + ``decomp_find_azim`` to accomodate ``numba`` support that requires ```numpy```-like objects + when compiling to machine code. Parameters ---------- - xmin : float - Minimum x axis value over which to calculate ``azim`` - xmax : float - Maximum x axis value over which to calculate ``azim`` + nbin : integer + Number of receiver function traces + nz : integer + Number of receiver function trace's points, corresponding to depth or time + indmax : integer + Index of maximum depth + indmin : integer + Index of minimum depth + baz_arr : :class:`~numpy.array` + Array containing each receiver function's back-azimuth + radialRF : :class:`~numpy.ndarray` + Array containing radial receiver functions + transvRF : :class:`~numpy.ndarray` + Array containing transversal receiver functions + progress_hook : :class:`~numba_progress.ProgressBar` + Progressbar hook to integrate ``numba_progress`` Attributes ---------- - hstream : :class:`~obspy.core.Stream` - Stream containing the 5 harmonics, oriented in direction ``azim`` - azim : float - Direction (azimuth) along which the B1 component of the stream - is minimized (between ``xmin`` and ``xmax``) - var : :class:`~numpy.ndarray` - Variance of the 5 harmonics between ``xmin`` and ``xmax`` - + C0 : :class:`~numpy.ndarray` + Array of C0 harmonic + C1 : :class:`~numpy.ndarray` + Array of C1 harmonic + C2 : :class:`~numpy.ndarray` + Array of C2 harmonic + C3 : :class:`~numpy.ndarray` + Array of C3 harmonic + C4 : :class:`~numpy.ndarray` + Array of C4 harmonic + C0var : :class:`~numpy.ndarray` + Array of C0 harmonic variance + C1var : :class:`~numpy.ndarray` + Array of C1 harmonic variance + C2var : :class:`~numpy.ndarray` + Array of C2 harmonic variance + C3var : :class:`~numpy.ndarray` + Array of C3 harmonic variance + C4var : :class:`~numpy.ndarray` + Array of C4 harmonic variance + indaz : integer + Index of minimum ``C1Var`` """ - import time - t0 = time.time() - - if not xmin: - xmin = self.xmin - if not xmax: - xmax = self.xmax - - print() - print('Decomposing receiver functions into baz harmonics') # Some integers - nbin = len(self.radialRF) - nz = len(self.radialRF[0].data) naz = 180 daz = np.float(360/naz) deg2rad = np.pi/180. - # Define depth range over which to calculate azimuth - indmin = int(xmin/self.radialRF[0].stats.delta) - indmax = int(xmax/self.radialRF[0].stats.delta) - - # Copy stream stats - str_stats = self.radialRF[0].stats - # Initialize work arrays C0 = np.zeros((nz, naz)) C1 = np.zeros((nz, naz)) @@ -346,10 +342,10 @@ def dcomp_find_azim(self, xmin=None, xmax=None): azim = iaz*daz # Radial component - for irow, trace in enumerate(self.radialRF): + for irow, dataR, baz in zip(range(len(radialRF_arr)),radialRF_arr,baz_arr): - baz = trace.stats.baz - OBS[irow] = trace.data[iz] + baz = baz + OBS[irow] = dataR[iz] H[irow, 0] = 1.0 H[irow, 1] = np.cos(deg2rad*(baz-azim)) H[irow, 2] = np.sin(deg2rad*(baz-azim)) @@ -359,10 +355,10 @@ def dcomp_find_azim(self, xmin=None, xmax=None): shift = 90. # Transverse component - for irow, trace in enumerate(self.transvRF): + for irow, dataT, baz in zip(range(len(transvRF_arr)),transvRF_arr,baz_arr): - baz = trace.stats.baz - OBS[irow+nbin] = trace.data[iz] + baz = baz + OBS[irow+nbin] = dataT[iz] H[irow+nbin, 0] = 0.0 H[irow+nbin, 1] = np.cos(deg2rad*(baz+shift-azim)) H[irow+nbin, 2] = np.sin(deg2rad*(baz+shift-azim)) @@ -372,7 +368,7 @@ def dcomp_find_azim(self, xmin=None, xmax=None): # Solve system of equations with truncated SVD u, s, v = np.linalg.svd(H) s[s < 0.001] = 0. - CC = np.linalg.solve(s[:, None] * v, u.T.dot(OBS)[:5]) + CC = np.linalg.solve(s.reshape(s.shape[0],1) * v, u.T.dot(OBS)[:5]) # Fill up arrays C0[iz, iaz] = np.float(CC[0]) @@ -381,6 +377,8 @@ def dcomp_find_azim(self, xmin=None, xmax=None): C3[iz, iaz] = np.float(CC[3]) C4[iz, iaz] = np.float(CC[4]) + progress_hook.update(1) + # Minimize variance of third component over specific depth range to # find azim C1var = np.zeros(naz) @@ -394,20 +392,7 @@ def dcomp_find_azim(self, xmin=None, xmax=None): C3var = np.sqrt(np.mean(np.square(C3[indmin:indmax, indaz]))) C4var = np.sqrt(np.mean(np.square(C4[indmin:indmax, indaz]))) - # Put back into traces - A = Trace(data=C0[:, indaz], header=str_stats) - B1 = Trace(data=C1[:, indaz], header=str_stats) - B2 = Trace(data=C2[:, indaz], header=str_stats) - C1 = Trace(data=C3[:, indaz], header=str_stats) - C2 = Trace(data=C4[:, indaz], header=str_stats) - - # Put all treaces into stream - self.hstream = Stream(traces=[A, B1, B2, C1, C2]) - self.azim = indaz*daz - self.var = [C0var, C1var, C2var, C3var, C4var] - - t1 = time.time() - print("Elapsed time ", t1-t0) + return C0,C1,C2,C3,C4,C0var,C1var,C2var,C3var,C4var,indaz def dcomp_fix_azim(self, azim=None): """ @@ -503,7 +488,7 @@ def dcomp_fix_azim(self, azim=None): # Put all traces into stream self.hstream = Stream(traces=[A, B1, B2, C1, C2]) - def forward(self, baz_list=None): + def forward(self, baz_arr=None): """ Method to forward calculate radial and transverse component receiver functions given the 5 pre-determined harmonics and @@ -513,7 +498,7 @@ def forward(self, baz_list=None): Parameters ---------- - baz_list : list + baz_arr : list List of back-azimuth directions over which to calculate the receiver functions. If no list is specified, the method will use the same back-azimuths as those in the original @@ -532,12 +517,12 @@ def forward(self, baz_list=None): if not hasattr(self, 'hstream'): raise(Exception("Decomposition has not been performed yet")) - if not baz_list: + if not baz_arr: print("Warning: no BAZ specified - using all baz from " + "stored streams") - baz_list = [tr.stats.baz for tr in self.radialRF] - if not isinstance(baz_list, list): - baz_list = [baz_list] + baz_arr = [tr.stats.baz for tr in self.radialRF] + if not isinstance(baz_arr, list): + baz_arr = [baz_arr] # Some constants nz = len(self.hstream[0].data) @@ -547,7 +532,7 @@ def forward(self, baz_list=None): self.radial_forward = Stream() self.transv_forward = Stream() - for baz in baz_list: + for baz in baz_arr: trR = Trace(header=self.hstream[0].stats) trT = Trace(header=self.hstream[0].stats) diff --git a/rfpy/scripts/rfpy_harmonics.py b/rfpy/scripts/rfpy_harmonics.py index a62cccb..15d283f 100644 --- a/rfpy/scripts/rfpy_harmonics.py +++ b/rfpy/scripts/rfpy_harmonics.py @@ -88,14 +88,6 @@ def get_harmonics_arguments(argv=None): "Default behaviour uses short key form (NET.STN) for the folder " + "names, regardless of the key type of the database." ) - parser.add_argument( - "-UN", "--use-numba", - action="store_true", - dest="use_numba", - default=False, - help="Use numba jit on calculation [Default False]" - ) - # Event Selection Criteria TimeGroup = parser.add_argument_group( @@ -217,6 +209,13 @@ def get_harmonics_arguments(argv=None): default=False, help="Set this option to save the Harmonics object " + "to a pickled file. [Default does not save object]") + HarmonicGroup.add_argument( + "--use-numba", + action="store_true", + dest="use_numba", + default=False, + help="Use numba jit to increase processing speed [Default does not use numba]" + ) PlotGroup = parser.add_argument_group( title='Settings for plotting results', @@ -510,8 +509,9 @@ def main(): # Stack with or without dip if args.find_azim: + # Check if the process use numba if args.use_numba: - harmonics.dcomp_find_azim_numba(xmin=args.trange[0], xmax=args.trange[1]) + harmonics.dcomp_find_azim(xmin=args.trange[0], xmax=args.trange[1], use_numba=True) else: harmonics.dcomp_find_azim(xmin=args.trange[0], xmax=args.trange[1]) From 1b4dfd1574d02ac31cb8a8c026f3e46c5558d18c Mon Sep 17 00:00:00 2001 From: Anang Sahroni Date: Sat, 25 Jun 2022 19:10:16 +0700 Subject: [PATCH 08/12] add progressbar to standard harmonics --- rfpy/harmonics.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/rfpy/harmonics.py b/rfpy/harmonics.py index 9f3b72c..8625d3b 100644 --- a/rfpy/harmonics.py +++ b/rfpy/harmonics.py @@ -30,6 +30,7 @@ import matplotlib.pyplot as plt import numba from numba_progress import ProgressBar +from tqdm import trange class Harmonics(object): @@ -168,7 +169,7 @@ def dcomp_find_azim(self, xmin=None, xmax=None, use_numba=False): C3 = np.zeros((nz, naz)) C4 = np.zeros((nz, naz)) # Loop over each depth step - for iz in range(nz): + for iz in trange(nz, ascii=True): # Build matrices OBS and H for each azimuth for iaz in range(naz): From bfd5f576f0e1f649da486347a56c32e2b63a71f1 Mon Sep 17 00:00:00 2001 From: Anang Sahroni Date: Sat, 25 Jun 2022 19:20:40 +0700 Subject: [PATCH 09/12] add use_numba to the documentation --- docs/rfpy.rst | 6 ++++++ docs/scripts.rst | 2 ++ docs/tutorials.rst | 9 ++++++++- 3 files changed, 16 insertions(+), 1 deletion(-) diff --git a/docs/rfpy.rst b/docs/rfpy.rst index ab1cb31..cfc546e 100644 --- a/docs/rfpy.rst +++ b/docs/rfpy.rst @@ -654,6 +654,12 @@ Or, alternatively, >>> harmonics.dcomp_find_azim() +To use `numba` compiler to increase speed, + +.. sourcecode:: python + + >>> harmonics.dcomp_find_azim(use_numba=True) + In either case the harmonic components are available as an attribute of type :class:`~obspy.core.Stream` (``harmonics.hstream``) and, if available, the azimuth of the dominant direction (``harmonics.azim``). diff --git a/docs/scripts.rst b/docs/scripts.rst index b353306..c4807fb 100644 --- a/docs/scripts.rst +++ b/docs/scripts.rst @@ -578,6 +578,8 @@ Usage set] --save Set this option to save the Harmonics object to a pickled file. [Default does not save object] + --use-numba Use numba jit to increase processing speed + [Default does not use numba] Settings for plotting results: Specify parameters for plotting the back-azimuth harmonics. diff --git a/docs/tutorials.rst b/docs/tutorials.rst index bbbc06c..6d2cad5 100644 --- a/docs/tutorials.rst +++ b/docs/tutorials.rst @@ -324,12 +324,19 @@ to 10 seconds (to avoid the large zero-lag pulse): .. note:: - Warning!! This command is particularly slow, especially for large data sets. + Warning!! This command is particularly slow, especially for large data sets. The speed can + be increase by using ``--use-numba`` argument. .. code-block:: $ rfpy_harmonics --no-outlier --find-azim --trange=2.,10. MMPY.pkl +To increase the processing speed use ``--use_numba`` to utilize ``numba`` just-in-time compiler: + +.. code-block:: + + $ rfpy_harmonics --use_numba --no-outlier --find-azim --trange=2.,10. MMPY.pkl + ################################################################################ # __ _ _ # # _ __ / _|_ __ _ _ | |__ __ _ _ __ _ __ ___ ___ _ __ (_) ___ ___ # From 5c0146a5dcd3117582648403ba0cc7f244c39234 Mon Sep 17 00:00:00 2001 From: Anang Sahroni Date: Sun, 26 Jun 2022 06:26:40 +0700 Subject: [PATCH 10/12] delete time testing print statements --- rfpy/harmonics.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/rfpy/harmonics.py b/rfpy/harmonics.py index 8625d3b..c92a07b 100644 --- a/rfpy/harmonics.py +++ b/rfpy/harmonics.py @@ -136,8 +136,6 @@ def dcomp_find_azim(self, xmin=None, xmax=None, use_numba=False): Variance of the 5 harmonics between ``xmin`` and ``xmax`` """ - import time - t0 = time.time() if not xmin: xmin = self.xmin @@ -260,9 +258,6 @@ def dcomp_find_azim(self, xmin=None, xmax=None, use_numba=False): self.azim = indaz*daz self.var = [C0var, C1var, C2var, C3var, C4var] - t1 = time.time() - print("Elapsed time ", t1-t0) - @staticmethod @numba.jit(nopython=True, nogil=True) def _dcomp_calculate_harmonics_isolated(nbin=None, nz=None, baz_arr=None, \ From 42b6859012b7d1d48788eb2b80428de916d10f7c Mon Sep 17 00:00:00 2001 From: Anang Sahroni Date: Sun, 26 Jun 2022 14:36:49 +0700 Subject: [PATCH 11/12] revert accidentally change variable name for baz_list --- rfpy/harmonics.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/rfpy/harmonics.py b/rfpy/harmonics.py index c92a07b..fb8a399 100644 --- a/rfpy/harmonics.py +++ b/rfpy/harmonics.py @@ -484,7 +484,7 @@ def dcomp_fix_azim(self, azim=None): # Put all traces into stream self.hstream = Stream(traces=[A, B1, B2, C1, C2]) - def forward(self, baz_arr=None): + def forward(self, baz_list=None): """ Method to forward calculate radial and transverse component receiver functions given the 5 pre-determined harmonics and @@ -494,7 +494,7 @@ def forward(self, baz_arr=None): Parameters ---------- - baz_arr : list + baz_list : list List of back-azimuth directions over which to calculate the receiver functions. If no list is specified, the method will use the same back-azimuths as those in the original @@ -513,12 +513,12 @@ def forward(self, baz_arr=None): if not hasattr(self, 'hstream'): raise(Exception("Decomposition has not been performed yet")) - if not baz_arr: + if not baz_list: print("Warning: no BAZ specified - using all baz from " + "stored streams") - baz_arr = [tr.stats.baz for tr in self.radialRF] - if not isinstance(baz_arr, list): - baz_arr = [baz_arr] + baz_list = [tr.stats.baz for tr in self.radialRF] + if not isinstance(baz_list, list): + baz_list = [baz_list] # Some constants nz = len(self.hstream[0].data) @@ -528,7 +528,7 @@ def forward(self, baz_arr=None): self.radial_forward = Stream() self.transv_forward = Stream() - for baz in baz_arr: + for baz in baz_list: trR = Trace(header=self.hstream[0].stats) trT = Trace(header=self.hstream[0].stats) From 61f86fa942c941d962ae37dd7288623a4e93d107 Mon Sep 17 00:00:00 2001 From: Anang Sahroni Date: Mon, 27 Jun 2022 17:24:19 +0700 Subject: [PATCH 12/12] add warning while using numba --- docs/tutorials.rst | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/docs/tutorials.rst b/docs/tutorials.rst index 6d2cad5..21619d6 100644 --- a/docs/tutorials.rst +++ b/docs/tutorials.rst @@ -324,18 +324,22 @@ to 10 seconds (to avoid the large zero-lag pulse): .. note:: - Warning!! This command is particularly slow, especially for large data sets. The speed can - be increase by using ``--use-numba`` argument. + Warning!! This command is particularly slow, especially for large data sets. .. code-block:: $ rfpy_harmonics --no-outlier --find-azim --trange=2.,10. MMPY.pkl -To increase the processing speed use ``--use_numba`` to utilize ``numba`` just-in-time compiler: +To increase the processing speed use ``--use-numba`` to utilize ``numba`` just-in-time compiler: + +.. note:: + + Warning!! While this option will significantly increase the performance, due to the current behavior of numba + generic CTRL+C command will not stop the process. .. code-block:: - $ rfpy_harmonics --use_numba --no-outlier --find-azim --trange=2.,10. MMPY.pkl + $ rfpy_harmonics --use-numba --no-outlier --find-azim --trange=2.,10. MMPY.pkl ################################################################################ # __ _ _ #