Update example.py, requirements.txt, Model.py, LossFunctions.py,...

Update example.py, requirements.txt, Model.py, LossFunctions.py, HelperFunctions.py, AdditionalPDFs.py, PeakOTron.py

AdditionalPDFs.py

+import numpy as np
+from scipy.stats import rv_discrete, rv_continuous, uniform
+import scipy.special as sc
+import matplotlib.pyplot as plt
+
+
+from scipy.stats._distn_infrastructure import (
+        rv_discrete, _ncx2_pdf, _ncx2_cdf, get_distribution_names)
+
+
+class gpd_gen(rv_discrete):
+    def _argcheck(self, mu, lbda):
+        return mu >= 0.0 and lbda >= 0.0 and lbda <= 1.0
+
+    def _rvs(self, mu, lbda):
+        population = np.asarray(
+            self._random_state.poisson(mu, self._size)
+        )
+        if population.shape == ():
+            population = population.reshape(-1)
+        offspring = population.copy()
+        while np.any(offspring > 0):
+            # probability dists are NOT ufuncs
+            # print("offspring", offspring)
+            offspring[:] = [
+                self._random_state.poisson(m)
+                for m in lbda*offspring
+            ]
+            population += offspring
+        return population
+
+    def _pmf(self, k, mu, lbda):
+        return np.exp(self._logpmf(k, mu, lbda))
+
+    def _logpmf(self, k, mu, lbda):
+        mu_pls_klmb = mu + lbda*k
+        return np.log(mu) + sc.xlogy(k-1, mu_pls_klmb) - mu_pls_klmb - sc.gammaln(k+1)
+
+    def _munp(self, n, mu, lbda):
+        if n == 1:
+            return mu/(1-lbda)
+        elif n == 2:
+            return (mu/(1-lbda))**2+mu/(1-lbda)**3
+
+
+gpoisson = gpd_gen(name='gpoisson')
+
+    
+class borel_gen(rv_discrete):
+    def _argcheck(self, mu):
+        return ((mu > 0) & (mu<1))
+
+    def _logpmf(self, k, mu):
+        n = k+1
+        Pk = sc.xlogy(n-1, mu*n) - sc.gammaln(n + 1) - mu*n
+        return Pk
+
+    def _pmf(self, k, mu):
+        return np.exp(self._logpmf(k, mu))
+
+#     def _rvs(self, mu, size=None, random_state=None):
+#         u = uniform.rvs(loc=0, scale = 1, size=size)
+#         cum = np.cumsum([self._pmf(_k, mu) for _k in range(0, 100)])
+#         print(cum)
+#         rnd = [ np.argmax( cum>=_u ) for _u in u ]
+#         return rnd
+    
+    def _rvs(self, mu, size=None, random_state=None, epsilon=1e-4):
+        _u = uniform.rvs(loc=0, scale = 1-epsilon, size=size)
+        _sum = 0
+        _k=0
+        _elem = []
+        _max_u = np.max(_u)
+        
+        while(_sum<_max_u):
+            _pmf = self._pmf(_k, mu)
+            _elem.append(_pmf)
+            _sum+=_pmf
+            _k+=1
+            
+        _cum = np.cumsum(_elem)
+        _rnd = [ np.argmax( _cum>=__u ) for __u in _u ]
+        
+        return _rnd
+
+
+    def _stats(self, mu):
+        _mu = 1/(1-mu)
+        _var = mu/(1-mu)**3
+        tmp = np.asarray(mu)
+        mu_nonzero = ((tmp > 0) & (tmp<1))
+        #g1 and g2: Lagrangian Probability Distributions, 978-0-8176-4365-2, page 159
+        g1 = scipy._lib._util._lazywhere(mu_nonzero, (tmp,), lambda x: (1+2*x)/scipy.sqrt(x*(1-x)), np.inf)
+        g2 = scipy._lib._util._lazywhere(mu_nonzero, (tmp,), lambda x: 3 + (1 + 8*x+6*x**2)/(x*(1-x)), np.inf)
+        return _mu, _var, g1, g2
+
+
+borel= borel_gen(name='borel')
+    
+  
+    
+class erlang_gen(rv_discrete):
+    
+    
+    
+    def _pdf(self, x, a):
+        # gamma.pdf(x, a) = x**(a-1) * exp(-x) / gamma(a)
+        return np.exp(self._logpdf(x, a))
+
+    def _logpdf(self, k, mu, nu):
+        return sc.xlogy(a-1.0, x) - x - sc.gammaln(a)
+
+   
+    
+    
+  
+
+#     def _rvs(self, mu, nu, size=None, random_state=None):
+#         u = scipy.stats.uniform.rvs(loc=0, scale = 1, size=size)
+#         cum = np.cumsum([self._pmf(_k, mu, nu) for _k in range(0, 100)])
+#         rnd = [ np.argmax( cum>=_u ) for _u in u ]
+#         return rnd
+
+pairs = list(globals().items())
+_distn_names, _distn_gen_names = get_distribution_names(pairs, rv_discrete)
+
+__all__ = _distn_names + _distn_gen_names

