Merge branch 'master' of ariadne.geophysik.ruhr-uni-bochum.de:/data/git/pylot

pylot/core/active/activeSeismoPick.py

@@ -4,11 +4,10 @@ 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):
         '''
-        The Survey Class contains all shots [type: seismicshot] of a survey
+        The Survey Class contains all shots [class: Seismicshot] of a survey
         as well as the aquisition geometry and the topography.
 
         It contains methods to pick all traces of all shots.
@@ -24,7 +23,7 @@ class Survey(object):
         self._generateSurvey()
         self._initiateFilenames()
         if useDefaultParas == True:
-            self.setParametersForShots()
+            self.setParametersForAllShots()
         self._removeAllEmptyTraces()
         self._updateShots()
 
@@ -51,32 +50,36 @@ class Survey(object):
         print ("Total number of traces: %d \n" % self.countAllTraces())
 
     def _removeAllEmptyTraces(self):
-        filename = 'removeEmptyTraces.out'
+        '''
+        Removes traces of the dataset that are not found in the input receiver files.
+        '''
+        logfile = 'removeEmptyTraces.out'
         count = 0
         for shot in self.data.values():
             removed = shot.removeEmptyTraces()
             if removed is not None:
-                if count == 0: outfile = open(filename, 'w')
+                if count == 0: outfile = open(logfile, 'w')
                 count += 1
                 outfile.writelines('shot: %s, removed empty traces: %s\n'
                                    % (shot.getShotnumber(), removed))
         print ("\nremoveEmptyTraces: Finished! Removed %d traces" % count)
         if count > 0:
             print ("See %s for more information "
-                   "on removed traces." % (filename))
+                   "on removed traces." % (logfile))
             outfile.close()
 
     def _updateShots(self):
         '''
-        Removes traces that do not exist in the dataset for any reason.
+        Removes traces that do not exist in the dataset for any reason,
+        but were set in the input files.
         '''
-        filename = 'updateShots.out'
+        logfile = 'updateShots.out'
         count = 0;
         countTraces = 0
         for shot in self.data.values():
             del_traceIDs = shot.updateTraceList()
             if len(del_traceIDs) > 0:
-                if count == 0: outfile = open(filename, 'w')
+                if count == 0: outfile = open(logfile, 'w')
                 count += 1
                 countTraces += len(del_traceIDs)
                 outfile.writelines("shot: %s, removed traceID(s) %s because "
@@ -87,36 +90,37 @@ class Survey(object):
                "%d traces" % (count, countTraces))
         if count > 0:
             print ("See %s for more information "
-                   "on removed traces." % (filename))
+                   "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.
-        (This can be used to create a pick with t = 0 at the source origin)
+        (Commonly used to generate a travel time t = 0 at the source origin)
         '''
         for shot in self.data.values():
             shot.setPick(traceID, pick)
 
-    def setParametersForShots(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 "
                    "setParamters. This might not be clever.")
-        # CHANGE this later. Parameters only needed for survey, not for each shot.
         for shot in self.data.values():
             shot.setCut(cutwindow)
             shot.setTmovwind(tmovwind)
             shot.setTsignal(tsignal)
             shot.setTgap(tgap)
             shot.setOrder(order=4)
