From c7d8d0ecf572224beaacba56cf8c13b1186abba9 Mon Sep 17 00:00:00 2001 From: sebastianp Date: Wed, 13 Jul 2016 14:15:44 +0200 Subject: [PATCH 1/2] 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 From 8433767b222729d5ad6f05eb632d8575cf4028f6 Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Wed, 13 Jul 2016 14:37:12 +0200 Subject: [PATCH 2/2] [change] consequently use new pdf evaluation concept throughout the entire code --- pylot/core/pick/compare.py | 4 +-- pylot/core/util/pdf.py | 54 +++++++++----------------------------- 2 files changed, 14 insertions(+), 44 deletions(-) diff --git a/pylot/core/pick/compare.py b/pylot/core/pick/compare.py index 0adaa388..73122008 100644 --- a/pylot/core/pick/compare.py +++ b/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), diff --git a/pylot/core/util/pdf.py b/pylot/core/util/pdf.py index 443bf072..7c26b2be 100644 --- a/pylot/core/util/pdf.py +++ b/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_other = np.array([other.data(x) for x in axis]) + pdf_self = np.array(self.data(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_other = np.array([other.data(x) for x in axis]) + pdf_self = np.array(self.data(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, - self.npts)]) + return str([self.data()]) @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 @@ -328,7 +329,7 @@ class ProbabilityDensityFunction(object): 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] + 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]) @@ -370,7 +371,7 @@ class ProbabilityDensityFunction(object): def qtile_dist_quot(self,x): - if x <= 0 or x >= 0.5: + if x < 0 or x > 0.5: raise ValueError('Value out of range.') return self.quantile_distance(0.5-x)/self.quantile_distance(x) @@ -378,9 +379,7 @@ class ProbabilityDensityFunction(object): def plot(self, label=None): import matplotlib.pyplot as plt - axis = self.axis - - plt.plot(axis, self.data(axis)) + plt.plot(self.axis, self.data()) plt.xlabel('x') plt.ylabel('f(x)') plt.autoscale(axis='x', tight=True) @@ -450,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