HelperFunctions.py

+from iminuit import Minuit
+from iminuit.util import describe
+from scipy.interpolate import interp1d
+from astropy.stats import bootstrap
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+#HELPER FUNCTIONS
+
+class FakeFuncCode:
+    def __init__(self, f, prmt=None, dock=0, append=None):
+        #f can either be tuple or function object
+        self.co_varnames = describe(f)
+        self.co_argcount = len(self.co_varnames)
+        self.co_argcount -= dock
+        self.co_varnames = self.co_varnames[dock:]
+
+        if prmt is not None:  #rename parameters from the front
+            for i, p in enumerate(prmt):
+                self.co_varnames[i] = p
+
+        if isinstance(append, str): append = [append]
+
+        if append is not None:
+            old_count = self.co_argcount
+            self.co_argcount += len(append)
+            self.co_varnames = tuple(
+                list(self.co_varnames[:old_count]) +
+                append +
+                list(self.co_varnames[old_count:]))
+
+
+
+            
+def LatexFormat(value, scirange=[0.01,1000]):
+    if(np.abs(value)>scirange[0] and np.abs(value)<scirange[1]):
+        float_str = r"${:3.3f}$".format(value)
+    else:
+        try:
+            float_str = "{:3.3E}".format(value)
+            base, exponent = float_str.split("E")
+            float_str = r"${0} \times 10^{{{1}}}$".format(base, int(exponent))
+        except:
+            float_str=str(value)
+    return float_str
+
+
+def FormatExponent(ax, axis='y'):
+
+    # Change the ticklabel format to scientific format
+    ax.ticklabel_format(axis=axis, style='sci', scilimits=(-2, 2))
+
+    # Get the appropriate axis
+    if axis == 'y':
+        ax_axis = ax.yaxis
+        x_pos = 0.0
+        y_pos = 1.0
+        horizontalalignment='left'
+        verticalalignment='bottom'
+    else:
+        ax_axis = ax.xaxis
+        x_pos = 1.0
+        y_pos = -0.05
+        horizontalalignment='right'
+        verticalalignment='top'
+
+    # Run plt.tight_layout() because otherwise the offset text doesn't update
+    plt.tight_layout()
+    ##### THIS IS A BUG 
+    ##### Well, at least it's sub-optimal because you might not
+    ##### want to use tight_layout(). If anyone has a better way of 
+    ##### ensuring the offset text is updated appropriately
+    ##### please comment!
+
+    # Get the offset value
+    offset = ax_axis.get_offset_text().get_text()
+
+    if len(offset) > 0:
+        # Get that exponent value and change it into latex format
+        minus_sign = u'\u2212'
+        expo = np.float(offset.replace(minus_sign, '-').split('e')[-1])
+        offset_text = r'x$\mathregular{10^{%d}}$' %expo
+
+        # Turn off the offset text that's calculated automatically
+        ax_axis.offsetText.set_visible(False)
+
+        # Add in a text box at the top of the y axis
+        ax.text(x_pos, y_pos, offset_text, transform=ax.transAxes,
+               horizontalalignment=horizontalalignment,
+               verticalalignment=verticalalignment)
+    return ax
+
+
+def Bootstrap(data, statistic, n_bootstrap, alpha=0.95):
+    if not (0 < alpha < 1):
+        raise ValueError("confidence level must be in (0, 1)")
+
+
+    if len(data) < 1:
+        raise ValueError("data must contain at least one measurement.")
+
+
+    boot_stat = bootstrap(data, n_bootstrap, bootfunc=statistic)
+
+    stat_data = statistic(data)
+    mean_stat = np.mean(boot_stat)
+    est_stat = 2*stat_data - mean_stat
+    std_err = np.std(boot_stat)
+    z_score = np.sqrt(2.0)*sc.erfinv(alpha)
+    conf_interval = est_stat + z_score*np.array((-std_err, std_err))
+
+
+    return est_stat, std_err, conf_interval  
+
+
+def Linear(x, m, c):
+    return m*x + c
+
+def HistMoment(counts, bin_centres, i, mu):
+    return np.sum(counts*(bin_centres-mu)**i)
+
+def HistMean(counts, bin_centres):
+    N = np.sum(counts)
+    mu = HistMoment(counts, bin_centres, 1, 0)/N
+    return mu
+    
+def HistStd(counts, bin_centres):
+    N = np.sum(counts)
+    mu = HistMean(counts, bin_centres)  
+    var = HistMoment(counts, bin_centres, 2, mu)/(N-1)
+    sigma = np.sqrt(var)
+    return sigma
+
+def HistSkew(counts, bin_centres):
+    N = np.sum(counts)
+    mu = HistMean(counts, bin_centres)
+    gamma = (N*np.sqrt((N-1))/(N-2))*HistMoment(counts, bin_centres, 3, mu)/(HistMoment(counts, bin_centres, 2, mu)**(1.5))
+    return gamma
+
+
+
+
+def GP_lbda(mu, sigma, gamma):
+    lbda = 0.5*(((mu*gamma)/(sigma))- 1)
+    return lbda
+
+def GP_gain(mu, sigma, gamma):
+    lbda = GP_lbda(mu, sigma, gamma)
+    gain = (sigma**2/mu)*((1 - lbda)**2)
+    return gain
+
+def GP_mu(mu, sigma, gamma):
+    lbda = GP_lbda(mu, sigma, gamma)
+    mu_gp = (1/(1-lbda))*(mu**2/sigma**2)
+    return mu_gp
+
+
+
+
+
+
+    
+
+