-        print ("setParametersForShots: Parameters set to:\n"
+        print ("setParametersForAllShots: Parameters set to:\n"
                "cutwindow = %s, tMovingWindow = %f, tsignal = %f, tgap = %f"
                % (cutwindow, tmovwind, tsignal, tgap))
 
     def setManualPicksFromFiles(self, directory='picks'):
         '''
         Read manual picks from *.pck files in a directory.
+        Can be used for comparison of automatic and manual picks.
+
         The * must be identical with the shotnumber.
         '''
         for shot in self.data.values():
@@ -125,6 +129,7 @@ class Survey(object):
     def getDiffsFromManual(self):
         '''
         Returns a dictionary with the differences between manual and automatic pick for all shots.
+        Key: Seismicshot [object]
         '''
         diffs = {}
         for shot in self.data.values():
@@ -136,6 +141,10 @@ class Survey(object):
         return diffs
 
     def plotDiffs(self):
+        '''
+        Creates a plot of all Picks colored by the
+        difference between automatic and manual pick.
+        '''
         import matplotlib.pyplot as plt
         diffs = [];
         dists = [];
@@ -163,8 +172,12 @@ class Survey(object):
         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')
+        plt.legend()
 
     def plotHist(self, nbins=20, ax=None):
+        '''
+        Plot a histogram of the difference between automatic and manual picks.
+        '''
         import matplotlib.pyplot as plt
         plt.interactive(True)
         diffs = []
@@ -180,9 +193,21 @@ class Survey(object):
         plt.xlabel('Difference in time (auto - manual) [s]')
         return diffs
 
-    def pickAllShots(self, windowsize, HosAic='hos', vmin=333, vmax=5500, folm=0.6):
+    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.
+
+        :param: vmin, vmax, minimum (maximum) permitted apparent velocity on direct path between src and rec
+        :type: real
+
+        :param: folm, fraction of local maximum for HOS pick (0.6 = 60% of the highest maximum)
+        :type: real
+
+        :param: HosAic, pick with hos only ('hos') or use AIC ('aic')
+        :type: string
+
+        :param: aicwindow, window around the initial pick to search for local AIC min (samples)
+        :type: tuple
         '''
         from datetime import datetime
         starttime = datetime.now()
@@ -207,7 +232,7 @@ class Survey(object):
                     pickwin_used = (pwleft, pwright)
 
                 shot.setPickwindow(traceID, pickwin_used)
-                shot.pickTraces(traceID, windowsize, folm, HosAic)  # picker
+                shot.pickTraces(traceID, folm, HosAic, aicwindow)  # picker
 
                 shot.setSNR(traceID)
                 # if shot.getSNR(traceID)[0] < snrthreshold:
@@ -231,6 +256,9 @@ class Survey(object):
               % (pickedtraces, ntraces, float(pickedtraces) / float(ntraces) * 100.))
 
     def cleanBySPE(self, maxSPE):
+        '''
+        Sets all picks as invalid if they exceed a certain value of the symmetric pick error.
+        '''
         for shot in self.data.values():
             for traceID in shot.getTraceIDlist():
                 if shot.getPickFlag(traceID) == 1:
@@ -238,6 +266,9 @@ class Survey(object):
                         shot.setPickFlag(traceID, 0)
 
     def plotSPE(self):
+        '''
+        Plots the symmetric pick error sorted by magnitude.
+        '''
         import matplotlib.pyplot as plt
         spe = []
         for shot in self.data.values():
@@ -251,7 +282,7 @@ class Survey(object):
 
     def recover(self):
         '''
-        Recovers all (accidently) removed picks. Still regards SNR threshold.
+        Recovers all manually removed picks. Still regards SNR threshold.
         '''
         print('Recovering survey...')
         numpicks = 0
@@ -266,18 +297,28 @@ class Survey(object):
         print('Recovered %d picks' % numpicks)
 
     def setArtificialPick(self, traceID, pick):
+        '''
+        Sets an artificial pick for a certain receiver (traceID) for all shots.
+        '''
         for shot in self.data.values():
             shot.setPick(traceID, pick)
             shot.setPickwindow(traceID, shot.getCut())
 
     def countAllTraces(self):
+        '''
+        Returns the number of traces in total.
+        '''
         numtraces = 0
         for shot in self.getShotlist():
-            for rec in self.getReceiverlist():  ### shot.getReceiverlist etc.
+            for rec in self.getReceiverlist():
                 numtraces += 1
+
         return numtraces
 
     def getShotlist(self):
+        '''
+        Returns a list of all shotnumbers contained in the set Sourcefile.
+        '''
         filename = self.getPath() + self.getSourcefile()
         srcfile = open(filename, 'r')
         shotlist = []
@@ -288,6 +329,9 @@ class Survey(object):
         return shotlist
 
     def getReceiverlist(self):
+        '''
+        Returns a list of all trace IDs contained in the set Receiverfile.
+        '''
         filename = self.getPath() + self.getReceiverfile()
         recfile = open(filename, 'r')
         reclist = []
@@ -313,6 +357,12 @@ class Survey(object):
         return self._obsdir
 
     def getStats(self):
+        '''
+        Generates and returns a dictionary containing statistical information
+        of the survey.
+        
+        Key: shotnumber
+        '''
         info_dict = {}
         for shot in self.data.values():
             pickedTraces = 0
@@ -334,11 +384,17 @@ class Survey(object):
         return info_dict
 
     def getShotForShotnumber(self, shotnumber):
+        '''
+        Returns Seismicshot [object] of a certain shotnumber if possible.
+        '''
         for shot in self.data.values():
             if shot.getShotnumber() == shotnumber:
                 return shot
 
     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.
@@ -354,13 +410,13 @@ class Survey(object):
         srcfile.writelines('%10s\n' % len(self.data))  # number of sources
         for shotnumber in self.getShotlist():
             shot = self.getShotForShotnumber(shotnumber)
-            ttfilename = str(shotnumber) + ttFileExtension
+            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))  # lat, lon, depth
+            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\n' % 1)
             srcfile.writelines('%10s %10s %10s\n' % (1, 1, ttfilename))
             ttfile = open(directory + '/' + ttfilename, 'w')
             traceIDlist = shot.getTraceIDlist()
@@ -393,6 +449,9 @@ class Survey(object):
         print(msg)
 
     def countPickedTraces(self, shot):
+        '''
+        Counts all picked traces of a shot (type Seismicshot).
+        '''
         count = 0
         for traceID in shot.getTraceIDlist():
             if shot.getPickFlag(traceID) is not 0:
@@ -400,6 +459,9 @@ class Survey(object):
         return count
 
     def countAllPickedTraces(self):
+        '''
+        Counts all picked traces of the survey.
+        '''
         count = 0
         for shot in self.data.values():
             for traceID in shot.getTraceIDlist():
@@ -422,20 +484,9 @@ class Survey(object):
         figPerSubplot = columns * rows
 
         index = 1
-        # shotnames = []
-        # shotnumbers = []
-
-        # for shot in self.data.values():
-        #     shotnames.append(shot.getShotname())
-        #     shotnumbers.append(shot.getShotnumber())
-
-        # shotnumbers = [shotnumbers for (shotnumbers, shotnames) in sorted(zip(shotnumbers, shotnames))]
 
         for shotnumber in self.getShotlist():
             if index <= figPerSubplot:
-                # ax = fig.add_subplot(3,3,i, projection = '3d', title = 'shot:'
-                # +str(shot_dict[shotnumber].getShotnumber()), xlabel = 'X', ylabel = 'Y', zlabel = 'traveltime')
-                # shot_dict[shotnumber].plot3dttc(ax = ax, plotpicks = True)
                 ax = fig.add_subplot(rows, columns, index)
                 if mode == '3d':
                     self.getShot(shotnumber).matshow(ax=ax, colorbar=False, annotations=True, legend=False)
@@ -516,6 +567,9 @@ class Survey(object):
             return ax
 
     def createPlot(self, dist, pick, inkByVal, label, ax=None, cbar=None):
+        '''
+        Used by plotAllPicks.
+        '''
         import matplotlib.pyplot as plt
         plt.interactive(True)
         cm = plt.cm.jet
@@ -547,6 +601,10 @@ class Survey(object):
         sys.stdout.flush()
 
     def saveSurvey(self, filename='survey.pickle'):
+        '''
+        Save Survey object to a file. 
+        Can be loaded by using Survey.from_pickle(filename).
+        '''
         import cPickle
         cleanUp(self)
         outfile = open(filename, 'wb')

pylot/core/active/fmtomoUtils.py

@@ -6,15 +6,28 @@ import datetime
 import numpy as np
 
 class Tomo3d(object):
-    def __init__(self, nproc, iterations):
+    def __init__(self, nproc, iterations, 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 = 0           # current iteration
+        self.citer = citer       # current iteration
         self.sources = self.readSrcFile()
         self.traces = self.readTraces()
         self.directories = []
+        self.overwrite = overwrite
         
     def defParas(self):
         self.defFMMParas()
@@ -29,6 +42,7 @@ class Tomo3d(object):
         self.pg = 'propgrid.in'
         self.rec = 'receivers.in'
         self.frech = 'frechet.in'
+        self.frechout = 'frechet.dat'
         self.ot = 'otimes.dat'
         self.ttim = 'arrivals.dat'
         self.mode = 'mode_set.in'
@@ -53,8 +67,8 @@ class Tomo3d(object):
         self.resout = 'residuals.dat'
 
     def copyRef(self):
-        # Attention: 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.ivg, self.cvg))
         os.system('cp %s %s'%(self.iig, self.cig))
         os.system('cp %s %s'%(self.isl, self.csl))
 
@@ -65,6 +79,30 @@ class Tomo3d(object):
     def runFrech(self):
         os.system(self.frechgen)
 
+    def runTOMO3D(self):
+        starttime = datetime.datetime.now()
+        print('Starting TOMO3D on %s parallel processes for %s iteration(s).'
+              %(self.nproc, self.iter))
+        if self.citer == 0:
+            self.makeInvIterDir()
+            self.startForward(self.cInvIterDir)
+            self.raiseIter()
+        
+        while self.citer <= self.iter:
+            self.makeInvIterDir()
+            self.startInversion()
+            self.saveVgrid()
+            self.startForward(self.cInvIterDir)
+            self.raiseIter()
+
+        if self.citer > self.iter:
+            self.removeDirectories()
+            self.unlink(os.path.join(self.cwd, self.frechout))
+            self.unlink(os.path.join(self.cwd, self.ttim))
+
+        tdelta = datetime.datetime.now() - starttime
+        print('runTOMO3D: Finished %s iterations after %s.'%(self.iter, tdelta))
+
     def runFmm(self, directory, logfile, processes):
         os.chdir(directory)
         fout = open(logfile, 'w')
@@ -73,6 +111,47 @@ class Tomo3d(object):
         os.chdir(self.cwd)
         return processes
 
+    def startForward(self, logdir):
+        self._printLine()
+        print('Starting forward simulation for iteration %s.'%(self.citer))
+
+        if self.citer == 0:
+            self.copyRef()
+            self.runFrech()
+            self.makeDirectories()
+
+        starttime = datetime.datetime.now()
+        processes = []
+
+        for procID in range(1, self.nproc + 1):
+            directory = self.getProcDir(procID)
+            logfn = 'fm3dlog_' + str(procID) + '.out'
+            log_out = os.path.join(logdir, logfn)
+            
+            self.writeSrcFile(procID, directory)
+            self.writeTracesFile(procID, directory)
+            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))
+            processes = self.runFmm(directory, log_out, processes)
+            
+        for p in processes:
+            p.wait()
+
+        self.mergeOutput(self.cInvIterDir)
+        self.clearDirectories()
+        self.copyArrivals()
+        if self.citer == 0:
+            self.copyArrivals(self.rtrav)
+
+        self.calcRes()
+        tdelta = datetime.datetime.now() - starttime
+        print('Finished Forward calculation after %s'%tdelta)
+
+    def startInversion(self):
+        print('Calling %s...'%self.inv)
+        os.system(self.inv)
+
     def calcRes(self):
         resout = os.path.join(self.cwd, self.resout)
         if self.citer == 0:
@@ -85,111 +164,91 @@ class Tomo3d(object):
         RMS, var, chi2 = residuals[-1].split()
         print('Residuals: RMS = %s, var = %s, Chi^2 = %s.'%(RMS, var, chi2))
 
-    def runTOMO3D(self):
-        starttime = datetime.datetime.now()
-        print('Starting TOMO3D on %s parallel processes for %s iteration(s).'
-              %(self.nproc, self.iter))
-        if self.citer == 0:
-            self.startForward()
-            self.raiseIter()
-        
-        while self.citer <= self.iter:
-            self.startInversion()
-            self.startForward()
-            self.raiseIter()
-
-        tdelta = datetime.datetime.now() - starttime
-        print('runTOMO3D: Finished %s iterations after %s.'%(self.iter, tdelta))
-
-    def _printLine(self):
-        print('----------------------------------------')
-
     def raiseIter(self):
         self.citer +=1
         self._printLine()
+        invfile = open(self.cwd + '/inviter.in', 'w')
+        invfile.write('%s'%self.citer)
+        invfile.close()
 
-    def startInversion(self):
+    def makeDir(self, directory):
-        print('Calling %s...'%self.inv)
-        os.system(self.inv)
-            
-    def startForward(self):
-        self._printLine()
-        print('Starting forward simulation for iteration %s.'%(self.citer))
-        self.makeDirectories()
-        starttime = datetime.datetime.now()
-        processes = []
-
-        if self.citer == 0:
-            self.copyRef()
-            self.runFrech()
-
-        for procID in range(1, self.nproc + 1):
-            directory = self.getProcDir(procID)
-            log_out = self.cwd + '/fm3dlog_' + str(procID) + '.out'
-            
-            self.writeSrcFile(procID, directory)
-            self.writeTracesFile(procID, directory)
-            os.system('cp {cvg} {cig} {mode} {pg} {frechout} {dest}'
-                      .format(cvg=self.cvg, cig=self.cig, frechout=self.frech,
-                              mode=self.mode, pg=self.pg, dest=directory))
-            processes = self.runFmm(directory, log_out, processes)
-            
-        for p in processes:
-            p.wait()
-
-        self.mergeOutput()
-        self.clearDirectories()
-        self.copyArrivals()
-        if self.citer == 0:
-            self.copyArrivals(self.rtrav)
-
-        self.calcRes()
-        tdelta = datetime.datetime.now() - starttime
-        print('Finished Forward calculation after %s'%tdelta)
-
-    def copyArrivals(self, target = None):
-        if target == None:
-            target = self.mtrav
-        os.system('cp %s %s'%(self.ttim, target))
-
-    def makeDIR(self, directory):
         err = os.system('mkdir %s'%directory)
+        if err is 0:
+            self.directories.append(directory)
+            return
         if err is 256:
-            response = raw_input('Warning: Directory already existing. Continue (y/n)?\n')
+            if self.overwrite == True:
-            if response == 'y':
                 print('Overwriting existing files.')
-            else:
+                self.clearDir(directory)
-                sys.exit('Aborted')
+                self.directories.append(directory)
+                return
+        raise RuntimeError('Could not create directory: %s'%directory)
 
-        self.directories.append(directory)
-
-    def clearDIR(self, directory):
-        err = os.system('rm -r %s'%directory)
-        
     def makeDirectories(self):
         for procID in range(1, self.nproc + 1):
             directory = self.getProcDir(procID)
-            self.makeDIR(directory) 
+            self.makeDir(directory) 
 
+    def makeInvIterDir(self):
+        invIterDir = self.cwd + '/it_%s'%(self.citer)
+        err = os.system('mkdir %s'%invIterDir)
+        if err is 256:
+            if self.overwrite:
+                self.clearDir(invIterDir)
+        elif err is not 0:
+            raise RuntimeError('Could not create directory: %s'%invIterDir)
+        self.cInvIterDir = invIterDir        
+
+    def clearDir(self, directory):
+        print('Wiping directory %s...'%directory)
+        for filename in os.listdir(directory):
+            filenp = os.path.join(directory, filename)
+            os.remove(filenp)
+        
     def clearDirectories(self):
         for directory in self.directories:
-            self.clearDIR(directory) 
+            self.clearDir(directory) 
+
+    def rmDir(self, directory):
+        print('Removing directory %s...'%directory)
+        return os.rmdir(directory)
+
+    def removeDirectories(self):
+        for directory in self.directories:
+            self.rmDir(directory) 
         self.directories = []
         
-    def readNsrc(self):
+    def getProcDir(self, procID):
-        srcfile = open(self.csl, 'r')
+        return os.path.join(self.cwd, self.folder) + str(procID)
-        nsrc = int(srcfile.readline())
-        srcfile.close()
-        return nsrc
 
-    def readNtraces(self):
+    def getTraceIDs4Sources(self, sourceIDs):
-        recfile = open(self.rec, 'r')
+        traceIDs = []
-        nrec = int(recfile.readline())
+        for traceID in self.traces.keys():
-        recfile.close()
+            if self.traces[traceID]['source'] in sourceIDs:
-        return nrec
+                traceIDs.append(traceID)
+        return traceIDs
+
+    def getTraceIDs4Source(self, sourceID):
+        traceIDs = []
+        for traceID in self.traces.keys():
+            if self.traces[traceID]['source'] == sourceID:
+                traceIDs.append(traceID)
+        return traceIDs
+
+    def copyArrivals(self, target = None):
+        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')
+        os.system('cp %s %s'%(vgpath, self.cInvIterDir))
 
     def calcSrcPerKernel(self):
         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):
@@ -207,6 +266,18 @@ class Tomo3d(object):
         elif proc > remain:
             start = (srcPK + 1) * remain + srcPK * (proc - remain) + 1
             return range(start, start + srcPK)
+    
+    def readNsrc(self):
+        srcfile = open(self.csl, 'r')
+        nsrc = int(srcfile.readline())
+        srcfile.close()
+        return nsrc
+
+    def readNtraces(self):
+        recfile = open(self.rec, 'r')
+        nrec = int(recfile.readline())
+        recfile.close()
+        return nrec
 
     def readSrcFile(self):
         nsrc = self.readNsrc()
@@ -258,6 +329,62 @@ class Tomo3d(object):
             
         return traces
 
+    def readArrivals(self, procID):
+        directory = self.getProcDir(procID)
+        arrfile = open(directory + '/arrivals.dat', 'r')
+        sourceIDs = self.srcIDs4Kernel(procID)
+        
+        arrivals = []
+        for sourceID in sourceIDs:
+            traceIDs = self.getTraceIDs4Source(sourceID)
+            for traceID in traceIDs:
+                line = arrfile.readline().split()
+                if line != []:
+                    # recID and srcID for the individual processor will not be needed
+                    recID_proc, srcID_proc, ray, normal, arrtime, diff, head = line
+                    arrivals.append([traceID, sourceID, ray, normal, arrtime, diff, head])
+                
+        return arrivals
+
+    def readRays(self, procID):
+        directory = self.getProcDir(procID)
+        raysfile = open(directory + '/rays.dat', 'r')
+        sourceIDs = self.srcIDs4Kernel(procID)
+        
+        rays = {}
+        for sourceID in sourceIDs:
+            traceIDs = self.getTraceIDs4Source(sourceID)
+            for traceID in traceIDs:
+                line1 = raysfile.readline().split()
+                if line1 != []:
+                    # recID and srcID for the individual processor will not be needed
+                    recID_proc, srcID_proc, ray, normal, nsec = line1
+                    raysecs = {}
+
+                    for sec in range(int(nsec)):
+                        line2 = raysfile.readline().split()
+                        npoints, region, diff, head = line2
+                        raypoints = []
+                    
+                        for j in range(int(npoints)):
+                            raypoints.append(raysfile.readline() + '\n')
+
+                        raysecs[sec] = {'npoints': npoints,
+                                        'region': region,
+                                        'diff': diff,
+                                        'head': head,
+                                        'raypoints': raypoints
+                                        }
+
+
+                    rays[traceID] = {'sourceID': sourceID,
+                                     'raypath': ray,
+                                     'normal': normal,
+                                     'nsec': nsec,
+                                     'raysections': raysecs
+                                     }
+        return rays
+
     def writeSrcFile(self, procID, directory):
         srcfile = open('%s/sources.in'%directory, 'w')
         sourceIDs = self.srcIDs4Kernel(procID)
@@ -289,154 +416,44 @@ class Tomo3d(object):
             recfile.write('%s\n'%source)
             recfile.write('%s\n'%trace['path'])
 
-    def getTraceIDs4Sources(self, sourceIDs):
+    def mergeArrivals(self, directory):
-        traceIDs = []
+        arrfn = os.path.join(directory, self.ttim)
-        for traceID in self.traces.keys():
+        arrivalsOut = open(arrfn, 'w')
-            if self.traces[traceID]['source'] in sourceIDs:
-                traceIDs.append(traceID)
-        return traceIDs
-
-    def getTraceIDs4Source(self, sourceID):
-        traceIDs = []
-        for traceID in self.traces.keys():
-            if self.traces[traceID]['source'] == sourceID:
-                traceIDs.append(traceID)
-        return traceIDs
-
-    def readArrivals(self, procID):
-        directory = self.getProcDir(procID)
-        arrfile = open(directory + '/arrivals.dat', 'r')
-        sourceIDs = self.srcIDs4Kernel(procID)
-        
-        arrivals = []
-        for sourceID in sourceIDs:
-            traceIDs = self.getTraceIDs4Source(sourceID)
-            for traceID in traceIDs:
-                line = arrfile.readline().split()
-                if line != []:
-                    # recID and srcID for the individual processor will not be needed
-                    recID_proc, srcID_proc, ray, normal, arrtime, diff, head = line
-                    arrivals.append([traceID, sourceID, ray, normal, arrtime, diff, head])
-                
-        return arrivals
-
-    def mergeArrivals(self):
-        arrivalsOut = open(os.path.join(self.cwd, self.ttim), 'w')
         print('Merging arrivals.dat...')
         for procID in range(1, self.nproc + 1):
             arrivals = self.readArrivals(procID)
             for line in arrivals:
                 arrivalsOut.write('%6s %6s %6s %6s %15s %5s %5s\n'%tuple(line))
 
-    # def mergeFrechet(self):
+        os.system('ln -fs %s %s'%(arrfn, os.path.join(self.cwd, self.ttim)))
-    #     print('Merging frechet.dat...')
-    #     frechetOut = open(self.cwd + '/frechet.dat', 'w')
-    #     for procID in range(1, self.nproc + 1):
-    #         frechet = self.readFrechet(procID)
-    #         traceIDs = frechet.keys()
-    #         traceIDs.sort()
-    #         for traceID in traceIDs:
-    #             frech = frechet[traceID]
-    #             frechetOut.write('%6s %6s %6s %6s %6s\n'%
-    #                              (traceID,
-    #                               frech['sourceID'],
-    #                               frech['raypath'],
-    #                               frech['normal'],
-    #                               frech['NPDEV']))
-    #             for pdev in frech['PDEV']:
-    #                 frechetOut.writelines(pdev)
 
-    def readRays(self, procID):
+    def mergeRays(self, directory):
-        directory = self.getProcDir(procID)
-        raysfile = open(directory + '/rays.dat', 'r')
-        sourceIDs = self.srcIDs4Kernel(procID)
-        
-        rays = {}
-        for sourceID in sourceIDs:
-            traceIDs = self.getTraceIDs4Source(sourceID)
-            for traceID in traceIDs:
-                line1 = raysfile.readline().split()
-                if line1 != []:
-                    # recID and srcID for the individual processor will not be needed
-                    recID_proc, srcID_proc, ray, normal, nsec = line1
-                    raysecs = {}
-
-                    for sec in range(int(nsec)):
-                        line2 = raysfile.readline.split()
-                        npoints, region, diff, head = line2
-                        raypoints = []
-                    
-                        for j in range(int(npoints)):
-                            raypoints.append(raysfile.readline() + '\n')
-
-                        raysecs[sec] = {'npoints': npoints,
-                                        'region': region,
-                                        'diff': diff,
-                                        'head': head,
-                                        'raypoints': raypoints
-                                        }
-
-
-                    rays[traceID] = {'sourceID': sourceID,
-                                     'raypath': ray,
-                                     'normal': normal,
-                                     'nsec': nsec,
-                                     'raysections': raysecs
-                                     }
-        return rays
-
-    def mergeRays(self):
         print('Merging rays.dat...')
-        with open(self.cwd + 'rays.dat', 'w') as outfile:
+        with open(directory + '/rays.dat', 'w') as outfile:
-        for procID in range(1, self.nproc + 1):
-            rays = self.readRays(procID)
-            for traceID in rays:
-                ray = rays[traceID]
-                outfile.write('%6s %6s %6s %6s %6s'%(traceID,
-                                                     ray['sourceID'],
-                                                     ray['raypath'],
-                                                     ray['normal'],
-                                                     ray['nsec']))
-                for sec in range(int(ray['nsec'])):
-                    raysec = ray['raysections'][sec]
-                    outfile.write('%6s %6s %6s %6s'%(raysec['npoints'],
-                                                     raysec['region'],
-                                                     raysec['diff'],
-                                                     raysec['head']))
-                    outfile.writelines(raysec['raypoints'])
-                                                  
-
-    # def readFrechet(self, procID):
-    #     directory = self.getProcDir(procID)
-    #     frechfile = open(directory + '/frechet.dat', 'r')
-    #     sourceIDs = self.srcIDs4Kernel(procID)
-        
-    #     frechet = {}
-    #     for sourceID in sourceIDs:
-    #         traceIDs = self.getTraceIDs4Source(sourceID)
-    #         for traceID in traceIDs:
-    #             line = frechfile.readline().split()
-    #             if line != []:
-    #                 # recID and srcID for the individual processor will not be needed
-    #                 PDEV = []
-    #                 recID_proc, srcID_proc, ray, normal, NPDEV = line
-    #                 for i in range(int(NPDEV)):
-    #                     PDEV.append(frechfile.readline())
-
-    #                 frechet[traceID] = {'sourceID': sourceID,
-    #                                     'raypath': ray,
-    #                                     'normal': normal,
-    #                                     'NPDEV': NPDEV,
-    #                                     'PDEV': PDEV
-    #                                     }
-    #     return frechet
-
-    def mergeFrechet(self):
-        print('Merging frechet.dat...')
-
-        with open(self.cwd + '/frechet.dat', 'w') as outfile:
             for procID in range(1, self.nproc + 1):
-                filename = self.getProcDir(procID) + '/frechet.dat'
+                rays = self.readRays(procID)
+                for traceID in rays:
+                    ray = rays[traceID]
+                    outfile.write('%6s %6s %6s %6s %6s\n'%(traceID,
+                                                           ray['sourceID'],
+                                                           ray['raypath'],
+                                                           ray['normal'],
+                                                           ray['nsec']))
+                    for sec in range(int(ray['nsec'])):
+                        raysec = ray['raysections'][sec]
+                        outfile.write('%6s %6s %6s %6s\n'%(raysec['npoints'],
+                                                           raysec['region'],
+                                                           raysec['diff'],
+                                                           raysec['head']))
+                        outfile.writelines(raysec['raypoints'])
+                                                  
+    def mergeFrechet(self, directory):
+        print('Merging frechet.dat...')
+        
+        frechfnout = os.path.join(directory, 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)
                 with open(filename) as infile:
                     for sourceID in self.srcIDs4Kernel(procID):
                         for traceID in self.getTraceIDs4Source(sourceID):
@@ -445,15 +462,18 @@ class Tomo3d(object):
                             for index in range(int(NPDEV)):
                                 outfile.write(infile.readline())
 
+        os.system('ln -fs %s %s'%(frechfnout, os.path.join(self.cwd, self.frechout)))
 
-    def getProcDir(self, procID):
+    def mergeOutput(self, directory):
-        return os.path.join(self.cwd, self.folder) + str(procID)
+        self.mergeArrivals(directory)
-
+        self.mergeFrechet(directory)
-    def mergeOutput(self):
+        self.mergeRays(directory)
-        self.mergeArrivals()
-        self.mergeFrechet()
-        self.mergeRays()
 
+    def unlink(self, filepath):
+        return os.system('unlink %s'%filepath)
+        
+    def _printLine(self):
+        print('----------------------------------------')
 
 def vgrids2VTK(inputfile='vgrids.in', outputfile='vgrids.vtk', absOrRel='abs', inputfileref='vgridsref.in'):
     '''

pylot/core/io/phases.py

@@ -193,13 +193,13 @@ def picksdict_from_obs(fn):
 
 
 def picksdict_from_picks(evt):
-    '''
+    """
     Takes an Event object and return the pick dictionary commonly used within
     PyLoT
     :param evt: Event object contain all available information
     :type evt: `~obspy.core.event.Event`
     :return: pick dictionary
-    '''
+    """
     picks = {}
     for pick in evt.picks:
         phase = {}
@@ -345,7 +345,7 @@ def reassess_pilot_event(root_dir, event_id, out_dir=None, fn_param=None, verbos
 
 
 def writephases(arrivals, fformat, filename):
-    '''
+    """
     Function of methods to write phases to the following standard file
     formats used for locating earthquakes:
 
@@ -363,7 +363,7 @@ def writephases(arrivals, fformat, filename):
 
     :param: filename, full path and name of phase file
     :type: string
-    '''
+    """
 
     if fformat == 'NLLoc':
         print ("Writing phases to %s for NLLoc" % filename)
@@ -425,7 +425,6 @@ def writephases(arrivals, fformat, filename):
                                                                                            sweight))
 
         fid.close()
-
     elif fformat == 'HYPO71':
         print ("Writing phases to %s for HYPO71" % filename)
         fid = open("%s" % filename, 'w')