Guest User

Untitled

a guest
Aug 19th, 2019
66
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