• API
• FAQ
• Tools
• Archive
SHARE
TWEET

Untitled

a guest Aug 19th, 2019 65 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. #!/usr/bin/env python2
2. # -*- coding: utf-8 -*-
3. """
4. Created on Wed Mar 15 11:16:40 2017
5.
6. @author: ltlustos
7. """
8. from scipy.optimize import curve_fit
9. import numpy as np
10. import pylab as pl
11. from scipy.signal import argrelmax, argrelmin
12. import scipy.ndimage  # @UnresolvedImport
13.
14.
15. def mov_median(data, span):
16.     return scipy.ndimage.median_filter(data, span)
17.
18.
19. def mov_mean(data, span):
20.     return np.convolve(data, np.ones((span,))/span, mode='same')
21.
22.
23. def moment1(x, y=None, nstd=None, warning=True):
24.     if y is None:
25.         y = x
26.         x = np.arange(len(y))
27.     if y.min() < 0:
28.         if warning:
29.             print('Warniong: negative weights: [%g, %g]' % (y.min(), y.max()))
30. #            raise Exception('neg weight')
31.     if x.dtype in (np.uint8, np.uint16, np.uint32, np.uint64):
32.         x = x.astype(np.double)
33.     if y.dtype in (np.uint8, np.uint16, np.uint32, np.uint64):
34.         y = y.astype(np.double)
35. #    offs = y.min()
36. #    y -= offs
37.     m0 = np.double(np.sum(y))
38.     m1 = np.sum(y * x)
39.     m2 = np.sum(y * x ** 2)
40.     amean = m1 / m0
41.     astd = np.sqrt((m2 - m1 ** 2 / m0) / m0)
42.     if nstd is not None:
43.         idxx = np.where(np.abs(x-amean) < nstd*astd)
44.         amean, astd = moment1(x[idxx], y[idxx])
45.     return [amean, astd]
46.
47.
48. def local_min(x, n=3):
49.     return argrelmin(x, order=n)
50.
51.
52. def local_max(x, n=3):
53.     return argrelmax(x, order=n)
54.
55.
56. def find_local_max(x, y, mincount=1, sigma_est=6):
57.     y /= y.max()
58.     y = mov_mean(mov_mean(y ** 3, sigma_est), sigma_est)
59.
60.     x = x[0:-3]
61.     y = y[0:-3]
62.
63.     idxx = np.where(y < mincount)[0]
64.     y[idxx] = 0
65.     return x[local_max(y)], y[local_max(y)]
66.
67. ###############################################################
68.
69.
70. class GOF(object):
71.
72.     def __init__(self, y, fy, weights=None):
73.         self.__y = y
74.         self.__fy = fy
75.         if weights is not None:
76.             weights /= np.sum(weights)
77.         else:
78.             weights = 1.0
79.
80.         ybar = np.mean(self.__y)
81.         # total sum of squares
82.         self.__SST = np.sum(weights * (self.__y - ybar) ** 2)
83.         # estimated sum of squares
84.         self.__SSE = np.sum(weights * (self.__fy - self.__y) ** 2)
85.
86.         if self.__SST != 0:
87.             self.__R2 = 1 - self.__SSE / self.__SST
88.         elif (self.__SST == 0) & (self.__SSE == 0):
89.             self.__R2 = np.NaN
90.         else:
91.             if np.sqrt(np.abs(self.__SSE)) < np.sqrt(np.finfo(np.double).eps) * np.mean(abs(self.__y)):
92.                 self.__R2 = np.NaN
93.             else:
94.                 self.__R2 = np.Inf
95.         if np.isnan(self.__R2):
96.             self.__R2 = -1.0
97.
98.     @property
99.     def R2(self):
100.         return self.__R2
101.
102.     @property
103.     def ESS(self):
104.         return self.__SSE
105.
106.     @property
107.     def TSS(self):
108.         return self.__SST
109.
110. ###############################################################
111.
112.
113. class Fit(object):  # dont foreget the "object" root class
114.     def __init__(self, func):
115.         self.__func__ = func
116.         self._coeff_ = None
117.         self._var_matrix_ = None
118.
119.     def __eval_func__(self, x, *p):
120.         return eval(self.__func__)
121.
122.     def eval_func(self, x, p):
123.         return self.__eval_func__(x, *p)
124.
125.     def fit(self, x, y, p0, weights=None):
126.         # from scipy.optimize import curve_fit
127.         self.__x__ = x
128.         self.__y__ = y
129.         self.__weights__ = weights
130.         if weights is not None:
131.             self.__weights__ = 1.0 / weights
132.         self._coeff_, self._var_matrix_ = curve_fit(self.__eval_func__, x, y, p0=p0, sigma=self.__weights__,
133.                                                     absolute_sigma=weights is not None)
134.         self.GOF = GOF(y, self.__eval_func__(x, *self._coeff_), weights=weights)
135.         return self._coeff_, self._var_matrix_, self.GOF
136.
137.     def refit_weighted(self, n_rep=1):
138.         for i in range(n_rep):
139.             self.fit(self.__x__, self.__y__, p0=self.p, weights=self.eval_func(self.__x__, self.p)**2)
140.         return self._coeff_, self._var_matrix_, self.GOF
141.
142.     def plot(self, x=None, p=None, legend=True):
143.         if x is None:
144.             x = np.linspace(self.__x__.min(), self.__x__.max(), 5*len(self.__x__))
145.             # decrease spacing between fit point by factor 5
146.         if p is None:
147.             p = self._coeff_
148.         if legend:
149.             pl.plot(x, self.__eval_func__(x, *p), label='%s\n%s, R2=%g' % (self.__func__, p, self.GOF.R2))
150.             pl.plot(self.__x__, self.__y__, 'o', label='data')
151.         else:
152.             pl.plot(x, self.__eval_func__(x, *p))
153.             pl.plot(self.__x__, self.__y__, 'o')
154.         if legend:
155.             pl.legend(loc='1', fontsize='small')
156.
157.     @property
158.     def p(self):
159.         return self._coeff_
160. ###############################################################
161.
162.
163. def gauss(x, p):
164.     return p[0] * np.exp(-(x-p[1])**2/(2.*p[2]**2))
165.
166.
167. def erf(x, p):
168.     import scipy.special
169.     return p[0] * scipy.special.erf((x - p[1]) / p[2] / np.sqrt(2.0))
170.
171.
172. def erf_r(x, p0, p1, p2):
173.     import scipy.special
174.     return p0 * (1 - scipy.special.erf((x - p1) / p2 / np.sqrt(2.0))) / 2.0
175.
176.
177. def erf_l(x, p0, p1, p2):
178.     import scipy.special
179.     return p0 * (1 + scipy.special.erf((x - p1) / p2 / np.sqrt(2.0))) / 2.0
180.
181.
182. def tanh_r(x, p0, p1, p2):
183.     return p0 * (1 - np.tanh((x - p1) / p2 / np.sqrt(2.0))) / 2.0
184.
185.
186. def tanh_l(x, p0, p1, p2):
187.     return p0 * (1 + np.tanh((x - p1) / p2 / np.sqrt(2.0))) / 2.0
188.
189.
190. class FitGauss(Fit):
191.     def __init__(self):
192.         # p[0] : amplitude
193.         # p[1] : mu
194.         # p[2]: sigma
195.         super(FitGauss, self).__init__('gauss(x, p)')
196.
197.     def fit(self, x, y, p0=None, weights=None):
198.         if p0 is None:
199.             mu, std = moment1(x, y)
200.             p0 = [y.max(), mu, std]
201.         super(FitGauss, self).fit(x, y, p0, weights=weights)
202.         self._coeff_[2] = np.abs(self._coeff_[2])
203.         return self._coeff_, self._var_matrix_, self.GOF
204.
205.
206. class FitGauss_Erf_right(Fit):
207.     def __init__(self):
208.         # p[0] : amplitude gauss
209.         # p[1] : mu
210.         # p[2]: sigma
211.         # p[3] : amplitude erf
212.         super(FitGauss_Erf_right, self).__init__('gauss(x, p[0:3]) + erf_r(x, p[3], p[1], p[2])')
213.
214.     def fit(self, x, y, p0=None, weights=None):
215.         if p0 is None:
216.             mu, std = moment1(x, y)
217.             p0 = [y.max(), mu, std, y.max()/10, mu, std]
218.         super(FitGauss_Erf_right, self).fit(x, y, p0, weights=weights)
219.         self._coeff_[2] = np.abs(self._coeff_[2])
220.         return self._coeff_, self._var_matrix_, self.GOF
221.
222.     def plot(self, x=None, p=None):
223.         super(FitGauss_Erf_right, self).plot(x, p)
224.         import pylab as pl
225.         if x is None:
226.             x = np.linspace(self.__x__.min(), self.__x__.max(), 5*len(self.__x__))
227.             # decrease spacing between fit point by factor 5
228.         if p is None:
229.             p = self._coeff_
230.         pl.plot(x, erf_r(x, p[3], p[1], p[2]), ':', zorder=0, color='0.75')
231.         pl.plot(x, gauss(x, p[0:3]), ':', zorder=0, color='0.75')
232.
233.
234. class FitGauss_tanh_right(Fit):
235.     def __init__(self):
236.         # p[0] : amplitude gauss
237.         # p[1] : mu
238.         # p[2]: sigma
239.         # p[3] : amplitude erf
240.         super(FitGauss_tanh_right, self).__init__('gauss(x, p[0:3]) + tanh_r(x, p[3], p[1], p[2])')
241.
242.     def fit(self, x, y, p0=None, weights=None):
243.         if p0 is None:
244.             mu, std = moment1(x, y)
245.             p0 = [y.max()-y[0], mu, std, y[0], mu, std]
246.             print('p0=', p0)
247.         super(FitGauss_tanh_right, self).fit(x, y, p0, weights=weights)
248.         self._coeff_[2] = np.abs(self._coeff_[2])
249.         return self._coeff_, self._var_matrix_, self.GOF
250.
251.     def plot(self, x=None, p=None):
252.         super(FitGauss_tanh_right, self).plot(x, p)
253.         import pylab as pl
254.         if x is None:
255.             x = np.linspace(self.__x__.min(), self.__x__.max(), 5*len(self.__x__))
256.             # decrease spacing between fit point by factor 5
257.         if p is None:
258.             p = self._coeff_
259.         pl.plot(x, tanh_r(x, p[3], p[1], p[2]), ':', zorder=0, color='0.75')
260.         pl.plot(x, gauss(x, p[0:3]), ':', zorder=0, color='0.75')
261.
262.
263. class Fit2Gauss_tanh_right(Fit):
264.     def __init__(self):
265.         # p[0] : amplitude gauss
266.         # p[1] : mu
267.         # p[2]: sigma
268.         # p[3] : amplitude erf
269.         super(Fit2Gauss_tanh_right, self).__init__('gauss(x, p[0:3]) + tanh_r(x, p[3], p[1], p[2]) + gauss(x, p[4:])')
270.
271.     def fit(self, x, y, p0=None, weights=None):
272.         if p0 is None:
273.             mu, std = moment1(x, y)
274.             p0 = [y.max(), mu, std, y.max()/10, mu, std]
275.         super(Fit2Gauss_tanh_right, self).fit(x, y, p0, weights=weights)
276.         self._coeff_[2] = np.abs(self._coeff_[2])
277.         return self._coeff_, self._var_matrix_, self.GOF
278.
279.     def plot(self, x=None, p=None):
280.         super(Fit2Gauss_tanh_right, self).plot(x, p)
281.         import pylab as pl
282.         if x is None:
283.             x = np.linspace(self.__x__.min(), self.__x__.max(), 5*len(self.__x__))
284.             # decrease spacing between fit point by factor 5
285.         if p is None:
286.             p = self._coeff_
287.         pl.plot(x, tanh_r(x, p[3], p[1], p[2]), ':', zorder=0, color='0.75')
288.         pl.plot(x, gauss(x, p[0:3]), ':', zorder=0, color='0.75')
289.         pl.plot(x, gauss(x, p[4:]), ':', zorder=0, color='0.75')
290.
291.
292. class Fit2Gauss(Fit):
293.     def __init__(self):
294.         # p[0] : amplitude
295.         # p[1] : mu
296.         # p[2]: sigma
297.         super(Fit2Gauss, self).__init__('gauss(x, p[0:3]) + gauss(x, p[3:])')
298.
299.     def fit(self, x, y, p0=None, weights=None):
300.         if p0 is None:
301.             # to be done
302.             mu, std = moment1(x, y)
303.             p0 = [y.max()/2, mu/2, std/2, y.max(), mu, std]
304.         super(Fit2Gauss, self).fit(x, y, p0, weights=weights)
305.         self._coeff_[2::3] = np.abs(self._coeff_[2::3])
306.         return self._coeff_, self._var_matrix_, self.GOF
307.
308.
309. class Fit3Gauss(Fit):
310.     def __init__(self):
311.         # p[0] : amplitude
312.         # p[1] : mu
313.         # p[2]: sigma
314.         super(Fit3Gauss, self).__init__('gauss(x, p[0:3]) + gauss(x, p[3:6]) + gauss(x, p[6:9])')
315.
316.     def fit(self, x, y, p0=None, weights=None):
317.         if p0 is None:
318.             # to be done
319.             mu, std = moment1(x, y)
320.             p0 = [y.max()/2, mu/2, std/2, y.max(), mu, std]
321.         super(Fit3Gauss, self).fit(x, y, p0, weights=weights)
322.         self._coeff_[2::3] = np.abs(self._coeff_[2::3])
323.         return self._coeff_, self._var_matrix_, self.GOF
324.
325.
326. class Fit4Gauss(Fit):
327.     def __init__(self):
328.         # p[0] : amplitude
329.         # p[1] : mu
330.         # p[2]: sigma
331.         super(Fit4Gauss, self).__init__('gauss(x, p[0:3]) + gauss(x, p[3:6]) + gauss(x, p[6:9]) + gauss(x, p[9:12])')
332.
333.     def fit(self, x, y, p0=None, weights=None):
334.         if p0 is None:
335.             # to be done
336.             mu, std = moment1(x, y)
337.             p0 = [y.max()/2, mu/2, std/2, y.max(), mu, std]
338.         super(Fit3Gauss, self).fit(x, y, p0, weights=weights)
339.         self._coeff_[2::3] = np.abs(self._coeff_[2::3])
340.         return self._coeff_, self._var_matrix_, self.GOF
341.
342. ###############################################################
343.
344.
345. class Fit_ToT_keV(Fit):
346.     def __init__(self):
347.         super(Fit_ToT_keV, self).__init__('(p[0] * x + p[1] - p[2] / ( x - p[3])) * np.array((x > '
348.                                           '((p[0]*p[3] - p[1] - np.sqrt(p[0]**2*p[3]**2 + 2*p[0]*p[1]*p[3] '
349.                                           '+ 4*p[0]*p[2] + p[1]**2))/(2*p[0]))))')
350.
351.     def __eval_func__(self, x, *p):
352.         # x0 = (p[0]*p[3] - p[1] + np.sqrt(p[0]**2*p[3]**2 + 2*p[0]*p[1]*p[3] + 4*p[0]*p[2] + p[1]**2))/(2*p[0])
353.         # print x0
354.         y = p[0] * x + p[1] - p[2] / (x - p[3])
355.         return y * np.array(y > 0)
356.
357.     def fit(self, x, y, p0=None, weights=None):
358.         if p0 is None:
359.             A = np.polyfit(x, y, 1)
360.             p0 = [A[1], A[0], 10, 3]
361.         return super(Fit_ToT_keV, self).fit(x, y, p0, weights=weights)
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.

Top