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

pylot/core/pick/compare.py

@@ -307,7 +307,7 @@ class PDFDictionary(object):
             pdfs = self.pdf_data[station]
             for l, phase in enumerate(pdfs.keys()):
                 try:
-                    axarr[n, l].plot(pdfs[phase].axis, pdfs[phase].data)
+                    axarr[n, l].plot(pdfs[phase].axis, pdfs[phase].data())
                     if n is 0:
                         axarr[n, l].set_title(phase)
                     if l is 0:
@@ -322,7 +322,7 @@ class PDFDictionary(object):
                 except IndexError as e:
                     print('trying aligned plotting\n{0}'.format(e))
                     hide_labels = False
-                    axarr[l].plot(pdfs[phase].axis, pdfs[phase].data)
+                    axarr[l].plot(pdfs[phase].axis, pdfs[phase].data())
                     axarr[l].set_title(phase)
                     if l is 0:
                         axann = axarr[l].annotate(station, xy=(.05, .5),
@@ -350,13 +350,15 @@ class PDFstatistics(object):
         self.stations = {}
         self.p_std = {}
         self.s_std = {}
-        self.makeDirlist()
+        self.theta015 = []
-        self.getData()
+        self.theta1 = []
-        self.getStatistics()
+        self.theta2 = []
-        self.arraylen = 0
+        self.dirlist = list()
-        self.Theta015 = []
+        self.evtdict = {}
-        self.Theta25 = []
+        #self.makeDirlist()
-        self.Theta485 = []
+        #self.getData()
+        #self.getStatistics()
+
         #self.showData()
 
 
@@ -367,8 +369,10 @@ class PDFstatistics(object):
             self.evtdict[rd] = glob.glob1((os.path.join(self.directory, rd)), '*.xml')
 
     def getData(self):
+        arraylen = 0
         for dir in self.dirlist:
             for evt in self.evtdict[dir]:
+                print evt
                 self.stations[evt] = []
                 self.p_std[evt] = []
                 self.s_std[evt] = []
@@ -377,16 +381,24 @@ class PDFstatistics(object):
                     # print station, pdfs
                     try:
                         p_std = pdfs['P'].standard_deviation()
+                        self.theta015.append(pdfs['P'].qtile_dist_quot(0.015))
+                        self.theta1.append(pdfs['P'].qtile_dist_quot(0.1))
+                        self.theta2.append(pdfs['P'].qtile_dist_quot(0.2))
                     except KeyError:
                         p_std = np.nan
                     try:
                         s_std = pdfs['S'].standard_deviation()
+                        self.theta015.append(pdfs['S'].qtile_dist_quot(0.015))
+                        self.theta1.append(pdfs['S'].qtile_dist_quot(0.1))
+                        self.theta2.append(pdfs['S'].qtile_dist_quot(0.2))
                     except KeyError:
                         s_std = np.nan
                     self.stations[evt].append(station)
                     self.p_std[evt].append(p_std)
                     self.s_std[evt].append(s_std)
-                    self.arraylen += 1
+                    arraylen += 1
+        self.arraylen = arraylen
+        self.makeArray()
 
 
     def makeArray(self):
@@ -426,5 +438,11 @@ class PDFstatistics(object):
 
 
 if __name__ == "__main__":
-    rootdir = '/home/sebastianp/Data/Reassessment/Insheim/'
+    rootdir = '/home/sebastianp/Data/Reassessment/Insheim/2012.10/'
-    Insheim = PDFstatistics(rootdir)
+    Insheim = PDFstatistics(rootdir)
+    Insheim.dirlist = [rootdir]
+    Insheim.evtdict[rootdir] = ['e0002.294.12.xml']
+    Insheim.getData()
+    print Insheim.theta015
+    print Insheim.theta1
+    print Insheim.theta2

pylot/core/util/pdf.py

@@ -149,8 +149,8 @@ class ProbabilityDensityFunction(object):
         x0, incr, npts = self.commonparameter(other)
 
         axis = create_axis(x0, incr, npts)
-        pdf_self = np.array([self.data(x) for x in axis])
+        pdf_self = np.array(self.data(axis))
-        pdf_other = np.array([other.data(x) for x in axis])
+        pdf_other = np.array(other.data(axis))
 
         pdf = np.convolve(pdf_self, pdf_other, 'full') * incr
 
@@ -174,8 +174,8 @@ class ProbabilityDensityFunction(object):
         x0, incr, npts = self.commonparameter(other)
 
         axis = create_axis(x0, incr, npts)
-        pdf_self = np.array([self.data(x) for x in axis])
+        pdf_self = np.array(self.data(axis))
-        pdf_other = np.array([other.data(x) for x in axis])
+        pdf_other = np.array(other.data(axis))
 
         pdf = np.correlate(pdf_self, pdf_other, 'full') * incr
 
@@ -193,21 +193,22 @@ class ProbabilityDensityFunction(object):
 
     def __nonzero__(self):
         prec = self.precision(self.incr)
-        data = np.array([self.data(t) for t in self.axis])
+        data = np.array(self.data())
         gtzero = np.all(data >= 0)
         probone = bool(np.round(self.prob_gt_val(self.axis[0]), prec) == 1.)
         return bool(gtzero and probone)
 
     def __str__(self):
-        return str([self.data(val) for val in create_axis(self.x0, self.incr,
+        return str([self.data()])
-                                                          self.npts)])
 
     @staticmethod
     def precision(incr):
         prec = int(np.ceil(np.abs(np.log10(incr)))) - 2
         return prec if prec >= 0 else 0
 
-    def data(self, value):
+    def data(self, value=None):
+        if value is None:
+            return self._pdf(self.axis, self.params)
         return self._pdf(value, self.params)
 
     @property
@@ -325,14 +326,15 @@ class ProbabilityDensityFunction(object):
             raise ValueError('value out of bounds: {0}'.format(value))
         return self.prob_limits((value, self.axis[-1]))
 
-    def prob_limits(self, limits):
+    def prob_limits(self, limits, oversampling=1.):
-        lim = np.arange(limits[0], limits[1], self.incr)
+        sampling = self.incr / oversampling
-        data = [self.data(t) for t in lim]
+        lim = np.arange(limits[0], limits[1], sampling)
+        data = self.data(lim)
         min_est, max_est = 0., 0.
         for n in range(len(data) - 1):
             min_est += min(data[n], data[n + 1])
             max_est += max(data[n], data[n + 1])
-        return (min_est + max_est) / 2. * self.incr
+        return (min_est + max_est) / 2. * sampling
 
     def prob_val(self, value):
         if not (self.axis[0] <= value <= self.axis[-1]):
@@ -352,7 +354,7 @@ class ProbabilityDensityFunction(object):
         m = (r + l) / 2
         diff = prob_value - self.prob_lt_val(m)
 
-        while abs(diff) > eps:
+        while abs(diff) > eps and ((r - l) > self.incr):
             if diff > 0:
                 l = m
             else:
@@ -367,12 +369,17 @@ class ProbabilityDensityFunction(object):
 
         return qu - ql
 
+
+    def qtile_dist_quot(self,x):
+        if x < 0 or x > 0.5:
+            raise ValueError('Value out of range.')
+        return self.quantile_distance(0.5-x)/self.quantile_distance(x)
+
+
     def plot(self, label=None):
         import matplotlib.pyplot as plt
 
-        axis = self.axis
+        plt.plot(self.axis, self.data())
-
-        plt.plot(axis, self.data(axis))
         plt.xlabel('x')
         plt.ylabel('f(x)')
         plt.autoscale(axis='x', tight=True)
@@ -442,32 +449,3 @@ class ProbabilityDensityFunction(object):
 
         return x0, incr, npts
 
-    def rearrange(self, other):
-        '''
-        Method rearrange takes another Probability Density Function and returns
-        a new axis with mid-point 0 and covering positive and negative range
-        of axis values, either containing the maximum value of both axis or
-        the sum of the maxima
-        :param other:
-        :return:
-        '''
-
-        x0 = self.x0
-        incr = self.incr
-        npts = self.npts
-
-
-        pdf_self = np.zeros(npts)
-        pdf_other = np.zeros(npts)
-
-        x = create_axis(x0, incr, npts)
-
-        sstart = find_nearest(x, self.x0)
-        s_end = sstart + self.data.size
-        ostart = find_nearest(x, other.x0)
-        o_end = ostart + other.data.size
-
-        pdf_self = self.broadcast(pdf_self, sstart, s_end, self.data)
-        pdf_other = self.broadcast(pdf_other, ostart, o_end, other.data)
-
-        return x0, incr, npts, pdf_self, pdf_other