LossFunctions.py

+import numpy as np
+from scipy.interpolate import interp1d
+from iminuit.util import describe
+from Model import epsilon
+from HelperFunctions import FakeFuncCode
+import matplotlib.pyplot as plt
+ 
+
+           
+                                              
+class BinnedLH:
+    
+
+    def __init__(self, f, bin_centres, counts, bw):
+        self.f = f
+        self.x = bin_centres
+        self.dx = bw
+        self.counts = counts
+        self.N = np.sum(counts)
+        
+        self.last_arg = None
+        self.func_code = FakeFuncCode(f, dock=True)
+        self.n_calls=0
+        self.eps = epsilon()
+        
+
+    def __call__(self, *arg):
+        self.last_arg = arg
+                
+        
+        y_hat = self.f(self.x, *arg)
+        
+        y_hat = np.nan_to_num(y_hat, nan=self.eps, posinf=self.eps, neginf=self.eps)
+        y_hat = np.where(y_hat<self.eps, self.eps, y_hat)
+        
+        E = y_hat*self.N*self.dx
+        h = self.counts
+        mask = (h>0)
+        E = E[mask]
+        h = h[mask]
+        
+        nlogL = -np.sum(h*(np.log(E) - np.log(h)) + (h-E))
+
+        self.n_calls+=1
+        
+        
+        return nlogL
+    
+    
+    
+    
+class Chi2Regression:
+    
+
+    def __init__(self, f, x, y, y_err, epsilon=1.35):
+        self.f = f
+        self.x = x
+        self.y = y
+        self.y_err = y_err
+        self.eps = np.finfo(np.float64).eps * 10
+        self.y_err[self.y_err<self.eps] = self.eps
+        self.last_arg = None
+        self.func_code = FakeFuncCode(f, dock=True)
+        self.ndof = len(self.y) - (self.func_code.co_argcount - 1)
+      
+
+    def __call__(self, *arg):
+        
+        self.last_arg = arg
+        
+        y_hat = self.f(self.x, *arg) 
+        
+        loss = ((y_hat - self.y)/(self.y_err))**2
+        
+        return np.sum(loss)
+    
+    
+    
+    
+class HuberRegression:
+    
+
+    def __init__(self, f, x, y, y_err, delta=1.345):
+        self.f = f
+        self.x = x
+        self.y = y
+        self.y_err = y_err
+        self.delta = delta
+        self.eps = np.finfo(np.float64).eps * 10
+        self.y_err[self.y_err<self.eps] = self.eps
+        self.last_arg = None
+        self.func_code = FakeFuncCode(f, dock=True)
+        self.ndof = len(self.y) - (self.func_code.co_argcount - 1)
+        
+
+    def __call__(self, *arg):
+        
+        self.last_arg = arg
+        
+        a = abs((self.y - self.f(self.x, *arg))/self.y_err)
+        cond_flag = (a>self.delta)
+        
+        loss = np.sum((~cond_flag) * (0.5 * a ** 2) - (cond_flag) * self.delta * (0.5 * self.delta - a), -1)
+        
+        return np.sum(loss)
\ No newline at end of file

