From a475b366d47c23a329698db3ce7e41a40bf21f27 Mon Sep 17 00:00:00 2001
From: Sebastianw Wehling-Benatelli <sebastian.wehling@rub.de>
Date: Thu, 31 Mar 2016 14:29:29 +0200
Subject: [PATCH] [adresses #195] comparison dictionaries can now be compiled
 from QuakeML files

---
 pylot/core/pick/compare.py | 11 ++++++++++-
 pylot/core/util/pdf.py     | 13 +++++++++++--
 2 files changed, 21 insertions(+), 3 deletions(-)

diff --git a/pylot/core/pick/compare.py b/pylot/core/pick/compare.py
index e996bffe..4e1025ea 100644
--- a/pylot/core/pick/compare.py
+++ b/pylot/core/pick/compare.py
@@ -49,4 +49,13 @@ def compare_picksets(a, b):
 
     compare_pdfs = dict()
 
-    return compare_pdfs
\ No newline at end of file
+    for station, phases in pdf_a.items():
+        if station in pdf_b.keys():
+            compare_pdf = dict()
+            for phase in phases:
+                if phase in pdf_b[station].keys():
+                    compare_pdf[phase] = phases[phase] - pdf_b[station][phase]
+            if compare_pdf is not None:
+                compare_pdfs[station] = compare_pdf
+
+    return compare_pdfs
diff --git a/pylot/core/util/pdf.py b/pylot/core/util/pdf.py
index c9e79e54..05104ee1 100644
--- a/pylot/core/util/pdf.py
+++ b/pylot/core/util/pdf.py
@@ -270,12 +270,21 @@ class ProbabilityDensityFunction(object):
     def prob_lt_val(self, value):
         if value <= self.axis[0] or value > self.axis[-1]:
             raise ValueError('value out of bounds: {0}'.format(value))
-        return self.data[self.axis <= value].sum() * self.incr
+        return self.prob_limits((self.axis[0], value))
 
     def prob_gt_val(self, value):
         if value < self.axis[0] or value >= self.axis[-1]:
             raise ValueError('value out of bounds: {0}'.format(value))
-        return self.data[self.axis >= value].sum() * self.incr
+        return self.prob_limits((value, self.axis[-1]))
+
+    def prob_limits(self, limits):
+        lim_ind = np.logical_and(limits[0] <= self.axis, self.axis <= limits[1])
+        data = self.data[lim_ind]
+        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
 
     def prob_val(self, value):
         if not (self.axis[0] <= value <= self.axis[-1]):