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. OK, I Understand
 
Top