From c7d8d0ecf572224beaacba56cf8c13b1186abba9 Mon Sep 17 00:00:00 2001
From: sebastianp <sebastian.priebe@rub.de>
Date: Wed, 13 Jul 2016 14:15:44 +0200
Subject: [PATCH] undid earlier changes in PDFStatistics

---
 pylot/core/pick/compare.py | 38 ++++++++++++++++++++++++++++----------
 pylot/core/util/pdf.py     | 16 ++++++++++++----
 2 files changed, 40 insertions(+), 14 deletions(-)

diff --git a/pylot/core/pick/compare.py b/pylot/core/pick/compare.py
index 37534660..0adaa388 100644
--- a/pylot/core/pick/compare.py
+++ b/pylot/core/pick/compare.py
@@ -350,13 +350,15 @@ class PDFstatistics(object):
         self.stations = {}
         self.p_std = {}
         self.s_std = {}
-        self.makeDirlist()
-        self.getData()
-        self.getStatistics()
-        self.arraylen = 0
-        self.Theta015 = []
-        self.Theta25 = []
-        self.Theta485 = []
+        self.theta015 = []
+        self.theta1 = []
+        self.theta2 = []
+        self.dirlist = list()
+        self.evtdict = {}
+        #self.makeDirlist()
+        #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/'
-    Insheim = PDFstatistics(rootdir)
\ No newline at end of file
+    rootdir = '/home/sebastianp/Data/Reassessment/Insheim/2012.10/'
+    Insheim = PDFstatistics(rootdir)
+    Insheim.dirlist = [rootdir]
+    Insheim.evtdict[rootdir] = ['e0002.294.12.xml']
+    Insheim.getData()
+    print Insheim.theta015
+    print Insheim.theta1
+    print Insheim.theta2
diff --git a/pylot/core/util/pdf.py b/pylot/core/util/pdf.py
index 9e1b0bee..443bf072 100644
--- a/pylot/core/util/pdf.py
+++ b/pylot/core/util/pdf.py
@@ -325,14 +325,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):
-        lim = np.arange(limits[0], limits[1], self.incr)
+    def prob_limits(self, limits, oversampling=1.):
+        sampling = self.incr / oversampling
+        lim = np.arange(limits[0], limits[1], sampling)
         data = [self.data(t) for t in 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 +353,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,6 +368,13 @@ 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