From dbaead4754c0fade758aadf30db3e393e320adc0 Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Wed, 18 May 2016 12:00:45 +0200 Subject: [PATCH 01/15] code cleanup + commenting --- pylot/core/active/surveyUtils.py | 72 +++++++++----------------------- 1 file changed, 19 insertions(+), 53 deletions(-) diff --git a/pylot/core/active/surveyUtils.py b/pylot/core/active/surveyUtils.py index ff637612..04c3509e 100644 --- a/pylot/core/active/surveyUtils.py +++ b/pylot/core/active/surveyUtils.py @@ -16,26 +16,13 @@ def readParameters(parfile, parameter): return value - -def setArtificialPick(shot_dict, traceID, pick): - """ - - :param shot_dict: - :param traceID: - :param pick: - :return: - """ - for shot in shot_dict.values(): - shot.setPick(traceID, pick) - shot.setPickwindow(traceID, shot.getCut()) - - def fitSNR4dist(shot_dict, shiftdist=30, shiftSNR=100): """ + Approach to fit the decreasing SNR with wave travel distance. - :param shot_dict: - :param shiftdist: - :param shiftSNR: + :param shot_dict: dictionary containing Seismicshot objects (e.g. survey.getShotDict()) + :param shiftdist: shift compensating curve by a certain distance to the left + :param shiftSNR: shift compensating curve by a certain SNR value to the bottom :return: """ import numpy as np @@ -87,6 +74,8 @@ def plotFittedSNR(dists, snrthresholds, snrs, snrBestFit): def setDynamicFittedSNR(shot_dict, shiftdist=30, shiftSNR=100, p1=0.004, p2=-0.0007): """ + Set SNR values for a dictionary containing Seismicshots (e.g. survey.getShotDict()) + by parameters calulated from fitSNR4dist. :param shot_dict: :type shot_dict: dict @@ -119,6 +108,7 @@ def setDynamicFittedSNR(shot_dict, shiftdist=30, shiftSNR=100, p1=0.004, p2=-0.0 def setConstantSNR(shot_dict, snrthreshold=2.5): """ + Set a constant SNR value to all Seismicshots in a dictionary (e.g. survey.getShotDict()). :param shot_dict: :param snrthreshold: @@ -158,42 +148,18 @@ def findTracesInRanges(shot_dict, distancebin, pickbin): def cleanUp(survey): """ - - :param survey: - :return: + Cleans up a Survey object by removing frontend information + which can not be saved in the pickle format. """ for shot in survey.data.values(): shot.traces4plot = {} - -# def plotScatterStats(survey, key, ax = None): -# import matplotlib.pyplot as plt -# x = []; y = []; value = [] -# stats = survey.getStats() -# for shotnumber in stats.keys(): -# if type(value) == list: -# value.append(stats[shotnumber][key][0]) -# else: -# value.append(stats[shotnumber][key]) -# x.append(survey.data[shotnumber].getSrcLoc()[0]) -# y.append(survey.data[shotnumber].getSrcLoc()[1]) - -# if ax == None: -# fig = plt.figure() -# ax = fig.add_subplot(111) - -# sc = ax.scatter(x, y, s = value, c = value) -# plt.xlabel('X') -# plt.ylabel('Y') -# cbar = plt.colorbar(sc) -# cbar.set_label(key) - -def plotScatterStats4Shots(survey, key): +def plotScatterStats4Shots(survey, variable): """ Statistics, scatter plot. - key can be 'mean SNR', 'median SNR', 'mean SPE', 'median SPE', or 'picked traces' + :param survey: - :param key: + :param variable: can be 'mean SNR', 'median SNR', 'mean SPE', 'median SPE', or 'picked traces' :return: """ import matplotlib.pyplot as plt @@ -225,7 +191,7 @@ def plotScatterStats4Shots(survey, key): for shot in statsShot.keys(): x.append(statsShot[shot]['x']) y.append(statsShot[shot]['y']) - value.append(statsShot[shot][key]) + value.append(statsShot[shot][variable]) fig = plt.figure() ax = fig.add_subplot(111) @@ -239,19 +205,19 @@ def plotScatterStats4Shots(survey, key): plt.xlabel('X') plt.ylabel('Y') cbar = plt.colorbar(sc) - cbar.set_label(key) + cbar.set_label(variable) for shot in statsShot.keys(): ax.annotate(' %s' % shot.getShotnumber(), xy=(shot.getSrcLoc()[0], shot.getSrcLoc()[1]), fontsize='x-small', color='k') -def plotScatterStats4Receivers(survey, key): +def plotScatterStats4Receivers(survey, variable): """ Statistics, scatter plot. - key can be 'mean SNR', 'median SNR', 'mean SPE', 'median SPE', or 'picked traces' + :param survey: - :param key: + :param variable: can be 'mean SNR', 'median SNR', 'mean SPE', 'median SPE', or 'picked traces' :return: """ import matplotlib.pyplot as plt @@ -283,7 +249,7 @@ def plotScatterStats4Receivers(survey, key): for traceID in statsRec.keys(): x.append(statsRec[traceID]['x']) y.append(statsRec[traceID]['y']) - value.append(statsRec[traceID][key]) + value.append(statsRec[traceID][variable]) fig = plt.figure() ax = fig.add_subplot(111) @@ -297,7 +263,7 @@ def plotScatterStats4Receivers(survey, key): plt.xlabel('X') plt.ylabel('Y') cbar = plt.colorbar(sc) - cbar.set_label(key) + cbar.set_label(variable) shot = survey.data.values()[0] for traceID in shot.getTraceIDlist(): From 3138bbfa9381d268853c9db704b203ead7ee5506 Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Thu, 19 May 2016 11:20:37 +0200 Subject: [PATCH 02/15] code cleanup and commenting --- pylot/core/active/fmtomoUtils.py | 170 ++++++++++++++++++++++++++----- 1 file changed, 146 insertions(+), 24 deletions(-) diff --git a/pylot/core/active/fmtomoUtils.py b/pylot/core/active/fmtomoUtils.py index f70b0b43..086732d2 100644 --- a/pylot/core/active/fmtomoUtils.py +++ b/pylot/core/active/fmtomoUtils.py @@ -6,23 +6,15 @@ import datetime import numpy as np class Tomo3d(object): - def __init__(self, nproc, iterations, citer = 0, overwrite = False): + def __init__(self, citer = 0, overwrite = False): ''' Class build from FMTOMO script tomo3d. Can be used to run several instances of FMM code in parallel. - - :param: nproc, number of parallel processes - :type: integer - - :param: iterations, number of iterations - :type: integer :param: citer, current iteration (default = 0: start new model) :type: integer ''' self.setCWD() self.defParas() - self.nproc = nproc - self.iter = iterations # number of iterations self.citer = citer # current iteration self.sources = self.readSrcFile() self.traces = self.readTraces() @@ -34,21 +26,37 @@ class Tomo3d(object): self.defInvParas() def defFMMParas(self): + ''' + Initiates parameters for the forward calculation. + ''' + # Name of fast marching program self.fmm = '{0}/fm3d'.format(self.cwd) + # Name of program calculating Frechet derivatives self.frechgen = '{0}/frechgen'.format(self.cwd) + # Name of current velocity/inversion grid self.cvg = 'vgrids.in' + # Name of current interfaces grid self.cig = 'interfaces.in' + # Name of file containing current source locations self.csl = 'sources.in' + # Name of file containing propagation grid self.pg = 'propgrid.in' + # Name of file containing receiver coordinates self.rec = 'receivers.in' self.frech = 'frechet.in' self.frechout = 'frechet.dat' + # Name of file containing measured data self.ot = 'otimes.dat' + # Name of file containing output velocity information self.ttim = 'arrivals.dat' self.mode = 'mode_set.in' + # Name of temporary folders created for each process self.folder = '.proc_' def defInvParas(self): + ''' + Initiates inversion parameters for FMTOMO. + ''' # Name of program for performing inversion self.inv = '{0}/invert3d'.format(self.cwd) # Name of file containing current model traveltimes @@ -67,19 +75,41 @@ class Tomo3d(object): self.resout = 'residuals.dat' def copyRef(self): - # Copies reference grids to used grids (e.g. sourcesref.in to sources.in) + ''' + Copies reference grids to used grids (e.g. sourcesref.in to sources.in) + ''' os.system('cp %s %s'%(self.ivg, self.cvg)) os.system('cp %s %s'%(self.iig, self.cig)) os.system('cp %s %s'%(self.isl, self.csl)) - def setCWD(self): - self.cwd = subprocess.check_output(['pwd'])[0:-1] - print('Working directory is pwd: %s'%self.cwd) + def setCWD(self, directory = None): + ''' + Set working directory containing all necessary files. + + Default: pwd + ''' + if directory == None: + directory = subprocess.check_output(['pwd'])[0:-1] + + self.cwd = directory + print('Working directory is: %s'%self.cwd) def runFrech(self): os.system(self.frechgen) - def runTOMO3D(self): + def runTOMO3D(self, nproc, iterations): + ''' + Starts up the FMTOMO code for the set number of iterations on nproc parallel processes. + + :param: nproc, number of parallel processes + :type: integer + + :param: iterations, number of iterations + :type: integer + ''' + self.nproc = nproc + self.iter = iterations # number of iterations + starttime = datetime.datetime.now() print('Starting TOMO3D on %s parallel processes for %s iteration(s).' %(self.nproc, self.iter)) @@ -104,6 +134,10 @@ class Tomo3d(object): print('runTOMO3D: Finished %s iterations after %s.'%(self.iter, tdelta)) def runFmm(self, directory, logfile, processes): + ''' + Calls an instance of the FMM code in the process directory. + Requires a list of all active processes and returns an updated list. + ''' os.chdir(directory) fout = open(logfile, 'w') processes.append(subprocess.Popen(self.fmm, stdout = fout)) @@ -112,6 +146,9 @@ class Tomo3d(object): return processes def startForward(self, logdir): + ''' + Runs an instance of the FMM code in the process directory. + ''' self._printLine() print('Starting forward simulation for iteration %s.'%(self.citer)) @@ -128,8 +165,8 @@ class Tomo3d(object): logfn = 'fm3dlog_' + str(procID) + '.out' log_out = os.path.join(logdir, logfn) - self.writeSrcFile(procID, directory) - self.writeTracesFile(procID, directory) + self.writeSrcFile(procID) + self.writeTracesFile(procID) os.system('cp {cvg} {cig} {mode} {pg} {frechin} {dest}' .format(cvg=self.cvg, cig=self.cig, frechin=self.frech, mode=self.mode, pg=self.pg, dest=directory)) @@ -149,10 +186,16 @@ class Tomo3d(object): print('Finished Forward calculation after %s'%tdelta) def startInversion(self): + ''' + Simply calls the inversion program. + ''' print('Calling %s...'%self.inv) os.system(self.inv) def calcRes(self): + ''' + Calls residual calculation program. + ''' resout = os.path.join(self.cwd, self.resout) if self.citer == 0: os.system('%s > %s'%(self.resid, resout)) @@ -185,11 +228,17 @@ class Tomo3d(object): raise RuntimeError('Could not create directory: %s'%directory) def makeDirectories(self): + ''' + Makes temporary directories for all processes. + ''' for procID in range(1, self.nproc + 1): directory = self.getProcDir(procID) self.makeDir(directory) def makeInvIterDir(self): + ''' + Makes directories for each iteration step for the output. + ''' invIterDir = self.cwd + '/it_%s'%(self.citer) err = os.system('mkdir %s'%invIterDir) if err is 256: @@ -200,12 +249,18 @@ class Tomo3d(object): self.cInvIterDir = invIterDir def clearDir(self, directory): + ''' + Wipes a certain directory. + ''' print('Wiping directory %s...'%directory) for filename in os.listdir(directory): filenp = os.path.join(directory, filename) os.remove(filenp) def clearDirectories(self): + ''' + Wipes all generated temporary directories. + ''' for directory in self.directories: self.clearDir(directory) @@ -214,14 +269,24 @@ class Tomo3d(object): return os.rmdir(directory) def removeDirectories(self): + ''' + Removes all generated temporary directories. + ''' for directory in self.directories: self.rmDir(directory) self.directories = [] def getProcDir(self, procID): + ''' + Returns the temporary directory for a certain process + with procID = process number. + ''' return os.path.join(self.cwd, self.folder) + str(procID) def getTraceIDs4Sources(self, sourceIDs): + ''' + Returns corresponding trace IDs for a set of given source IDs. + ''' traceIDs = [] for traceID in self.traces.keys(): if self.traces[traceID]['source'] in sourceIDs: @@ -229,6 +294,9 @@ class Tomo3d(object): return traceIDs def getTraceIDs4Source(self, sourceID): + ''' + Returns corresponding trace IDs for a source ID. + ''' traceIDs = [] for traceID in self.traces.keys(): if self.traces[traceID]['source'] == sourceID: @@ -236,22 +304,38 @@ class Tomo3d(object): return traceIDs def copyArrivals(self, target = None): + ''' + Copies the FMM output file (self.ttim) to a specific target file. + Default target is self.mtrav (model travel times). + ''' if target == None: target = os.path.join(self.cwd, self.mtrav) os.system('cp %s %s'%(os.path.join( self.cInvIterDir, self.ttim), target)) def saveVgrid(self): - vgpath = os.path.join(self.cwd, 'vgrids.in') + ''' + Saves the current velocity grid for the current iteration step. + ''' + vgpath = os.path.join(self.cwd, self.cvg) os.system('cp %s %s'%(vgpath, self.cInvIterDir)) def calcSrcPerKernel(self): + ''' + (Equally) distributes all sources depending on the number of processes (kernels). + Returns two integer values. + First: minimum number of sources for each process + Second: remaining sources (always less than number of processes) + ''' nsrc = self.readNsrc() if self.nproc > nsrc: raise ValueError('Number of spawned processes higher than number of sources') return nsrc/self.nproc, nsrc%self.nproc def srcIDs4Kernel(self, procID): + ''' + Calculates and returns all source IDs for a given process ID. + ''' proc = procID - 1 nsrc = self.readNsrc() srcPK, remain = self.calcSrcPerKernel() @@ -274,12 +358,18 @@ class Tomo3d(object): return nsrc def readNtraces(self): + ''' + Reads the total number of traces from self.rec header. + ''' recfile = open(self.rec, 'r') nrec = int(recfile.readline()) recfile.close() return nrec def readSrcFile(self): + ''' + Reads the whole sourcefile and returns structured information in a dictionary. + ''' nsrc = self.readNsrc() srcfile = open(self.csl, 'r') @@ -309,6 +399,10 @@ class Tomo3d(object): return sources def readTraces(self): + ''' + Reads the receiver input file and returns the information + in a structured dictionary. + ''' recfile = open(self.rec, 'r') ntraces = self.readNtraces() @@ -330,8 +424,13 @@ class Tomo3d(object): return traces def readArrivals(self, procID): + ''' + Reads the arrival times from a temporary process directory, + changes local to global sourceIDs and traceIDs and returns + a list of arrival times. + ''' directory = self.getProcDir(procID) - arrfile = open(directory + '/arrivals.dat', 'r') + arrfile = open(os.path.join(directory, self.ttime), 'r') sourceIDs = self.srcIDs4Kernel(procID) arrivals = [] @@ -347,6 +446,10 @@ class Tomo3d(object): return arrivals def readRays(self, procID): + ''' + Reads rays output from a temporary process directory and returns + the information in a structured dictionary. + ''' directory = self.getProcDir(procID) raysfile = open(directory + '/rays.dat', 'r') sourceIDs = self.srcIDs4Kernel(procID) @@ -385,8 +488,12 @@ class Tomo3d(object): } return rays - def writeSrcFile(self, procID, directory): - srcfile = open('%s/sources.in'%directory, 'w') + def writeSrcFile(self, procID): + ''' + Writes a source input file for a process with ID = procID. + ''' + directory = self.getProcDir(procID) + srcfile = open(os.path,join(directory, self.csl), 'w') sourceIDs = self.srcIDs4Kernel(procID) srcfile.write('%s\n'%len(sourceIDs)) @@ -401,7 +508,11 @@ class Tomo3d(object): srcfile.write('%s %s\n'%(int(interactions[0]), int(interactions[1]))) srcfile.write('%s\n'%source['veltype']) - def writeTracesFile(self, procID, directory): + def writeTracesFile(self, procID): + ''' + Writes a receiver input file for a process with ID = procID. + ''' + directory = self.getProcDir(procID) recfile = open('%s/receivers.in'%directory, 'w') sourceIDs = self.srcIDs4Kernel(procID) traceIDs = self.getTraceIDs4Sources(sourceIDs) @@ -417,9 +528,12 @@ class Tomo3d(object): recfile.write('%s\n'%trace['path']) def mergeArrivals(self, directory): + ''' + Merges the arrival times for all processes to self.cInvIterDir. + ''' arrfn = os.path.join(directory, self.ttim) arrivalsOut = open(arrfn, 'w') - print('Merging arrivals.dat...') + print('Merging %s...'%self.ttim) for procID in range(1, self.nproc + 1): arrivals = self.readArrivals(procID) for line in arrivals: @@ -428,6 +542,9 @@ class Tomo3d(object): os.system('ln -fs %s %s'%(arrfn, os.path.join(self.cwd, self.ttim))) def mergeRays(self, directory): + ''' + Merges the ray paths for all processes to self.cInvIterDir. + ''' print('Merging rays.dat...') with open(directory + '/rays.dat', 'w') as outfile: for procID in range(1, self.nproc + 1): @@ -448,9 +565,11 @@ class Tomo3d(object): outfile.writelines(raysec['raypoints']) def mergeFrechet(self, directory): - print('Merging frechet.dat...') - + ''' + Merges the frechet derivatives for all processes to self.cInvIterDir. + ''' frechfnout = os.path.join(directory, self.frechout) + print('Merging %s...'%self.frechout) with open(frechfnout, 'w') as outfile: for procID in range(1, self.nproc + 1): filename = os.path.join(self.getProcDir(procID), self.frechout) @@ -465,6 +584,9 @@ class Tomo3d(object): os.system('ln -fs %s %s'%(frechfnout, os.path.join(self.cwd, self.frechout))) def mergeOutput(self, directory): + ''' + Calls self.mergeArrivals, self.mergeFrechet and self.mergeRays. + ''' self.mergeArrivals(directory) self.mergeFrechet(directory) self.mergeRays(directory) From 4edd5f8e52eb461bc28846bcb6f490aa10fda55d Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Thu, 19 May 2016 11:21:24 +0200 Subject: [PATCH 03/15] code cleanup --- pylot/core/active/seismicshot.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pylot/core/active/seismicshot.py b/pylot/core/active/seismicshot.py index 68087561..20b8fb38 100644 --- a/pylot/core/active/seismicshot.py +++ b/pylot/core/active/seismicshot.py @@ -283,7 +283,7 @@ class SeismicShot(object): # raise ValueError('ambigious or empty traceID: %s' % traceID) - def pickTraces(self, traceID, windowsize, folm, HosAic='hos'): ########## input variables ########## + def pickTraces(self, traceID, folm, HosAic='hos', windowsize = (10, 0)): ########## input variables ########## # LOCALMAX NOT IMPLEMENTED! ''' Intitiate picking for a trace. From db17cb4f8d96924ef1b8e4dfbb8d3884f02e451c Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Thu, 19 May 2016 14:24:48 +0200 Subject: [PATCH 04/15] code cleanup --- pylot/core/active/seismicshot.py | 66 ++------------------------------ 1 file changed, 3 insertions(+), 63 deletions(-) diff --git a/pylot/core/active/seismicshot.py b/pylot/core/active/seismicshot.py index 20b8fb38..a5c0d26b 100644 --- a/pylot/core/active/seismicshot.py +++ b/pylot/core/active/seismicshot.py @@ -64,11 +64,6 @@ class SeismicShot(object): if traceID == trace.stats.channel: self.traces.remove(trace) - # for traceID in TraceIDs: - # traces = [trace for trace in self.traces if int(trace.stats.channel) == traceID] - # if len(traces) is not 1: - # self.traces.remove(trace) - def updateTraceList(self): ''' Looks for empty traces, returns a list of deleted traceIDs. @@ -151,7 +146,6 @@ class SeismicShot(object): if not self.getPickFlag(traceID) == 0: return self.picks[traceID]['mpp'] if returnRemoved == True: - # print('getPick: Returned removed pick for shot %d, traceID %d' %(self.getShotnumber(), traceID)) return self.picks[traceID]['mpp'] def getPickIncludeRemoved(self, traceID): @@ -281,10 +275,7 @@ class SeismicShot(object): self.setPick(traceID, None) print 'Warning: ambigious or empty traceID: %s' % traceID - # raise ValueError('ambigious or empty traceID: %s' % traceID) - def pickTraces(self, traceID, folm, HosAic='hos', windowsize = (10, 0)): ########## input variables ########## - # LOCALMAX NOT IMPLEMENTED! ''' Intitiate picking for a trace. @@ -335,11 +326,6 @@ class SeismicShot(object): if self.picks[traceID]['epp'] < 0: self.picks[traceID]['epp'] - # print('setEarllatepick: Set epp to 0 because it was < 0') - - # TEST OF 1/2 PICKERROR - # self.picks[traceID]['spe'] *= 0.5 - # TEST OF 1/2 PICKERROR def threshold(self, hoscf, aiccf, windowsize, pickwindow, folm): ''' @@ -368,12 +354,8 @@ class SeismicShot(object): leftb = int(pickwindow[0] / self.getCut()[1] * len(hoscflist)) rightb = int(pickwindow[1] / self.getCut()[1] * len(hoscflist)) - # threshold = folm * max(hoscflist[leftb : rightb]) # combination of local maximum and threshold - - ### TEST TEST threshold = folm * (max(hoscflist[leftb: rightb]) - min(hoscflist[leftb: rightb])) + min( hoscflist[leftb: rightb]) # combination of local maximum and threshold - ### TEST TEST m = leftb @@ -411,7 +393,6 @@ class SeismicShot(object): raise ValueError("Distance is NaN for traceID %s" % traceID) return dist - # return abs(float(self.getSrcLoc(traceID))-float(self.getRecLoc(traceID))) def getRecLoc(self, traceID): ########## input FILENAME ########## ''' @@ -432,9 +413,7 @@ class SeismicShot(object): z = coordlist[i].split()[3] return float(x), float(y), float(z) - # print "WARNING: traceID %s not found" % traceID raise ValueError("traceID %s not found" % traceID) - # return float(self.getSingleStream(traceID)[0].stats.seg2['RECEIVER_LOCATION']) def getSrcLoc(self): ########## input FILENAME ########## ''' @@ -477,26 +456,6 @@ class SeismicShot(object): if len(traceID_list) > 0: return traceID_list - # def setManualPicks(self, traceID, picklist): ########## picklist momentan nicht allgemein, nur testweise benutzt ########## - # ''' - # Sets the manual picks for a receiver with the ID == traceID for comparison. - - # :param: traceID - # :type: int - - # :param: picklist, list containing the manual picks (mostlikely, earliest, latest). - # :type: list - # ''' - # picks = picklist[traceID - 1].split() - # mostlikely = float(picks[1]) - # earliest = float(picks[2]) - # latest = float(picks[3]) - - # if not self.manualpicks.has_key(traceID): - # self.manualpicks[traceID] = (mostlikely, earliest, latest) - # else: - # raise KeyError('MANUAL pick to be set more than once for traceID %s' % traceID) - def setManualPicksFromFile(self, directory='picks'): ''' Read manual picks from *.pck file. @@ -565,7 +524,6 @@ class SeismicShot(object): :param: (tnoise, tgap, tsignal), as used in pylot SNR ''' - from pylot.core.pick.utils import getSNR tgap = self.getTgap() @@ -591,7 +549,7 @@ class SeismicShot(object): return distancearray - def plot2dttc(self, ax=None): ########## 2D ########## + def plot2dttc(self, ax=None): ''' Function to plot the traveltime curve for automated picks of a shot. 2d only! ATM: X DIRECTION!! ''' @@ -617,7 +575,7 @@ class SeismicShot(object): ax.text(0.5, 0.9, 'shot: %s' % self.getShotnumber(), transform=ax.transAxes , horizontalalignment='center') - def plotmanual2dttc(self, ax=None): ########## 2D ########## + def plotmanual2dttc(self, ax=None): ''' Function to plot the traveltime curve for manual picks of a shot. 2D only! ''' @@ -643,24 +601,6 @@ class SeismicShot(object): y.append(point[1]) ax.plot(x, y, 'b', label="Manual Picks") - # def plotpickwindow(self): ########## 2D ########## - # ''' - # Plots the pickwindow of a shot for the 2nd iteration step of picking. 2D only! - # ''' - # import matplotlib.pyplot as plt - # plt.interactive('True') - # pickwindowarray_lowerb = [] - # pickwindowarray_upperb = [] - - # i = 1 - # for traceID in self.pwindow.keys(): - # pickwindowarray_lowerb.append(self.pwindow[traceID][0]) - # pickwindowarray_upperb.append(self.pwindow[traceID][1]) - # i += 1 - - # plt.plot(self.getDistArray4ttcPlot(), pickwindowarray_lowerb, ':k') - # plt.plot(self.getDistArray4ttcPlot(), pickwindowarray_upperb, ':k') - def plotTrace(self, traceID, plotSNR=True, lw=1): fig = plt.figure() ax = fig.add_subplot(111) @@ -679,7 +619,7 @@ class SeismicShot(object): ax.legend() ax.text(0.05, 0.9, 'SNR: %s' % snr, transform=ax.transAxes) - def plot_traces(self, traceID): ########## 2D, muss noch mehr verbessert werden ########## + def plot_traces(self, traceID): from matplotlib.widgets import Button def onclick(event): From 1f47f3dd850114b2ecd3aec6f7a7e098ba57757f Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Mon, 23 May 2016 11:22:39 +0200 Subject: [PATCH 05/15] parallization of picking algorithm --- pylot/core/active/seismicshot.py | 73 +++++++++++++++++++++++++++++--- 1 file changed, 68 insertions(+), 5 deletions(-) diff --git a/pylot/core/active/seismicshot.py b/pylot/core/active/seismicshot.py index a5c0d26b..827b12ac 100644 --- a/pylot/core/active/seismicshot.py +++ b/pylot/core/active/seismicshot.py @@ -78,6 +78,12 @@ class SeismicShot(object): def setParameters(self, name, value): self.paras[name] = value + def setVmin(self, vmin): + self.setParameters('vmin', vmin) + + def setVmax(self, vmax): + self.setParameters('vmax', vmax) + def setCut(self, cut): self.setParameters('cut', cut) @@ -102,6 +108,43 @@ class SeismicShot(object): def setSourcefile(self, sourcefile): self.setParameters('sourcefile', sourcefile) + def setMethod(self, method): + self.setParameters('method', method) + + def setAicwindow(self, aicwindow): + self.setParameters('aicwindow', aicwindow) + + def setFolm(self, folm): + self.setParameters('folm', folm) + + def setDynPickwindow(self, traceID, cutdist = 5.): + distance = self.getDistance(traceID) # receive distance + + vmin = self.getVmin() + vmax = self.getVmax() + + pickwin_used = self.getCut() + cutwindow = self.getCut() + + # for higher distances use a linear vmin/vmax to cut out late/early regions with high noise + if distance > cutdist: + pwleft = distance / vmax + pwright = distance / vmin + if pwright > cutwindow[1]: + pwright = cutwindow[1] + pickwin_used = (pwleft, pwright) + + self.setPickwindow(traceID, pickwin_used) + + def getMethod(self): + return self.paras['method'] + + def getAicwindow(self): + return self.paras['aicwindow'] + + def getFolm(self): + return self.paras['folm'] + def getParas(self): return self.paras @@ -132,6 +175,12 @@ class SeismicShot(object): def getSourcefile(self): return self.paras['sourcefile'] + def getVmin(self): + return self.paras['vmin'] + + def getVmax(self): + return self.paras['vmax'] + def getManualPick(self, traceID): if not self.getManualPickFlag(traceID) == 0: return self.manualpicks[traceID]['mpp'] @@ -275,7 +324,21 @@ class SeismicShot(object): self.setPick(traceID, None) print 'Warning: ambigious or empty traceID: %s' % traceID - def pickTraces(self, traceID, folm, HosAic='hos', windowsize = (10, 0)): ########## input variables ########## + def pickParallel(self, folm, method = 'hos', aicwindow = (10, 0)): + import multiprocessing + + self.setFolm(folm) + self.setMethod(method) + self.setAicwindow(aicwindow) + + maxthreads = multiprocessing.cpu_count() + pool = multiprocessing.Pool(maxthreads) + + traceIDs = self.getTraceIDlist() + + pool.map(self.pickTrace, traceIDs) + + def pickTrace(self, traceID): ''' Intitiate picking for a trace. @@ -300,17 +363,17 @@ class SeismicShot(object): :param: HosAic, get hos or aic pick (can be 'hos'(default) or 'aic') :type: 'string' ''' + self.setDynPickwindow(traceID) + hoscf = self.getHOScf(traceID) ### determination of both, HOS and AIC (need to change threshold-picker) ### aiccf = self.getAICcf(traceID) - self.folm = folm - self.timeArray[traceID] = hoscf.getTimeArray() - aiccftime, hoscftime = self.threshold(hoscf, aiccf, windowsize, self.getPickwindow(traceID), folm) + aiccftime, hoscftime = self.threshold(hoscf, aiccf, self.getAicwindow(), self.getPickwindow(traceID), self.getFolm()) setHosAic = {'hos': hoscftime, 'aic': aiccftime} - self.setPick(traceID, setHosAic[HosAic]) + self.setPick(traceID, setHosAic[self.getMethod()]) def setEarllatepick(self, traceID, nfac=1.5): tgap = self.getTgap() From 73d71a61d5058edaf970bc26520b2ffa08b44d1a Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Mon, 23 May 2016 11:23:23 +0200 Subject: [PATCH 06/15] restructuring for parallization --- pylot/core/active/activeSeismoPick.py | 55 +++++++++++++-------------- 1 file changed, 26 insertions(+), 29 deletions(-) diff --git a/pylot/core/active/activeSeismoPick.py b/pylot/core/active/activeSeismoPick.py index 2d02496d..44381d5b 100644 --- a/pylot/core/active/activeSeismoPick.py +++ b/pylot/core/active/activeSeismoPick.py @@ -216,45 +216,42 @@ class Survey(object): for shot in self.data.values(): tstartpick = datetime.now(); + shot.setVmin(vmin) + shot.setVmax(vmax) count += 1 - for traceID in shot.getTraceIDlist(): - distance = shot.getDistance(traceID) # receive distance + shot.pickParallel(folm) - pickwin_used = shot.getCut() - cutwindow = shot.getCut() + # tpicksum += (datetime.now() - tstartpick); + # tpick = tpicksum / count + # tremain = (tpick * (len(self.getShotDict()) - count)) + # tend = datetime.now() + tremain + # progress = float(count) / float(len(self.getShotDict())) * 100 + # self._update_progress(shot.getShotname(), tend, progress) - # for higher distances use a linear vmin/vmax to cut out late/early regions with high noise - if distance > 5.: - pwleft = distance / vmax ################## TEST - pwright = distance / vmin - if pwright > cutwindow[1]: - pwright = cutwindow[1] - pickwin_used = (pwleft, pwright) + self.setSNR() + self.setEarllate() - shot.setPickwindow(traceID, pickwin_used) - shot.pickTraces(traceID, folm, HosAic, aicwindow) # picker - - shot.setSNR(traceID) - # if shot.getSNR(traceID)[0] < snrthreshold: - if shot.getSNR(traceID)[0] < shot.getSNRthreshold(traceID): - shot.removePick(traceID) - - # set epp and lpp if SNR > 1 (else earllatepicker cant set values) - if shot.getSNR(traceID)[0] > 1: - shot.setEarllatepick(traceID) - - tpicksum += (datetime.now() - tstartpick); - tpick = tpicksum / count - tremain = (tpick * (len(self.getShotDict()) - count)) - tend = datetime.now() + tremain - progress = float(count) / float(len(self.getShotDict())) * 100 - self._update_progress(shot.getShotname(), tend, progress) print('\npickAllShots: Finished\n') ntraces = self.countAllTraces() pickedtraces = self.countAllPickedTraces() print('Picked %s / %s traces (%d %%)\n' % (pickedtraces, ntraces, float(pickedtraces) / float(ntraces) * 100.)) + def setSNR(self): + for shot in self.data.values(): + for traceID in shot.getTraceIDlist(): + shot.setSNR(traceID) + # if shot.getSNR(traceID)[0] < snrthreshold: + if shot.getSNR(traceID)[0] < shot.getSNRthreshold(traceID): + shot.removePick(traceID) + + def setEarllate(self): + for shot in self.data.values(): + for traceID in shot.getTraceIDlist(): + # set epp and lpp if SNR > 1 (else earllatepicker cant set values) + if shot.getSNR(traceID)[0] > 1: + shot.setEarllatepick(traceID) + def cleanBySPE(self, maxSPE): ''' Sets all picks as invalid if they exceed a certain value of the symmetric pick error. From 3cc77f4868d9284e79814b557cf7dcaae5d51191 Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Mon, 23 May 2016 11:24:01 +0200 Subject: [PATCH 07/15] bugfixes --- pylot/core/active/fmtomoUtils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pylot/core/active/fmtomoUtils.py b/pylot/core/active/fmtomoUtils.py index 086732d2..6db9705d 100644 --- a/pylot/core/active/fmtomoUtils.py +++ b/pylot/core/active/fmtomoUtils.py @@ -430,7 +430,7 @@ class Tomo3d(object): a list of arrival times. ''' directory = self.getProcDir(procID) - arrfile = open(os.path.join(directory, self.ttime), 'r') + arrfile = open(os.path.join(directory, self.ttim), 'r') sourceIDs = self.srcIDs4Kernel(procID) arrivals = [] @@ -493,7 +493,7 @@ class Tomo3d(object): Writes a source input file for a process with ID = procID. ''' directory = self.getProcDir(procID) - srcfile = open(os.path,join(directory, self.csl), 'w') + srcfile = open(os.path.join(directory, self.csl), 'w') sourceIDs = self.srcIDs4Kernel(procID) srcfile.write('%s\n'%len(sourceIDs)) From 41b7ca6968a77bca06d8ca386277dabaeec0a6ba Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Mon, 23 May 2016 11:53:22 +0200 Subject: [PATCH 08/15] [task] reformatting activeSeismoPick and editing pool mapping to work properly --- pylot/core/active/activeSeismoPick.py | 140 +++++++++++++++----------- pylot/core/active/seismicshot.py | 11 +- 2 files changed, 92 insertions(+), 59 deletions(-) diff --git a/pylot/core/active/activeSeismoPick.py b/pylot/core/active/activeSeismoPick.py index 44381d5b..f2c18d83 100644 --- a/pylot/core/active/activeSeismoPick.py +++ b/pylot/core/active/activeSeismoPick.py @@ -4,6 +4,7 @@ import numpy as np from pylot.core.active import seismicshot from pylot.core.active.surveyUtils import cleanUp + class Survey(object): def __init__(self, path, sourcefile, receiverfile, useDefaultParas=False): ''' @@ -21,13 +22,13 @@ class Survey(object): self._sourcefile = sourcefile self._obsdir = path self._generateSurvey() - self._initiateFilenames() + self._initiate_fnames() if useDefaultParas == True: self.setParametersForAllShots() self._removeAllEmptyTraces() self._updateShots() - def _initiateFilenames(self): + def _initiate_fnames(self): for shot in self.data.values(): shot.setRecfile(self.getPath() + self.getReceiverfile()) shot.setSourcefile(self.getPath() + self.getSourcefile()) @@ -74,7 +75,7 @@ class Survey(object): but were set in the input files. ''' logfile = 'updateShots.out' - count = 0; + count = 0 countTraces = 0 for shot in self.data.values(): del_traceIDs = shot.updateTraceList() @@ -93,15 +94,8 @@ class Survey(object): "on removed traces." % (logfile)) outfile.close() - def setArtificialPick(self, traceID, pick): - ''' - Sets an artificial pick for a traceID of all shots in the survey object. - (Commonly used to generate a travel time t = 0 at the source origin) - ''' - for shot in self.data.values(): - shot.setPick(traceID, pick) - - def setParametersForAllShots(self, cutwindow=(0, 0.2), tmovwind=0.3, tsignal=0.03, tgap=0.0007): + def setParametersForAllShots(self, cutwindow=(0, 0.2), tmovwind=0.3, + tsignal=0.03, tgap=0.0007): if (cutwindow == (0, 0.2) and tmovwind == 0.3 and tsignal == 0.03 and tgap == 0.0007): print ("Warning: Standard values used for " @@ -136,8 +130,10 @@ class Survey(object): if not shot in diffs.keys(): diffs[shot] = {} for traceID in shot.getTraceIDlist(): - if shot.getPickFlag(traceID) == 1 and shot.getManualPickFlag(traceID) == 1: - diffs[shot][traceID] = shot.getPick(traceID) - shot.getManualPick(traceID) + if shot.getPickFlag(traceID) == 1 and shot.getManualPickFlag( + traceID) == 1: + diffs[shot][traceID] = shot.getPick( + traceID) - shot.getManualPick(traceID) return diffs def plotDiffs(self): @@ -146,14 +142,15 @@ class Survey(object): difference between automatic and manual pick. ''' import matplotlib.pyplot as plt - diffs = []; - dists = []; - mpicks = []; + diffs = [] + dists = [] + mpicks = [] picks = [] diffsDic = self.getDiffsFromManual() for shot in self.data.values(): for traceID in shot.getTraceIDlist(): - if shot.getPickFlag(traceID) == 1 and shot.getManualPickFlag(traceID) == 1: + if shot.getPickFlag(traceID) == 1 and shot.getManualPickFlag( + traceID) == 1: dists.append(shot.getDistance(traceID)) mpicks.append(shot.getManualPick(traceID)) picks.append(shot.getPick(traceID)) @@ -165,13 +162,16 @@ class Survey(object): fig = plt.figure() ax = fig.add_subplot(111) - sc_a = ax.scatter(dists, picks, c='0.5', s=10, edgecolors='none', label=labela, alpha=0.3) - sc = ax.scatter(dists, mpicks, c=diffs, s=5, edgecolors='none', label=labelm) + sc_a = ax.scatter(dists, picks, c='0.5', s=10, edgecolors='none', + label=labela, alpha=0.3) + sc = ax.scatter(dists, mpicks, c=diffs, s=5, edgecolors='none', + label=labelm) cbar = plt.colorbar(sc, fraction=0.05) cbar.set_label(labelm) ax.set_xlabel('Distance [m]') ax.set_ylabel('Time [s]') - ax.text(0.5, 0.95, 'Plot of all MANUAL picks', transform=ax.transAxes, horizontalalignment='center') + ax.text(0.5, 0.95, 'Plot of all MANUAL picks', transform=ax.transAxes, + horizontalalignment='center') plt.legend() def plotHist(self, nbins=20, ax=None): @@ -186,14 +186,18 @@ class Survey(object): ax = fig.add_subplot(111) for shot in self.data.values(): for traceID in shot.getTraceIDlist(): - if shot.getPickFlag(traceID) == 1 and shot.getManualPickFlag(traceID) == 1: + if shot.getPickFlag(traceID) == 1 and shot.getManualPickFlag( + traceID) == 1: diffs.append(self.getDiffsFromManual()[shot][traceID]) - hist = plt.hist(diffs, nbins, histtype='step', normed=True, stacked=True) - plt.title('Histogram of the differences between automatic and manual pick') + hist = plt.hist(diffs, nbins, histtype='step', normed=True, + stacked=True) + plt.title( + 'Histogram of the differences between automatic and manual pick') plt.xlabel('Difference in time (auto - manual) [s]') return diffs - def pickAllShots(self, vmin=333, vmax=5500, folm=0.6, HosAic='hos', aicwindow=(10, 0)): + def pickAllShots(self, vmin=333, vmax=5500, folm=0.6, HosAic='hos', + aicwindow=(10, 0)): ''' Automatically pick all traces of all shots of the survey. @@ -211,11 +215,11 @@ class Survey(object): ''' from datetime import datetime starttime = datetime.now() - count = 0; + count = 0 tpicksum = starttime - starttime for shot in self.data.values(): - tstartpick = datetime.now(); + tstartpick = datetime.now() shot.setVmin(vmin) shot.setVmax(vmax) count += 1 @@ -235,7 +239,8 @@ class Survey(object): ntraces = self.countAllTraces() pickedtraces = self.countAllPickedTraces() print('Picked %s / %s traces (%d %%)\n' - % (pickedtraces, ntraces, float(pickedtraces) / float(ntraces) * 100.)) + % (pickedtraces, ntraces, + float(pickedtraces) / float(ntraces) * 100.)) def setSNR(self): for shot in self.data.values(): @@ -373,8 +378,11 @@ class Survey(object): pickedTraces += 1 info_dict[shot.getShotnumber()] = {'numtraces': numtraces, 'picked traces': [pickedTraces, - '%2.2f %%' % (float(pickedTraces) / - float(numtraces) * 100)], + '%2.2f %%' % ( + float( + pickedTraces) / + float( + numtraces) * 100)], 'mean SNR': np.mean(snrlist), 'mean distance': np.mean(dist)} @@ -388,10 +396,12 @@ class Survey(object): if shot.getShotnumber() == shotnumber: return shot - def exportFMTOMO(self, directory='FMTOMO_export', sourcefile='input_sf.in', ttFileExtension='.tt'): + def exportFMTOMO(self, directory='FMTOMO_export', sourcefile='input_sf.in', + ttFileExtension='.tt'): ''' Exports all picks into a directory as travel time files readable by FMTOMO obsdata. ''' + def getAngle(distance): PI = np.pi R = 6371. @@ -400,18 +410,20 @@ class Survey(object): count = 0 fmtomo_factor = 1000 # transforming [m/s] -> [km/s] - LatAll = []; - LonAll = []; + LatAll = [] + LonAll = [] DepthAll = [] srcfile = open(directory + '/' + sourcefile, 'w') srcfile.writelines('%10s\n' % len(self.data)) # number of sources for shotnumber in self.getShotlist(): shot = self.getShotForShotnumber(shotnumber) - ttfilename = str(shotnumber) + ttFileExtension # filename of travel time file for this shot + ttfilename = str( + shotnumber) + ttFileExtension # filename of travel time file for this shot (x, y, z) = shot.getSrcLoc() # getSrcLoc returns (x, y, z) - srcfile.writelines('%10s %10s %10s\n' % (getAngle(y), getAngle(x), (-1) * z)) # transform to lat, lon, depth - LatAll.append(getAngle(y)); - LonAll.append(getAngle(x)); + srcfile.writelines('%10s %10s %10s\n' % ( + getAngle(y), getAngle(x), (-1) * z)) # transform to lat, lon, depth + LatAll.append(getAngle(y)) + LonAll.append(getAngle(x)) DepthAll.append((-1) * z) srcfile.writelines('%10s\n' % 1) srcfile.writelines('%10s %10s %10s\n' % (1, 1, ttfilename)) @@ -424,9 +436,10 @@ class Survey(object): pick = shot.getPick(traceID) * fmtomo_factor delta = shot.getSymmetricPickError(traceID) * fmtomo_factor (x, y, z) = shot.getRecLoc(traceID) - ttfile.writelines('%20s %20s %20s %10s %10s\n' % (getAngle(y), getAngle(x), (-1) * z, pick, delta)) - LatAll.append(getAngle(y)); - LonAll.append(getAngle(x)); + ttfile.writelines('%20s %20s %20s %10s %10s\n' % ( + getAngle(y), getAngle(x), (-1) * z, pick, delta)) + LatAll.append(getAngle(y)) + LonAll.append(getAngle(x)) DepthAll.append((-1) * z) count += 1 ttfile.close() @@ -486,19 +499,24 @@ class Survey(object): if index <= figPerSubplot: ax = fig.add_subplot(rows, columns, index) if mode == '3d': - self.getShot(shotnumber).matshow(ax=ax, colorbar=False, annotations=True, legend=False) + self.getShot(shotnumber).matshow(ax=ax, colorbar=False, + annotations=True, + legend=False) elif mode == '2d': self.getShot(shotnumber).plot2dttc(ax) self.getShot(shotnumber).plotmanual2dttc(ax) index += 1 if index > figPerSubplot: - fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0) + fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, + hspace=0) fig = plt.figure() index = 1 - fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0) + fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, + hspace=0) - def plotAllPicks(self, plotRemoved=False, colorByVal='log10SNR', ax=None, cbar=None, refreshPlot=False): + def plotAllPicks(self, plotRemoved=False, colorByVal='log10SNR', ax=None, + cbar=None, refreshPlot=False): ''' Plots all picks over the distance between source and receiver. Returns (ax, region). Picks can be checked and removed by using region class (pylot.core.active.surveyPlotTools.regions) @@ -542,7 +560,8 @@ class Survey(object): for shot in self.data.values(): for traceID in shot.getTraceIDlist(): if plotRemoved == False: - if shot.getPickFlag(traceID) is not 0 or plotRemoved == True: + if shot.getPickFlag( + traceID) is not 0 or plotRemoved == True: dist.append(shot.getDistance(traceID)) pick.append(shot.getPick(traceID)) snrlog.append(math.log10(shot.getSNR(traceID)[0])) @@ -554,12 +573,15 @@ class Survey(object): 'spe': spe} self.color = color if refreshPlot is False: - ax, cbar = self.createPlot(dist, pick, color[colorByVal], label='%s' % colorByVal) + ax, cbar = self.createPlot(dist, pick, color[colorByVal], + label='%s' % colorByVal) region = regions(ax, cbar, self) ax.legend() return (ax, region) if refreshPlot is True: - ax, cbar = self.createPlot(dist, pick, color[colorByVal], label='%s' % colorByVal, ax=ax, cbar=cbar) + ax, cbar = self.createPlot(dist, pick, color[colorByVal], + label='%s' % colorByVal, ax=ax, + cbar=cbar) ax.legend() return ax @@ -574,27 +596,33 @@ class Survey(object): print('Generating new plot...') fig = plt.figure() ax = fig.add_subplot(111) - sc = ax.scatter(dist, pick, cmap=cm, c=inkByVal, s=5, edgecolors='none', label=label) + sc = ax.scatter(dist, pick, cmap=cm, c=inkByVal, s=5, + edgecolors='none', label=label) cbar = plt.colorbar(sc, fraction=0.05) cbar.set_label(label) ax.set_xlabel('Distance [m]') ax.set_ylabel('Time [s]') - ax.text(0.5, 0.95, 'Plot of all picks', transform=ax.transAxes, horizontalalignment='center') + ax.text(0.5, 0.95, 'Plot of all picks', transform=ax.transAxes, + horizontalalignment='center') else: - sc = ax.scatter(dist, pick, cmap=cm, c=inkByVal, s=5, edgecolors='none', label=label) + sc = ax.scatter(dist, pick, cmap=cm, c=inkByVal, s=5, + edgecolors='none', label=label) cbar = plt.colorbar(sc, cax=cbar.ax) cbar.set_label(label) ax.set_xlabel('Distance [m]') ax.set_ylabel('Time [s]') - ax.text(0.5, 0.95, 'Plot of all picks', transform=ax.transAxes, horizontalalignment='center') + ax.text(0.5, 0.95, 'Plot of all picks', transform=ax.transAxes, + horizontalalignment='center') return (ax, cbar) def _update_progress(self, shotname, tend, progress): - sys.stdout.write('Working on shot %s. ETC is %02d:%02d:%02d [%2.2f %%]\r' % (shotname, - tend.hour, - tend.minute, - tend.second, - progress)) + sys.stdout.write( + 'Working on shot %s. ETC is %02d:%02d:%02d [%2.2f %%]\r' % ( + shotname, + tend.hour, + tend.minute, + tend.second, + progress)) sys.stdout.flush() def saveSurvey(self, filename='survey.pickle'): diff --git a/pylot/core/active/seismicshot.py b/pylot/core/active/seismicshot.py index 827b12ac..3195af4a 100644 --- a/pylot/core/active/seismicshot.py +++ b/pylot/core/active/seismicshot.py @@ -11,6 +11,7 @@ from pylot.core.pick.charfuns import AICcf from pylot.core.pick.utils import getSNR from pylot.core.pick.utils import earllatepicker import matplotlib.pyplot as plt +import warnings plt.interactive('True') @@ -322,7 +323,7 @@ class SeismicShot(object): if len(traces) == 1: return Stream(traces) self.setPick(traceID, None) - print 'Warning: ambigious or empty traceID: %s' % traceID + warnings.warn('ambigious or empty traceID: %s' % traceID) def pickParallel(self, folm, method = 'hos', aicwindow = (10, 0)): import multiprocessing @@ -336,7 +337,11 @@ class SeismicShot(object): traceIDs = self.getTraceIDlist() - pool.map(self.pickTrace, traceIDs) + picks = pool.map(self.pickTrace, traceIDs) + + for traceID, pick in picks: + self.setPick(traceID, pick) + def pickTrace(self, traceID): ''' @@ -373,7 +378,7 @@ class SeismicShot(object): setHosAic = {'hos': hoscftime, 'aic': aiccftime} - self.setPick(traceID, setHosAic[self.getMethod()]) + return traceID, setHosAic[self.getMethod()] def setEarllatepick(self, traceID, nfac=1.5): tgap = self.getTgap() From 093f750aa1ea5234dbe4f85207740a841b6a8770 Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Mon, 23 May 2016 12:06:55 +0200 Subject: [PATCH 09/15] tried worker function --- pylot/core/active/seismicshot.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pylot/core/active/seismicshot.py b/pylot/core/active/seismicshot.py index 3195af4a..63f1911c 100644 --- a/pylot/core/active/seismicshot.py +++ b/pylot/core/active/seismicshot.py @@ -327,6 +327,7 @@ class SeismicShot(object): def pickParallel(self, folm, method = 'hos', aicwindow = (10, 0)): import multiprocessing + from pylot.core.util.utils import worker self.setFolm(folm) self.setMethod(method) @@ -337,11 +338,10 @@ class SeismicShot(object): traceIDs = self.getTraceIDlist() - picks = pool.map(self.pickTrace, traceIDs) + picks = worker(self.pickTrace, traceIDs, maxthreads) for traceID, pick in picks: self.setPick(traceID, pick) - def pickTrace(self, traceID): ''' From 25ca11f5724e8fdacc4bae98e52f17f3d36dc640 Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Mon, 23 May 2016 14:22:52 +0200 Subject: [PATCH 10/15] minor tweaks --- pylot/core/active/fmtomoUtils.py | 31 +++++++++++++++++++++++++++++-- 1 file changed, 29 insertions(+), 2 deletions(-) diff --git a/pylot/core/active/fmtomoUtils.py b/pylot/core/active/fmtomoUtils.py index 6db9705d..eac6e8c4 100644 --- a/pylot/core/active/fmtomoUtils.py +++ b/pylot/core/active/fmtomoUtils.py @@ -6,15 +6,19 @@ import datetime import numpy as np class Tomo3d(object): - def __init__(self, citer = 0, overwrite = False): + def __init__(self, fmtomodir, simuldir = 'fmtomo_simulation', citer = 0, overwrite = False): ''' Class build from FMTOMO script tomo3d. Can be used to run several instances of FMM code in parallel. :param: citer, current iteration (default = 0: start new model) :type: integer ''' + self.simuldir = simuldir self.setCWD() + self.buildFmtomodir(fmtomodir) + self.buildObsdata() self.defParas() + self.copyRef() self.citer = citer # current iteration self.sources = self.readSrcFile() self.traces = self.readTraces() @@ -25,6 +29,28 @@ class Tomo3d(object): self.defFMMParas() self.defInvParas() + def buildFmtomodir(self, directory): + tomo_files = ['fm3d', + 'frechgen', + 'frechgen.in', + 'invert3d', + 'invert3d.in', + 'mode_set.in', + 'obsdata', + 'obsdata.in', + 'residuals', + 'residuals.in', + 'tomo3d', + 'tomo3d.in'] + + for name in tomo_files: + filename = os.path.join(directory, name) + os.system('cp %s %s'%(filename, self.cwd)) + + def buildObsdata(self): + os.system('obsdata') + os.system('mv sources.in sourcesref.in') + def defFMMParas(self): ''' Initiates parameters for the forward calculation. @@ -89,8 +115,9 @@ class Tomo3d(object): Default: pwd ''' if directory == None: - directory = subprocess.check_output(['pwd'])[0:-1] + directory = os.path.join(os.getcwd(), self.simuldir) + os.chdir(directory) self.cwd = directory print('Working directory is: %s'%self.cwd) From 42cbfeb7872160b2467566da869579fabd632c75 Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Mon, 23 May 2016 14:24:02 +0200 Subject: [PATCH 11/15] temporary retaining of parallel tests --- pylot/core/active/seismicshot.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/pylot/core/active/seismicshot.py b/pylot/core/active/seismicshot.py index 63f1911c..3d2e7bc6 100644 --- a/pylot/core/active/seismicshot.py +++ b/pylot/core/active/seismicshot.py @@ -338,11 +338,15 @@ class SeismicShot(object): traceIDs = self.getTraceIDlist() - picks = worker(self.pickTrace, traceIDs, maxthreads) + # picks = worker(self.pickTrace, traceIDs, maxthreads) - for traceID, pick in picks: + # for traceID, pick in picks: + # self.setPick(traceID, pick) + + for traceID in traceIDs: + trID, pick = self.pickTrace(traceID) self.setPick(traceID, pick) - + def pickTrace(self, traceID): ''' Intitiate picking for a trace. From d9844fff1752986ec7806f4f688f84f3391f063d Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Mon, 23 May 2016 14:25:06 +0200 Subject: [PATCH 12/15] added worker --- pylot/core/util/utils.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/pylot/core/util/utils.py b/pylot/core/util/utils.py index c0eb0185..24282319 100644 --- a/pylot/core/util/utils.py +++ b/pylot/core/util/utils.py @@ -450,6 +450,12 @@ def runProgram(cmd, parameter=None): output = subprocess.check_output('{} | tee /dev/stderr'.format(cmd), shell=True) +def worker(func, input, cores): + from multiprocessing import Pool + pool = Pool(cores) + result = pool.map(func, input) + pool.close() + return result if __name__ == "__main__": import doctest From 46ccd44e16c61e82970729270c522c5c9972550d Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Mon, 23 May 2016 14:25:24 +0200 Subject: [PATCH 13/15] changed to generation of reference grids --- pylot/core/active/seismicArrayPreparation.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pylot/core/active/seismicArrayPreparation.py b/pylot/core/active/seismicArrayPreparation.py index 3fc10b9b..2e836f7d 100644 --- a/pylot/core/active/seismicArrayPreparation.py +++ b/pylot/core/active/seismicArrayPreparation.py @@ -401,7 +401,7 @@ class SeisArray(object): customgrid='mygrid.in', writeVTK=True): ''' Generate FMTOMO input files from the SeisArray dimensions. - Generates: vgrids.in, interfaces.in, propgrid.in + Generates: vgridsref.in, interfacesref.in, propgrid.in :param: nPointsPropgrid, number of points in each direction of the propagation grid (z, y, x) :type: tuple @@ -473,10 +473,10 @@ class SeisArray(object): outfile.close() def generateInterfaces(self, nTheta, nPhi, depthmax, cushionfactor=0.1, - outfilename='interfaces.in', method='linear', + outfilename='interfacesref.in', method='linear', returnInterfaces=False): ''' - Create an interfaces.in file for FMTOMO from the SeisArray boundaries. + Create an interfacesref.in file for FMTOMO from the SeisArray boundaries. :param: nTheta, number of points in Theta type: int @@ -610,7 +610,7 @@ class SeisArray(object): def generateVgrid(self, nTheta, nPhi, nR, Rbt, thetaSN=None, phiWE=None, cushionfactor=0.1, - outfilename='vgrids.in', method='linear', + outfilename='vgridsref.in', method='linear', infilename='mygrid.in', returnTopo=False): ''' Generate a velocity grid for fmtomo regarding topography with a linear gradient starting at the topography with 0.34 [km/s]. From a6eaac6c339b67ed3b2494438b8eb97743f11f94 Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Tue, 24 May 2016 14:19:37 +0200 Subject: [PATCH 14/15] Changes during parallelization tests of autopicker --- pylot/core/active/fmtomoUtils.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pylot/core/active/fmtomoUtils.py b/pylot/core/active/fmtomoUtils.py index eac6e8c4..d1c065bf 100644 --- a/pylot/core/active/fmtomoUtils.py +++ b/pylot/core/active/fmtomoUtils.py @@ -45,7 +45,8 @@ class Tomo3d(object): for name in tomo_files: filename = os.path.join(directory, name) - os.system('cp %s %s'%(filename, self.cwd)) + linkname = os.path.join(self.cwd, name) + os.system('ln -s %s %s'%(filename, linkname)) def buildObsdata(self): os.system('obsdata') From 8ca87bc7772625196c15db17c398f0ff926c41d3 Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Tue, 24 May 2016 14:20:59 +0200 Subject: [PATCH 15/15] changes while testing parallelization of autopicker --- pylot/core/active/activeSeismoPick.py | 24 +++++++++++++++++-- pylot/core/active/fmtomoUtils.py | 2 +- pylot/core/active/seismicshot.py | 33 ++++++++++++++------------- pylot/core/util/utils.py | 26 +++++++++++++++------ 4 files changed, 59 insertions(+), 26 deletions(-) diff --git a/pylot/core/active/activeSeismoPick.py b/pylot/core/active/activeSeismoPick.py index f2c18d83..3cd098b5 100644 --- a/pylot/core/active/activeSeismoPick.py +++ b/pylot/core/active/activeSeismoPick.py @@ -3,7 +3,11 @@ import sys import numpy as np from pylot.core.active import seismicshot from pylot.core.active.surveyUtils import cleanUp +import copy_reg +import types +from pylot.core.util.utils import worker, _pickle_method +copy_reg.pickle(types.MethodType, _pickle_method) class Survey(object): def __init__(self, path, sourcefile, receiverfile, useDefaultParas=False): @@ -196,6 +200,12 @@ class Survey(object): plt.xlabel('Difference in time (auto - manual) [s]') return diffs + def pickShot(self, shTr): + shotnumber, traceID = shTr + shot = self.getShotForShotnumber(shotnumber) + traceID, pick = shot.pickTrace(traceID) + return shotnumber, traceID, pick + def pickAllShots(self, vmin=333, vmax=5500, folm=0.6, HosAic='hos', aicwindow=(10, 0)): ''' @@ -218,20 +228,30 @@ class Survey(object): count = 0 tpicksum = starttime - starttime + shTr = [] + for shot in self.data.values(): tstartpick = datetime.now() shot.setVmin(vmin) shot.setVmax(vmax) count += 1 - shot.pickParallel(folm) + #shot.pickParallel(folm) + shot.setPickParameters(folm = folm, method = HosAic, aicwindow = aicwindow) + for traceID in shot.getTraceIDlist(): + shTr.append((shot.getShotnumber(), traceID)) + picks = worker(self.pickShot, shTr, async = True) + + for shotnumber, traceID, pick in picks.get(): + self.getShotForShotnumber(shotnumber).setPick(traceID, pick) + # tpicksum += (datetime.now() - tstartpick); # tpick = tpicksum / count # tremain = (tpick * (len(self.getShotDict()) - count)) # tend = datetime.now() + tremain # progress = float(count) / float(len(self.getShotDict())) * 100 # self._update_progress(shot.getShotname(), tend, progress) - + self.setSNR() self.setEarllate() diff --git a/pylot/core/active/fmtomoUtils.py b/pylot/core/active/fmtomoUtils.py index d1c065bf..94465a34 100644 --- a/pylot/core/active/fmtomoUtils.py +++ b/pylot/core/active/fmtomoUtils.py @@ -280,7 +280,7 @@ class Tomo3d(object): ''' Wipes a certain directory. ''' - print('Wiping directory %s...'%directory) + #print('Wiping directory %s...'%directory) for filename in os.listdir(directory): filenp = os.path.join(directory, filename) os.remove(filenp) diff --git a/pylot/core/active/seismicshot.py b/pylot/core/active/seismicshot.py index 3d2e7bc6..cbf1deb8 100644 --- a/pylot/core/active/seismicshot.py +++ b/pylot/core/active/seismicshot.py @@ -12,10 +12,14 @@ from pylot.core.pick.utils import getSNR from pylot.core.pick.utils import earllatepicker import matplotlib.pyplot as plt import warnings +import copy_reg +import types +from pylot.core.util.utils import worker, _pickle_method + +copy_reg.pickle(types.MethodType, _pickle_method) plt.interactive('True') - class SeismicShot(object): ''' SuperClass for a seismic shot object. @@ -325,27 +329,24 @@ class SeismicShot(object): self.setPick(traceID, None) warnings.warn('ambigious or empty traceID: %s' % traceID) - def pickParallel(self, folm, method = 'hos', aicwindow = (10, 0)): - import multiprocessing - from pylot.core.util.utils import worker - + def setPickParameters(self, folm, method = 'hos', aicwindow = (10, 0)): self.setFolm(folm) self.setMethod(method) self.setAicwindow(aicwindow) - maxthreads = multiprocessing.cpu_count() - pool = multiprocessing.Pool(maxthreads) - - traceIDs = self.getTraceIDlist() + # def pickParallel(self): + # traceIDs = self.getTraceIDlist() + # picks = [] + # #picks = worker(self.pickTrace, traceIDs) - # picks = worker(self.pickTrace, traceIDs, maxthreads) + # # for traceID, pick in picks: + # # self.setPick(traceID, pick) - # for traceID, pick in picks: - # self.setPick(traceID, pick) - - for traceID in traceIDs: - trID, pick = self.pickTrace(traceID) - self.setPick(traceID, pick) + # for traceID in traceIDs: + # trID, pick = self.pickTrace(traceID) + # picks.append([trID, pick]) + # #self.setPick(traceID, pick) + # return picks def pickTrace(self, traceID): ''' diff --git a/pylot/core/util/utils.py b/pylot/core/util/utils.py index 24282319..f2b5aae5 100644 --- a/pylot/core/util/utils.py +++ b/pylot/core/util/utils.py @@ -10,6 +10,25 @@ import numpy as np from obspy.core import UTCDateTime import obspy.core.event as ope +def _pickle_method(m): + if m.im_self is None: + return getattr, (m.im_class, m.im_func.func_name) + else: + return getattr, (m.im_self, m.im_func.func_name) + +def worker(func, input, cores = 'max', async = False): + import multiprocessing + + if cores == 'max': + cores = multiprocessing.cpu_count() + + pool = multiprocessing.Pool(cores) + if async == True: + result = pool.map_async(func, input) + else: + result = pool.map(func, input) + pool.close() + return result def createAmplitude(pickID, amp, unit, category, cinfo): ''' @@ -450,13 +469,6 @@ def runProgram(cmd, parameter=None): output = subprocess.check_output('{} | tee /dev/stderr'.format(cmd), shell=True) -def worker(func, input, cores): - from multiprocessing import Pool - pool = Pool(cores) - result = pool.map(func, input) - pool.close() - return result - if __name__ == "__main__": import doctest