example.py

+
+import os
+import sys
+import numpy as np
+import pandas as pd
+from PeakOTron import PeakOTron
+from joblib import dump
+import time
+
+print("--------------------")
+print('EXAMPLE SIPM CALIBRATION RUN')
+print("--------------------")
+
+
+out_dict = {}
+files_to_fit = []
+
+## Find all histograms in directory 
+for root, dirs, files in os.walk("./data/hamamatsu_pcb6/Light"):
+     
+        for file in files:
+            
+            if file.endswith(".txt"):
+                files_to_fit.append([file, os.path.join(root, file)])
+                
+
+## Print files.
+print("Files to fit:")
+for i, (file, _) in enumerate(files_to_fit):
+    print('File {0}: {1}'.format(i, file))
+    
+    
+    
+SiPM = "PM1125NS_SBO"
+
+
+## Loop thorough files
+for i, (file, path) in enumerate(files_to_fit):
+    items = file.split('_')
+
+    V = float(items[2].replace('V', '').replace('p', '.'))
+
+    print("\n\n")
+    print("===============================================================")
+    print("FIT {:d} - {:s}".format(i, file))
+    print("===============================================================")
+    print("\n\n")
+        
+
+    ## Load files. 
+    data = np.loadtxt(path, skiprows=8)
+
+    ## Create a PeakOTron Fit Object. 
+    f_data = PeakOTron(verbose=True)
+
+    ## Perform fit. 
+    f_data.Fit(data, 
+          tau=21.953,  #SLOW PULSE COMPONENT TIME CONSTANT (ns)
+          t_gate=104, #GATE LENGTH (ns)
+          t_0 = 64, #INTEGRATION TIME BEFORE GATE (ns)               
+          tau_R=21.953
+              
+    ) #BINNING RULE "knuth", "freedman", "scott" - use bw= #### to override. it is still required for DCR calculation. 
+
+    f_data.PlotFit(xlabel="ADC", display=False, save_directory="./Results/{0}_fit.png".format(os.path.splitext(file)[0]))
+    
+
+    
+    dump(f_data, "./Results/{0}.joblib".format(os.path.splitext(file)[0]))
+    
+    
+    
+    fit_out = {}
+    fit_val, fit_err = f_data.GetFitResults()
+    for key, val in fit_val.items():
+        print("{:s} : {:3.3E}".format(key, val))
+    
+    fit_out["SiPM"] = SiPM
+    fit_out["V"] = V
+
+    for key, value in fit_err.items():
+        fit_out["d_{:s}".format(key)] = value
+
+    fit_out.update(fit_val)
+    out_dict.update()
+    if out_dict == {}:
+        for key in fit_out.keys():
+            out_dict[key] = []
+
+    for key in fit_out.keys():
+        out_dict[key].append(fit_out[key])
+
+    print("===============================================================")
+    print("\n\n")
+        
+
+
+
+df = pd.DataFrame.from_dict(out_dict)
+df.to_csv("./fit_results_{:s}.csv".format(SiPM))

requirements.txt

+astropy==3.0.2
+iminuit==2.17.0
+joblib==0.13.2
+matplotlib==3.1.3
+numba==0.45.1
+numpy==1.17.2
+pandas==1.1.3
+scipy==1.3.1