From a8b7eff561f4b849b392ef24ccdd37a39c393c70 Mon Sep 17 00:00:00 2001
From: sebastianp <sebastian.priebe@rub.de>
Date: Wed, 29 Jun 2016 15:29:58 +0200
Subject: [PATCH] [task] implementing new methods for pdf comparison.

---
 pylot/RELEASE-VERSION      |  2 +-
 pylot/core/io/phases.py    |  2 +-
 pylot/core/pick/compare.py | 80 +++++++++++++++++++++++++++++++++++++-
 pylot/core/util/pdf.py     | 23 +++++++++++
 4 files changed, 103 insertions(+), 4 deletions(-)

diff --git a/pylot/RELEASE-VERSION b/pylot/RELEASE-VERSION
index 671949c9..e760bf7c 100644
--- a/pylot/RELEASE-VERSION
+++ b/pylot/RELEASE-VERSION
@@ -1 +1 @@
-7c5a-dirty
+8714-dirty
diff --git a/pylot/core/io/phases.py b/pylot/core/io/phases.py
index e3c7e90a..60d4d180 100644
--- a/pylot/core/io/phases.py
+++ b/pylot/core/io/phases.py
@@ -180,7 +180,7 @@ def picksdict_from_picks(evt):
         try:
             onsets = picks[station]
         except KeyError as e:
-            print(e)
+            #print(e)
             onsets = {}
         mpp = pick.time
         spe = pick.time_errors.uncertainty
diff --git a/pylot/core/pick/compare.py b/pylot/core/pick/compare.py
index b39f8ef8..b2d93a1c 100644
--- a/pylot/core/pick/compare.py
+++ b/pylot/core/pick/compare.py
@@ -2,9 +2,9 @@
 # -*- coding: utf-8 -*-
 
 import copy
-
+import os
 import numpy as np
-
+import glob
 import matplotlib.pyplot as plt
 
 from obspy import read_events
@@ -336,3 +336,79 @@ class PDFDictionary(object):
             plt.setp([a.get_yticklabels() for a in axarr[:, 0]], visible=False)
 
         return fig
+
+
+class PDFstatistics(object):
+    def __init__(self, directory):
+        self.directory = directory
+        self.stations = {}
+        self.p_std = {}
+        self.s_std = {}
+        self.makeDirlist()
+        self.getData()
+        self.getStatistics()
+        #self.showData()
+
+    def makeDirlist(self, fn_pattern='*'):
+        self.dirlist = glob.glob1(self.directory, fn_pattern)
+        #print self.dirlist
+        self.evtdict = {}
+        for rd in self.dirlist:
+            #print rd
+            self.evtdict[rd] = glob.glob1((os.path.join(self.directory, rd)), '*.xml')
+            #print rd, self.evtdict[rd]
+
+    def getData(self):
+        arraylen = 0
+        for dir in self.dirlist:
+            for evt in self.evtdict[dir]:
+                self.stations[evt] = []
+                self.p_std[evt] = []
+                self.s_std[evt] = []
+                self.getPDFDict(dir, evt)
+                for station, pdfs in self.pdfdict.pdf_data.items():
+                    #print station, pdfs
+                    try:
+                        p_std = pdfs['P'].standard_deviation()
+                    except KeyError:
+                        p_std = 0
+                    try:
+                        s_std = pdfs['S'].standard_deviation()
+                    except KeyError:
+                        s_std = 0
+                    self.stations[evt].append(station)
+                    self.p_std[evt].append(p_std)
+                    self.s_std[evt].append(s_std)
+                    arraylen += 1
+        self.makeArray(arraylen)
+
+    def makeArray(self, arraylen):
+        index = 0
+        self.p_stdarray = np.zeros(arraylen)
+        self.s_stdarray = np.zeros(arraylen)
+        for evt in self.p_std.keys():
+            for n in range(len(self.p_std[evt])):
+                self.p_stdarray[index] = self.p_std[evt][n]
+                self.s_stdarray[index] = self.s_std[evt][n]
+                index += 1
+
+    def showData(self):
+        #figure(p, s) = plt.pyplot.subplots(2)
+        #p.hist(self.p_stdarray, Bins=100, range=[min(self.p_stdarray),max(self.p_stdarray)/256])
+        #s.hist(self.s_stdarray, Bins=100, range=[min(self.s_stdarray),max(self.s_stdarray)/256])
+        #p.set_title('Histogramm der P-Wellen-Standartabweichung')
+        #s.set_title('Histogramm der S-Wellen-Standartabweichung')
+        # plt.show()
+        pass
+
+
+    def getPDFDict(self, month, evt):
+        self.pdfdict = PDFDictionary.from_quakeml(os.path.join(self.directory,month,evt))
+
+    def getStatistics(self):
+        self.p_mean = self.p_stdarray.mean()
+        self.p_std_std = self.p_stdarray.std()
+        self.p_median = np.median(self.p_stdarray)
+        self.s_mean = self.s_stdarray.mean()
+        self.s_std_std = self.s_stdarray.std()
+        self.s_median = np.median(self.s_stdarray)
diff --git a/pylot/core/util/pdf.py b/pylot/core/util/pdf.py
index c2956e90..ed6ea70c 100644
--- a/pylot/core/util/pdf.py
+++ b/pylot/core/util/pdf.py
@@ -294,6 +294,29 @@ class ProbabilityDensityFunction(object):
             return None
         return self.data[find_nearest(self.axis, value)] * self.incr
 
+    def quantile(self, prob_value, eps=0.01):
+        l = self.axis[0]
+        r = self.axis[-1]
+        m = (r - l) / 2
+        diff = prob_value - self.prob_lt_val(m)
+        while abs(diff) > eps:
+            if diff > 0:
+                l = m
+            else:
+                r = m
+            m = (r - l) / 2
+            diff = prob_value - self.prob_lt_val(m)
+        return m
+
+
+        pass
+
+    def quantile_distance(self, prob_value):
+        ql = self.quantile(prob_value)
+        qu = self.quantile(1 - prob_value)
+
+        return qu - ql
+
     def plot(self, label=None):
         import matplotlib.pyplot as plt