diff options
Diffstat (limited to 'plot_llh')
| -rw-r--r-- | plot_llh/LVCPT.py | 431 | ||||
| -rw-r--r-- | plot_llh/MinimalTools.py | 166 | ||||
| -rw-r--r-- | plot_llh/PhysConst.py | 390 | ||||
| -rw-r--r-- | plot_llh/make_plots.py | 117 | ||||
| -rw-r--r-- | plot_llh/mcmc_scan.py | 174 |
5 files changed, 1278 insertions, 0 deletions
diff --git a/plot_llh/LVCPT.py b/plot_llh/LVCPT.py new file mode 100644 index 0000000..f5b02e6 --- /dev/null +++ b/plot_llh/LVCPT.py @@ -0,0 +1,431 @@ + +# coding: utf-8 + +## The Theory + +import numpy +import MinimalTools as MT +import PhysConst as PC +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.collections as mco +import matplotlib as mpl +import scipy.interpolate as interpolate +import scipy.integrate as integrate +import scipy as sp +from numpy import linalg as LA + +use_cython = False + +if use_cython: + import cython.cLVCPT as clv + +mpl.rc('font', family='serif', size=20) + +pc = PC.PhysicsConstants() + +degree = np.pi/180.0 +pc.th12 = 33.36*degree#33.48*degree +pc.th23 = 45.*degree#42.3*degree +pc.th13 = 8.66*degree#8.5*degree +pc.delta1 = 0.0#300.0*degree#306.*degree # perhaps better just set to 0. +pc.dm21sq = 7.5e-5 +pc.dm31sq = 2.47e-3#2.457e-3 +pc.Refresh() + +MT.calcU(pc) +DELTAM2 = MT.flavorM2(pc) + +def Hamiltonian(Enu, LVATERM = np.zeros((3,3), dtype=numpy.complex), + LVCTERM = np.zeros((3,3), dtype=numpy.complex)): + return DELTAM2/(2.0*Enu) + LVATERM + Enu*LVCTERM + +def OscProbFromMixingMatrix(alpha, beta, MixMatrix): + return sum([(np.absolute(MixMatrix[i][alpha])*np.absolute(MixMatrix[i][beta]))**2 for i in range(pc.numneu)] ) + #return sum([(np.absolute(MixMatrix[i][alpha]))**2*(np.absolute(MixMatrix[i][beta]))**2 for i in range(pc.numneu)] ) + #prob = 0.0; + #for i in range(pc.numneu) : + # prob += (np.absolute(MixMatrix[i][alpha]))**2*(np.absolute(MixMatrix[i][beta]))**2 + #return prob + +def OscProb(alpha, Enu, LVATERM = np.zeros((3,3), dtype=numpy.complex), + LVCTERM = np.zeros((3,3), dtype=numpy.complex)): + eigvals, eigvec = MT.eigenvectors(Hamiltonian(Enu, LVATERM=LVATERM, LVCTERM=LVCTERM)) + #print eigvec.dtype + if use_cython: + return [ clv.OscProbFromMixingMatrix(alpha,beta,eigvec) for beta in range(pc.numneu)] + else: + return [ OscProbFromMixingMatrix(alpha,beta,eigvec) for beta in range(pc.numneu)] + +def FlavorRatio(initial_flavor_ratio, Enu, LVATERM = np.zeros((3,3), dtype=numpy.complex), + LVCTERM = np.zeros((3,3), dtype=numpy.complex)): + final_flavor_ratio = [0.0]*pc.numneu + osc_prob_array = [OscProb(beta,Enu,LVATERM=LVATERM,LVCTERM=LVCTERM) for beta in range(pc.numneu)] + + for alpha in range(pc.numneu): + for beta,phi in enumerate(initial_flavor_ratio): + final_flavor_ratio[alpha] += osc_prob_array[beta][alpha]*phi + return final_flavor_ratio + +def RRR(initial_flavor_ratio, Enu, LVATERM = np.zeros((3,3), dtype=numpy.complex), + LVCTERM = np.zeros((3,3), dtype=numpy.complex)): + ffr = FlavorRatio(initial_flavor_ratio,Enu,LVATERM=LVATERM,LVCTERM=LVCTERM) + return ffr[1]/ffr[0] +def SSS(initial_flavor_ratio, Enu, LVATERM = np.zeros((3,3), dtype=numpy.complex), + LVCTERM = np.zeros((3,3), dtype=numpy.complex)): + ffr = FlavorRatio(initial_flavor_ratio,Enu,LVATERM=LVATERM,LVCTERM=LVCTERM) + return ffr[2]/ffr[1] + +def PointToList(p1,p2): + return [[p1[0],p2[0]],[p1[1],p2[1]]] + +def PointFromFlavor(origin,scale,flavor_ratio_list): + nu_e_vec = np.array([1.,0.])*scale + nu_mu_vec = np.array([1./2.,np.sqrt(3.)/2.])*scale + nu_tau_vec = np.array([-1./2.,np.sqrt(3.)/2.])*scale + fpos = origin + flavor_ratio_list[0]*nu_e_vec + flavor_ratio_list[1]*nu_mu_vec + return [fpos[0],fpos[1]] + +def MakeFlavorTriangle(list_of_flavor_ratios, scale = 8, + p = np.array([0.,0.]), save_file = False, PlotPoints = False, PlotTrayectories = False, figure = None, alpha = 1.0, + filename = "triangle",icolor = "green", icolormap = "Greens", divisions = 5, initial_flavor_ratio = [1,0,0], + term = "a", subdivisions = False, triangle_collection = None, color_scale = "lin", return_fig = True, addtext = "", + add_default_text = True, ilw = 1., slw = 0.75, output_format = "eps", inner_line_color = "k", plot_color_bar = False): + # i will be nice ... + list_of_flavor_ratios = np.array(list_of_flavor_ratios) + + if figure == None: + fig = plt.figure(figsize=(scale,scale), frameon = False) + else: + fig = figure + + ax = fig.add_axes([0, 0, 1, 1]) + ax.axis('off') + + # delete extra lines + frame = plt.gca() + frame.axes.get_xaxis().set_visible(False) + frame.axes.get_yaxis().set_visible(False) + + s0 = np.array([1.,0.])*scale + s1 = np.array([1./2.,np.sqrt(3.)/2.])*scale + s2 = np.array([1./2.,-np.sqrt(3.)/2.])*scale + + # make triangle outer frame + + plt.plot(*PointToList(p, p+s0), color = "k", lw = 4) + plt.plot(*PointToList(p, p+s1), color = "k", lw = 2) + plt.plot(*PointToList(p+s0, p+s1), color = "k", lw = 2) + + # put outer triangle labels + + # ax.text((p+s0*0.5+s0*0.025)[0], (p+s0*0.5-np.array([0,0.15*scale]))[1],r"$\alpha^{\oplus}_{e}$", + # horizontalalignment="center",fontsize = 50, zorder = 10) + # ax.text((p+s1*0.5-0.2*s0)[0], (p+s1*0.5+0.1*s0)[1]+scale*0.1,r"$\alpha^{\oplus}_{\tau}$", + # horizontalalignment="center",fontsize = 50, zorder = 10, rotation = 60.) + # ax.text((p+s1*0.5 + 0.7*s0)[0], (p+s1*0.5 + 0.6*s0)[1]+0.05*scale,r"$\alpha^{\oplus}_{\mu}$", + # horizontalalignment="center",fontsize = 50, zorder = 10, rotation = -60) + + ax.text((p+s0*0.5+s0*0.025)[0], (p+s0*0.5-np.array([0,0.15*scale]))[1],r"$f_{e,\oplus}$", + horizontalalignment="center",fontsize = 50, zorder = 10) + ax.text((p+s1*0.5-0.2*s0)[0], (p+s1*0.5+0.1*s0)[1]+scale*0.1,r"$f_{\tau, \oplus}$", + horizontalalignment="center",fontsize = 50, zorder = 10, rotation = 60.) + ax.text((p+s1*0.5 + 0.7*s0)[0], (p+s1*0.5 + 0.6*s0)[1]+0.05*scale,r"$f_{\mu, \oplus}$", + horizontalalignment="center",fontsize = 50, zorder = 10, rotation = -60) + + # construct triangle grid + fsl = 35 + for i in range(divisions+1): + subsize = 1./float(divisions) + + ax.text((p+s0*subsize*float(i))[0], (p+s0*subsize*float(i)-np.array([0,0.05*scale]))[1],str(i*subsize), + horizontalalignment="center",fontsize = fsl) + plt.plot(*PointToList(p+s0*subsize*float(i), p+s1+s2*subsize*float(i)), color = inner_line_color, lw = ilw, ls = "dashed", zorder = -1) + ax.text((p+s1-s1*subsize*float(i)-np.array([0.06*scale,0.0]))[0], (p+s1-s1*subsize*float(i))[1],str(i*subsize), + horizontalalignment="center",fontsize = fsl) + plt.plot(*PointToList(p+s0*subsize*float(divisions-i), p+s1-s1*subsize*float(i)), color = inner_line_color, lw = ilw, ls = "dashed", zorder = -1) + + ax.text((p+s1+s2*subsize*float(i)+np.array([0.05*scale,0.0]))[0], (p+s1+s2*subsize*float(i))[1],str((divisions-i)*subsize), + horizontalalignment="center",fontsize = fsl) + plt.plot(*PointToList(p+s1*subsize*float(divisions-i), p+s1+s2*subsize*float(i)), color = inner_line_color, lw = ilw, ls = "dashed", zorder = -1) + + if subdivisions and i < divisions: + plt.plot(*PointToList(p+s0*subsize*float(i+0.5), p+s1+s2*subsize*float(i+0.5)), color = inner_line_color, lw = slw, ls = "dotted", zorder = -1) + if subdivisions and i > 0: + plt.plot(*PointToList(p+s0*subsize*float(divisions-(i-0.5)), p+s1-s1*subsize*float(i-0.5)), color = inner_line_color, lw = slw, ls = "dotted", zorder = -1) + plt.plot(*PointToList(p+s1*subsize*float(divisions-(i-0.5)), p+s1+s2*subsize*float(i-0.5)), color = inner_line_color, lw = slw, ls = "dotted", zorder = -1) + + + # plot triangle collection + if (triangle_collection != None): + # get total number of points + total_points = float(sum([ triangle.number_of_points for triangle in triangle_collection])) + max_points = float(max([ triangle.number_of_points for triangle in triangle_collection])) + color_map = plt.get_cmap(icolormap) + for triangle in triangle_collection: + if triangle.number_of_points > 0: + xx,yy = zip(*triangle.coordinates) + if color_scale == "lin": + plt.fill(xx,yy,lw = 0., zorder = -0.8, color = color_map(0.75), alpha = np.sqrt(float(triangle.number_of_points)/max_points)) + #plt.fill(xx,yy,lw = 0., zorder = -0.8, color = color_map(float(triangle.number_of_points)/max_points), alpha = alpha) + elif color_scale == "log": + plt.fill(xx,yy,lw = 0., zorder = -0.8, color = color_map(0.75), alpha = (np.log10(float(triangle.number_of_points))/np.log10(max_points))) + #plt.fill(xx,yy,lw = 0., zorder = -0.8, color = color_map(0.7), alpha = (np.log10(float(triangle.number_of_points))/np.log10(max_points))**(2./3.)) + #plt.fill(xx,yy,lw = 0., zorder = -0.8, color = color_map(np.log10(float(triangle.number_of_points))/np.log10(max_points)), alpha = alpha) + #plt.fill(xx,yy,lw = 0., zorder = -0.8, color = color_map(np.log10(float(triangle.number_of_points)))) + else : + raise NameError('Error. Love CA.') + + if(plot_color_bar): + # the color bar magic + # location set on 0 to 1 scales. + left = 0.1 + bottom = -0.25 + width = 0.8 + height = 0.025 + cbaxes = fig.add_axes([left,bottom,width,height]) + if color_scale == "lin": + norm = mpl.colors.Normalize(vmin = 0., vmax = max_points) + elif color_scale == "log": + norm = mpl.colors.Normalize(vmin = 0., vmax = 1.0) + else : + raise NameError('Error. Love CA.') + mpl.rcParams.update({'font.size': 10}) + triangle_colorbar = mpl.colorbar.ColorbarBase(cbaxes, cmap = color_map, norm = norm, + orientation = "horizontal", spacing = "proportional", + # ) + format ='%.0e') + cbaxes.set_xlabel("Likelihood", fontsize = 12) + + # plot flavor ratio points + if PlotTrayectories : + if len(list_of_flavor_ratios.shape) == 3 : + for flavor_ratio_l in list_of_flavor_ratios: + flv_ratio_coords = map(lambda f:PointFromFlavor(p,scale,np.array(f)),flavor_ratio_l) + xc, yc = zip(*flv_ratio_coords) + plt.plot(xc,yc, color = icolor, + ms = 10, linewidth = 4, zorder = 0) + elif len(list_of_flavor_ratios.shape) == 2 : + flv_ratio_coords = map(lambda f:PointFromFlavor(p,scale,np.array(f)),list_of_flavor_ratios) + xc, yc = zip(*flv_ratio_coords) + + plt.plot(xc,yc, color = icolor, + ms = 10, linewidth = 4, zorder = 0) + else: + raise NameError('Check your input flavor list array and the joined flag. Love CA.') + elif PlotPoints: + if len(list_of_flavor_ratios.shape) !=2 : + print len(list_of_flavor_ratios.shape) + raise NameError('Check your input flavor list array and the joined flag. Love CA.') + for i,flavor_ratio in enumerate(list_of_flavor_ratios): + if len(icolor) != len(list_of_flavor_ratios): + icolor_ = icolor + else: + icolor_ = icolor[i] + plt.plot(*PointFromFlavor(p,scale,np.array(flavor_ratio)), color = icolor_, + marker = 'o', ms = 10, linewidth = 0, + markeredgecolor=None,markeredgewidth=0.1, zorder = 1000) + + # put back color scale axis + if add_default_text: + ax.text((s0/5.+0.9*s1)[0],(s0/5.+0.9*s1)[1], + "LV "+term+"-term scan with\n $\ \phi_e:\ \phi_\\mu:\ \phi_\\tau = "+str(initial_flavor_ratio[0])+":\ "+str(initial_flavor_ratio[1])+":\ "+str(initial_flavor_ratio[2])+"$"+" \n "+addtext, + fontsize = 20) + + if(save_file): + # save figure + plt.savefig("./plots/"+filename+"."+output_format, dpi = 600, bbox_inches='tight') + else: + # show figure + plt.show() + if return_fig: + return fig + + +def s_bario(p,p0,p1,p2): + return (p0[1]*p2[0] - p0[0]*p2[1] + (p2[1] - p0[1])*p[0] + (p0[0] - p2[0])*p[1]) + +def t_bario(p,p0,p1,p2): + return (p0[0]*p1[1] - p0[1]*p1[0] + (p0[1] - p1[1])*p[0] + (p1[0] - p0[0])*p[1]) + +def IsInTriangle(p,p0,p1,p2,area): + s = s_bario(p,p0,p1,p2) + t = t_bario(p,p0,p1,p2) + #print s,t,2.0*area - s - t + return s >= -1.e-15 and t >= -1.0e-15 and s+t <= 2.0*area + + +class Triangle: + coordinates = [] + area = 0.0 + number_of_points = 0.0 + n_t = 0 + i = 0 + j = 0 + orientation = "" + + def IsPointIn(self,point): + p0 = self.coordinates[0] + p1 = self.coordinates[1] + p2 = self.coordinates[2] + return IsInTriangle(point,p0,p1,p2,self.area) + + +def GenerateTriangles(scale, divisions, p = np.array([0.,0.])): + s0 = np.array([1.,0.])*scale/float(divisions) + s1 = np.array([1./2.,np.sqrt(3.)/2.])*scale/float(divisions) + s2 = np.array([1./2.,-np.sqrt(3.)/2.])*scale/float(divisions) + + area = np.sqrt(3)*(LA.norm(s0)/2.0)**2 + + n_t = 0 + + triangle_collection = [] + for i in range(divisions): + for j in range(divisions-i): + lower_triangle = Triangle() + + p0_l = p + i*s0 + j*s1 + p1_l = p0_l + s0 + p2_l = p0_l + s1 + + lower_triangle.coordinates = [p0_l,p1_l,p2_l] + lower_triangle.n_t = n_t + lower_triangle.i = i + lower_triangle.j = j + lower_triangle.orientation = "L" + lower_triangle.area = area + + n_t += 1 + # append to triangle collection + triangle_collection.append(lower_triangle) + + upper_triangle = Triangle() + + p0_u = p2_l + p1_u = p1_l + p2_u = p1_l + s1 + + upper_triangle.coordinates = [p0_u,p1_u,p2_u] + upper_triangle.n_t = n_t + upper_triangle.i = i + upper_triangle.j = j + upper_triangle.orientation = "U" + upper_triangle.area = area + + n_t += 1 + # append to triangle collection + triangle_collection.append(upper_triangle) + return triangle_collection + +def AddPointToTriangleCollectionLegacy(flavor_ratio, triangle_collection, + p = np.array([0.,0.]), scale = 8, divisions = 10): + point = PointFromFlavor(p,scale,np.array(flavor_ratio)) + electron = 0; tau = 2; + # the silly way + for triangle in triangle_collection: + if(triangle.IsPointIn(point)): + triangle.number_of_points += 1. + +def AddPointToTriangleCollection(flavor_ratio, triangle_collection, + p = np.array([0.,0.]), scale = 8, divisions = 10): + point = PointFromFlavor(p,scale,np.array(flavor_ratio)) + electron = 0; muon = 1; tau = 2; + # the smart way + u_i = int(flavor_ratio[electron]*float(divisions)) + u_j = int(flavor_ratio[muon]*float(divisions)) + index = u_i*(2*divisions-u_i+1) + 2*u_j + if triangle_collection[index].IsPointIn(point): + triangle_collection[index].number_of_points += 1. + else: + triangle_collection[index+1].number_of_points += 1. +# legacy + #elif triangle_collection[index+1].IsPointIn(point): + # triangle_collection[index+1].number_of_points += 1. + #else: + # print "Math error." + # print point, "\n",u_i, u_j, "\n", triangle_collection[index].coordinates, "\n", triangle_collection[index+1].coordinates + # raise NameError("Error triangle location math") + +class AnarchySampling: + def __init__(self, n_sample, LV_scale_1, LV_scale_2, term): + self.n_sample = n_sample + self.th12_sample = np.arcsin(np.sqrt(np.random.uniform(0.,1., size=n_sample))) + self.th13_sample = np.arccos(np.sqrt(np.sqrt(np.random.uniform(0.,1., size=n_sample)))) + self.th23_sample = np.arcsin(np.sqrt(np.random.uniform(0.,1., size=n_sample))) + self.delta_sample = np.random.uniform(0.,2.*np.pi, size=n_sample) + + self.LV_scale_1 = LV_scale_1 + self.LV_scale_2 = LV_scale_2 + + self.term = term + +def GenerateFlavorRatioPoints(Initial_Flavor_Ratio, SamplingObject, gamma = 2.0, + Log10Emax = 7., Log10Emin = 4.0, Epoints = 30, + save_list = False, save_avg = True): + flavor_tray_list = [] + flavor_avg_list = [] + + # energy things + + Erange = np.logspace(Log10Emin,Log10Emax,Epoints) # in GeV + Emin = Erange[0] + Emax = Erange[-1] + + if gamma == 1 or gamma == 1.0: + spectral_normalization = np.log(Emax)-np.log(Emin) + else: + spectral_normalization = (Emax**(1.-gamma) - Emin**(1.-gamma))/(1.-gamma) + + spectral_function = lambda Enu: Enu**(-gamma)/spectral_normalization + + # loop over random parameters + for i in range(SamplingObject.n_sample): + lv_term = MT.LVP() + + lv_term.th12 = SamplingObject.th12_sample[i] + lv_term.th13 = SamplingObject.th13_sample[i] + lv_term.th23 = SamplingObject.th23_sample[i] + lv_term.delta1 = SamplingObject.delta_sample[i] + + lv_term.LVS21 = SamplingObject.LV_scale_1 + lv_term.LVS31 = SamplingObject.LV_scale_2 + + lv_term.Refresh() + + LVTERM = MT.LVTerm(lv_term); + + if SamplingObject.term == "a": + flavor_ratio_list = np.array(map(lambda Enu : FlavorRatio(Initial_Flavor_Ratio, Enu*pc.GeV, LVATERM = LVTERM), Erange)) + elif SamplingObject.term == "c": + flavor_ratio_list = np.array(map(lambda Enu : FlavorRatio(Initial_Flavor_Ratio, Enu*pc.GeV, LVCTERM = LVTERM), Erange)) + else : + raise NameError('Only a or c term.'+ str(term)) + + if save_avg: + if Epoints != 1: + flavor_avg = [0.]*lv_term.numneu + for alpha in range(lv_term.numneu): + #inter = interpolate.interp1d(Erange,flavor_ratio_list[:,alpha]) + inter = interpolate.UnivariateSpline(Erange,flavor_ratio_list[:,alpha]) + flavor_avg[alpha] = integrate.quad(lambda Enu : inter(Enu)*spectral_function(Enu), + Emin,Emax, limit = 75, epsrel = 1e-2, epsabs = 1.0e-2)[0] + #flavor_avg[alpha] = integrate.quadrature(lambda Enu : inter(Enu)*spectral_function(Enu), + # Emin,Emax, maxiter = 75, rtol = 1e-3, tol = 1.e-3)[0] + flavor_avg_list.append(flavor_avg) + else: + flavor_avg = flavor_ratio_list[0] + flavor_avg_list.append(flavor_avg) + + if save_list: + flavor_tray_list.append(flavor_ratio_list) + + if save_list and save_avg: + return flavor_tray_list, flavor_avg_list + elif save_list: + return flavor_tray_list + elif save_avg: + return flavor_avg_list + else : + print "Math is broken." + return None diff --git a/plot_llh/MinimalTools.py b/plot_llh/MinimalTools.py new file mode 100644 index 0000000..4ae8360 --- /dev/null +++ b/plot_llh/MinimalTools.py @@ -0,0 +1,166 @@ +import numpy as np +from PhysConst import PhysicsConstants + +def eigenvectors(M): + """ Calculates the eigenvectors and eigenvalues ordered by eigenvalue size + + @type M : matrix + @param M : matrix M + + @rtype : list + @return : [eigenvalues list, eigenvector list] + """ + D,V = np.linalg.eig(M) + DV = [] + VT = V.T + for i,eigenvalue in enumerate(D): + DV.append([eigenvalue,VT[i]]) + + DV = sorted(DV,key = lambda x : x[0].real)#np.abs(x[0].real)) + + V2 = [] + D2 = [] + for e in DV: + V2.append(e[1]) + D2.append(e[0]) + return D2,V2 + +# General Rotation Matrix +def R(i,j,cp,param): + """ Rotation Matrix + Calculates the R_ij rotations. Also incorporates CP-phases when necesary. + @type i : int + @param i : i-column. + @type j : int + @param j : j-row. + @type cp : int + @param cp : if cp = 0 : no CP-phase. else CP-phase = CP_array[cp] + + @rtype : numpy.array + @return : returns the R_ij rotation matrix. + """ + # if cp = 0 -> no complex phase + # R_ij, i<j + if(j<i): + # no funny business + l = i + i = j + j = l + + # rotation matrix generator + R = np.zeros([param.numneu,param.numneu],complex) + # diagonal terms + for k in np.arange(0,param.numneu,1): + if(k != i-1 and k != j-1): + R[k,k] = 1.0 + else : + R[k,k] = param.c[i,j] + # non-diagonal terms + if(cp != 0): + sd = np.sin(param.dcp[cp]) + cd = np.cos(param.dcp[cp]) + faseCP = complex(cd,sd) + else: + faseCP = complex(1.0,0.0) + + R[i-1,j-1] = param.s[i,j]*faseCP.conjugate() + R[j-1,i-1] = -param.s[i,j]*faseCP + return R + +def calcU(param): + """ Defining the mixing matrix parametrization. + @type param : PhysicsConstants + @param param : set of physical parameters to be used. + + @rtype : None + @return : Sets mixing matrix. + """ + if(param.numneu == 2): + return self.R(1,2,0,param) + elif(param.numneu == 3): + return np.dot(R(2,3,0,param),np.dot(R(1,3,1,param),R(1,2,0,param))) + elif(param.numneu == 4): + return np.dot(R(3,4,0,param),np.dot(R(2,4,2,param),np.dot(R(1,4,0,param),np.dot(R(2,3,0,param),np.dot(R(1,3,1,param),R(1,2,0,param)))))) + elif(param.numneu == 5): + return np.dot(R(4,5,0,param),np.dot(R(3,5,0,param),np.dot(R(2,5,0,param),np.dot(R(1,5,3,param),np.dot(R(3,4,0,param),np.dot(R(2,4,0,param),np.dot(R(1,4,2,param),np.dot(R(2,3,0,param),np.dot(R(1,3,1,param),R(1,2,0,param)))))))))) + elif(param.numneu == 6): + # 3+3 twin-sterile-neutrino model + return np.dot(R(3,6,0,param),np.dot(R(2,5,0,param),np.dot(R(1,4,0,param),np.dot(R(2,3,0,param),np.dot(R(1,3,1,param),R(1,2,0,param)))))) + else: + print "Sorry, too many neutrinos. Not yet implemented! =(." + quit() + + # antineutrino case + if param.neutype == "antineutrino" : + return self.U.conjugate() + + + +def massM2(param): + """ Mass term in the neutrino mass basis. + + @type param : PhysicsConstants + @param param : set of physical parameters to be used. + + @rtype : numpy array + @return : mass matrix in mass basis. + """ + M2 = np.zeros([param.numneu,param.numneu],complex) + for k in np.arange(1,param.numneu,1): + M2[k,k] = param.dmsq[k+1] + return M2 + +def flavorM2(param): + """ Mass term in the neutrino flavor basis. + + @type param : PhysicsConstants + @param param : set of physical parameters to be used. + + @rtype : numpy array + @return : mass matrix in flavor basis. + """ + U = calcU(param) + return np.dot(U,np.dot(massM2(param),U.conjugate().T)) + +class LVP(PhysicsConstants): + def __init__(self): + super(LVP,self).__init__() + self.th12 = 0.0 + self.th13 = 0.0 + self.th23 = 0.0 + self.delta1 = 0.0 + self.deltaCP = 0.0 + super(LVP,self).Refresh() + + # LVS + self.LVS21 = 0.0 # + self.LVS31 = 0.0 # + self.LVS41 = 0.0 # + self.LVS51 = 0.0 # + self.LVS61 = 0.0 # + # SQUARED MASS DIFFERENCE MATRIX + self.LVS = np.zeros([self.numneumax+2],float) + self.LVS[2] = self.LVS21 + self.LVS[3] = self.LVS31 + self.LVS[4] = self.LVS41 + self.LVS[5] = self.LVS51 + self.LVS[6] = self.LVS61 + + def Refresh(self): + super(LVP,self).Refresh() + LVS = self.LVS + LVS[2] = self.LVS21 + LVS[3] = self.LVS31 + LVS[4] = self.LVS41 + LVS[5] = self.LVS51 + LVS[6] = self.LVS61 + +def DiagonalMatrixLV(param): + DD = np.zeros([param.numneu,param.numneu],complex) + for k in np.arange(1,param.numneu,1): + DD[k,k] = param.LVS[k+1] + return DD + +def LVTerm(LVparam): + U = calcU(LVparam) + return np.dot(U,np.dot(DiagonalMatrixLV(LVparam),U.conjugate().T)) diff --git a/plot_llh/PhysConst.py b/plot_llh/PhysConst.py new file mode 100644 index 0000000..89a0be0 --- /dev/null +++ b/plot_llh/PhysConst.py @@ -0,0 +1,390 @@ +""" +Author : C.A. Arguelles +Date : 10/MAY/2011 + +Contains Physics constants and global variables. + +Log : +- Modified on 23/ABR/2012 by C.Arguelles + + Changed the definition of PhysicsConstants to + include an __init__ to separate the class global + properties from its instances. +""" + +# python standard modules +import numpy as np + +class PhysicsConstants(object): + + def __init__(self): + ## PHYSICS CONSTANTS + #=========================================================================== + # NAME + #=========================================================================== + + self.name = "STD" # Default values + self.linestyle = "solid" # Default linestyle in plots + self.markerstyle = "*" # Default marker style + self.colorstyle = "red" # Default color style + self.savefilename = "output.dat" # Default color style + + #=============================================================================== + # ## MATH + #=============================================================================== + self.PI=3.14159265 # Pi + self.PIby2=1.5707963268 # Pi/2 + self.sqr2=1.4142135624 # Sqrt[2] + self.ln2 = np.log(2.0) + + #=============================================================================== + # ## EARTH + #=============================================================================== + self.EARTHRADIUS = 6371.0 # [km] Earth radius + #=============================================================================== + # ## SUN + #=============================================================================== + self.SUNRADIUS = 109*self.EARTHRADIUS # [km] Sun radius + + #=============================================================================== + # # PHYSICAL CONSTANTS + #=============================================================================== + self.GF = 1.16639e-23 # [eV^-2] Fermi Constant + self.Na = 6.0221415e+23 # [mol cm^-3] Avogadro Number + self.sw_sq = 0.2312 # [dimensionless] sin(th_weinberg) ^2 + self.G = 6.67300e-11 # [m^3 kg^-1 s^-2] + self.alpha = 1.0/137.0 # [dimensionless] fine-structure constant + + #=============================================================================== + # ## UNIT CONVERSION FACTORS + #=============================================================================== + # Energy + self.TeV = 1.0e12 # [eV/TeV] + self.GeV = 1.0e9 # [eV/GeV] + self.MeV = 1.0e6 # [eV/MeV] + self.keV = 1.0e3 # [eV/keV] + self.Joule = 1/1.60225e-19 # [eV/J] + # Mass + self.kg = 5.62e35 # [eV/kg] + self.gr = 1e-3*self.kg # [eV/g] + # Time + self.sec = 1.523e15 # [eV^-1/s] + self.hour = 3600.0*self.sec # [eV^-1/h] + self.day = 24.0*self.hour # [eV^-1/d] + self.year = 365.0*self.day # [eV^-1/yr] + self.yearstosec = self.sec/self.year # [s/yr] + # Distance + self.meter = 5.076e6 # [eV^-1/m] + self.cm = 1.0e-2*self.meter # [eV^-1/cm] + self.km = 1.0e3*self.meter # [eV^-1/km] + self.fermi = 1.0e-15*self.meter # [eV^-1/fm] + self.angstrom = 1.0e-10*self.meter # [eV^-1/A] + self.AU = 149.60e9*self.meter # [eV^-1/AU] + self.parsec = 3.08568025e16*self.meter# [eV^-1/parsec] + # Integrated Luminocity # review + self.picobarn = 1.0e-36*self.cm**2 # [eV^-2/pb] + self.femtobarn = 1.0e-39*self.cm**2 # [eV^-2/fb] + # Presure + self.Pascal = self.Joule/self.meter**3 # [eV^4/Pa] + self.hPascal = 100.0*self.Pascal # [eV^4/hPa] + self.atm = 101325.0*self.Pascal # [eV^4/atm] + self.psi = 6893.0*self.Pascal # [eV^4/psi] + # Temperature + self.kelvin = 1/1.1604505e4 # [eV/K] + # Angle + self.degree = self.PI/180.0 # [rad/degree] + # magnetic field + self.T = 0.000692445 # [eV^2/T] + + # old notation + self.cm3toev3 = 7.68351405e-15 # cm^3-> ev^3 + self.KmtoEv =5.0677288532e+9 # km -> eV + self.yearstosec = 31536.0e3 # years -> sec + + #=============================================================================== + # ## NEUTRINO OSCILLATION PARAMETERS ## + #=============================================================================== + + self.numneu = 3 # number of neutrinos + self.numneumax = 6 # maximum neutrino number + self.neutype = 'neutrino' + #neutype = 'antineutrino' + + # values updated according to 1209.3023 Table 1 FreeFluxes + RSBL + + # MIXING ANGLES + + self.th12 = 0.579639 + self.th13 = 0.150098 + self.th23 = self.PIby2/2.0 + self.th14 = 0.0 + self.th24 = 0.0 + self.th34 = 0.0 + self.th15 = 0.0 + self.th25 = 0.0 + self.th35 = 0.0 + self.th45 = 0.0 + self.th16 = 0.0 + self.th26 = 0.0 + self.th36 = 0.0 + self.th46 = 0.0 + self.th56 = 0.0 + + # mixing angles matrix array + self.th = np.zeros([self.numneumax+1,self.numneumax+1],float) + self.th[1,2] = self.th12 + self.th[1,3] = self.th13 + self.th[2,3] = self.th23 + self.th[1,4] = self.th14 + self.th[2,4] = self.th24 + self.th[3,4] = self.th34 + self.th[1,5] = self.th15 + self.th[2,5] = self.th25 + self.th[3,5] = self.th35 + self.th[4,5] = self.th45 + self.th[1,6] = self.th16 + self.th[2,6] = self.th26 + self.th[3,6] = self.th36 + self.th[4,6] = self.th46 + self.th[5,6] = self.th56 + + self.s12 = np.sin(self.th12) + self.c12 = np.cos(self.th12) + self.s13 = np.sin(self.th13) + self.c13 = np.cos(self.th13) + self.s23 = np.sin(self.th23) + self.c23 = np.cos(self.th23) + self.s14 = np.sin(self.th14) + self.c14 = np.cos(self.th14) + self.s24 = np.sin(self.th24) + self.c24 = np.cos(self.th24) + self.s34 = np.sin(self.th34) + self.c34 = np.cos(self.th34) + self.s15 = np.sin(self.th15) + self.c15 = np.cos(self.th15) + self.s25 = np.sin(self.th25) + self.c25 = np.cos(self.th25) + self.s35 = np.sin(self.th35) + self.c35 = np.cos(self.th35) + self.s45 = np.sin(self.th45) + self.c45 = np.cos(self.th45) + self.s16 = np.sin(self.th16) + self.c16 = np.cos(self.th16) + self.s26 = np.sin(self.th26) + self.c26 = np.cos(self.th26) + self.s36 = np.sin(self.th36) + self.c36 = np.cos(self.th36) + self.s46 = np.sin(self.th46) + self.c46 = np.cos(self.th46) + self.s56 = np.sin(self.th56) + self.c56 = np.cos(self.th56) + + # cos(th_ij) matrix array + self.c = np.zeros([self.numneumax+1,self.numneumax+1],float) + self.c[1,2] = self.c12 + self.c[1,3] = self.c13 + self.c[1,4] = self.c14 + self.c[2,3] = self.c23 + self.c[2,4] = self.c24 + self.c[3,4] = self.c34 + self.c[1,5] = self.c15 + self.c[2,5] = self.c25 + self.c[3,5] = self.c35 + self.c[4,5] = self.c45 + self.c[1,6] = self.c16 + self.c[2,6] = self.c26 + self.c[3,6] = self.c36 + self.c[4,6] = self.c46 + self.c[5,6] = self.c56 + + # sin(th_ij) matrix array + self.s = np.zeros([self.numneumax+1,self.numneumax+1],float) + self.s[1,2] = self.s12 + self.s[1,3] = self.s13 + self.s[1,4] = self.s14 + self.s[2,3] = self.s23 + self.s[2,4] = self.s24 + self.s[3,4] = self.s34 + self.s[1,5] = self.s15 + self.s[2,5] = self.s25 + self.s[3,5] = self.s35 + self.s[4,5] = self.s45 + self.s[1,6] = self.s16 + self.s[2,6] = self.s26 + self.s[3,6] = self.s36 + self.s[4,6] = self.s46 + self.s[5,6] = self.s56 + + # CP PHASES + #self.delta21=3.3e-5 + #self.delta32=3.1e-3 + #self.delta31=self.delta32+self.delta21 + #self.deltaCP=self.PIby2 + + # CP Phases + self.deltaCP = 5.235987 + self.delta1 = self.deltaCP + self.delta2 = 0.0 + self.delta3 = 0.0 + + # d-CP phases + self.dcp = np.zeros([self.numneumax-2+1],complex) + self.dcp[0] = 1.0 + self.dcp[1] = self.delta1 + self.dcp[2] = self.delta2 + self.dcp[3] = self.delta3 + + # SQUARED MASS DIFFERENCE + self.dm21sq = 7.50e-5 # [eV^2] + self.dm31sq = 2.47e-3 # [eV^2] + self.dm32sq = -2.43e-3 # [eV^2] + # STERILE + self.dm41sq = 0.0 # [eV^2] + self.dm51sq = 0.0 # [eV^2] + self.dm61sq = 0.0 # [eV^2] + # SQUARED MASS DIFFERENCE MATRIX + self.dmsq = np.zeros([self.numneumax+2],float) + self.dmsq[2] = self.dm21sq + self.dmsq[3] = self.dm31sq + self.dmsq[4] = self.dm41sq + self.dmsq[5] = self.dm51sq + self.dmsq[6] = self.dm61sq + + self.dm2 = np.zeros([self.numneumax+1,self.numneumax+1],float) + self.dm2[1,2] = self.dm21sq + self.dm2[1,3] = self.dm31sq + self.dm2[2,3] = self.dm32sq + self.dm2[1,4] = self.dm41sq + self.dm2[1,5] = self.dm51sq + self.dm2[1,6] = self.dm61sq + + # MIXING MATRIX + self.U = None + + #=============================================================================== + # # PARTICLE MASSES + #=============================================================================== + self.muon_mass = 0.10565 # [GeV] + self.neutron_mass = 0.939565 # [GeV] + self.proton_mass = 0.938272 # [GeV] + self.electron_mass = 0.510998910e-3 # [GeV] + + self.atomic_mass_unit = 1.660538e-24 # [g] + + ## names + self.electron = 0 + self.muon = 1 + self.tau = 2 + self.sterile1 = 3 + self.sterile2 = 4 + self.sterile3 = 5 + + #=============================================================================== + # REFRESH + #=============================================================================== + + def Refresh(self): + # Refresh angles + self.s12 = np.sin(self.th12) + self.c12 = np.cos(self.th12) + self.s13 = np.sin(self.th13) + self.c13 = np.cos(self.th13) + self.s23 = np.sin(self.th23) + self.c23 = np.cos(self.th23) + self.s14 = np.sin(self.th14) + self.c14 = np.cos(self.th14) + self.s24 = np.sin(self.th24) + self.c24 = np.cos(self.th24) + self.s34 = np.sin(self.th34) + self.c34 = np.cos(self.th34) + self.s15 = np.sin(self.th15) + self.c15 = np.cos(self.th15) + self.s25 = np.sin(self.th25) + self.c25 = np.cos(self.th25) + self.s35 = np.sin(self.th35) + self.c35 = np.cos(self.th35) + self.s45 = np.sin(self.th45) + self.c45 = np.cos(self.th45) + self.s16 = np.sin(self.th16) + self.c16 = np.cos(self.th16) + self.s26 = np.sin(self.th26) + self.c26 = np.cos(self.th26) + self.s36 = np.sin(self.th36) + self.c36 = np.cos(self.th36) + self.s46 = np.sin(self.th46) + self.c46 = np.cos(self.th46) + self.s56 = np.sin(self.th56) + self.c56 = np.cos(self.th56) + + th = self.th + th[1,2] = self.th12 + th[1,3] = self.th13 + th[2,3] = self.th23 + th[1,4] = self.th14 + th[2,4] = self.th24 + th[3,4] = self.th34 + th[1,5] = self.th15 + th[2,5] = self.th25 + th[3,5] = self.th35 + th[4,5] = self.th45 + th[1,6] = self.th16 + th[2,6] = self.th26 + th[3,6] = self.th36 + th[4,6] = self.th46 + th[5,6] = self.th56 + # Refresh cos(th_ij) + c = self.c + c[1,2] = self.c12 + c[1,3] = self.c13 + c[1,4] = self.c14 + c[2,3] = self.c23 + c[2,4] = self.c24 + c[3,4] = self.c34 + c[1,5] = self.c15 + c[2,5] = self.c25 + c[3,5] = self.c35 + c[4,5] = self.c45 + c[1,6] = self.c16 + c[2,6] = self.c26 + c[3,6] = self.c36 + c[4,6] = self.c46 + c[5,6] = self.c56 + # Refresh sin(th_ij) + s = self.s + self.s[1,2] = self.s12 + self.s[1,3] = self.s13 + self.s[1,4] = self.s14 + self.s[2,3] = self.s23 + self.s[2,4] = self.s24 + self.s[3,4] = self.s34 + self.s[1,5] = self.s15 + self.s[2,5] = self.s25 + self.s[3,5] = self.s35 + self.s[4,5] = self.s45 + self.s[1,6] = self.s16 + self.s[2,6] = self.s26 + self.s[3,6] = self.s36 + self.s[4,6] = self.s46 + self.s[5,6] = self.s56 + # Refresh CP-Phases + dcp = self.dcp + dcp[0] = 1.0 + dcp[1] = self.delta1 + dcp[2] = self.delta2 + dcp[3] = self.delta3 + #dcp[4] = self.delta2 + # Refresh Square mass differences + dmsq = self.dmsq + dmsq[2] = self.dm21sq + dmsq[3] = self.dm31sq + dmsq[4] = self.dm41sq + dmsq[5] = self.dm51sq + dmsq[6] = self.dm61sq + + dm2 = self.dm2 + dm2[1,2] = self.dm21sq + dm2[1,3] = self.dm31sq + dm2[2,3] = self.dm32sq + dm2[1,4] = self.dm41sq + dm2[1,5] = self.dm51sq + dm2[1,6] = self.dm61sq + diff --git a/plot_llh/make_plots.py b/plot_llh/make_plots.py new file mode 100644 index 0000000..9216396 --- /dev/null +++ b/plot_llh/make_plots.py @@ -0,0 +1,117 @@ +#!/usr/bin/python + +import sys +sys.path.append('/Users/Teps/Theory') +#import header as h +#sys.path.append('/Users/Teps/Theory/HESE') +#import anarchy_header as ah +sys.path.append('/Users/Teps/Theory/HESE/Carlos') +import matplotlib as mpl +import matplotlib.pyplot as plt +from matplotlib import rc, rcParams +import MinimalTools as MT +import PhysConst as PC +import LVCPT as LVCPT +import numpy as np + +import sys,os + +rc('text', usetex=True) +rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18}) +cols = ['#29A2C6','#FF6D31','#FFCB18','#73B66B','#EF597B', '#333333'] + +font = {'family' : 'serif', + 'weight' : 'bold', + 'size' : 18} + +if len(sys.argv)< 2: + print sys.argv + print "Usage: make_plots.py input_filepath." + exit(1) + +#colors_schemes = ["Greens","Reds","Blues","PuRd","summer"] +colors_schemes = ["Greens","Reds","Blues","spring","summer"] +#colors_schemes = ["Greens","Reds","Blues","cool","summer"] +#colors_schemes = ["Greens","Reds","Blues","PuRd","summer"] +#colors_schemes = ["Blues","Greens","Reds","PuRd","summer"] +#colors_schemes = ["Greys","Greens","Reds","PuRd","summer"] +#colors_schemes = ["Greys","Greys","Greys","Greys","summer"] +#colors_schemes = ["PuRd","summer"] +output_format = "pdf" + +# if True then will plot all files in the same figure +use_same_canvas = True +figure = None + +for i in range(len(sys.argv)-1): + infile = sys.argv[i+1] + print "Load data: " + str(infile) + if infile[-3:] == 'txt': + flavor_list = np.genfromtxt(infile) + else: + flavor_list = np.load(infile) + if len(flavor_list[~np.isfinite(flavor_list)]) != 0: + fl = [] + for x in flavor_list: + if np.sum(~np.isfinite(x)) == 0: + fl.append(x.tolist()) + flavor_list = np.array(fl) + print flavor_list + print "Done loading data" + + if not use_same_canvas : + filename = "triangle_plot_"+os.path.splitext(sys.argv[i+1])[0] + else : + filename = "triangle_plot" + + # plots scale and diviions + scale = 8 + divisions = 40 + + print "Begin making plot ..." + triangle_collection = LVCPT.GenerateTriangles(scale,divisions*2) + map(lambda f : LVCPT.AddPointToTriangleCollection(f,triangle_collection, scale = scale, divisions = divisions*2),flavor_list) + + if use_same_canvas: + figure = LVCPT.MakeFlavorTriangle(flavor_list, divisions = 5, save_file=True, + filename = filename + "_hist", icolor = "g", icolormap = colors_schemes[i], + triangle_collection=triangle_collection, figure = figure, alpha = 0.8, + initial_flavor_ratio = [0,1,0], subdivisions = True, color_scale = "lin", + output_format = output_format, inner_line_color ="silver",add_default_text = False, + plot_color_bar =True) + + else: + figure = LVCPT.MakeFlavorTriangle(flavor_list, divisions = 5, save_file=True, + filename = filename + "_hist", icolor = "g", icolormap = colors_schemes[i], + triangle_collection=triangle_collection, alpha = 0.8, + initial_flavor_ratio = [0,1,0], subdivisions = True, color_scale = "lin", + output_format = output_format, inner_line_color = "silver",add_default_text = False, + plot_color_bar =True) + + print "Done making plot: " + filename + "_hist."+output_format + +ax = figure.get_axes()[0] +#ax.plot(6.5-0.35,5.5+0.3+0.35,"o", color = "grey", ms = 10, markeredgewidth = 0.1, alpha = 0.9) +#ax.text(6.7-0.35,5.44+0.3+0.35,r"$(1-x:x:0)$", fontsize = 16) +#ax.axvline(x = 7.9) +fsz = 32 +ax.plot(6.5-0.1,5.5+0.3+0.35+0.2+0.2,"o", color = "gold", ms = 25, markeredgewidth = 0.1, alpha = 0.9) +ax.text(6.7-0.1,5.44+0.3+0.35+0.2+0.2,r"$(1:2:0)$", fontsize = fsz) + +ax.plot(6.5-0.1,5.5+0.35+0.2,"o", color = "#2B653E", ms = 25, markeredgewidth = 0.1, alpha = 0.9) +ax.text(6.7-0.1,5.44+0.35+0.2,r"$(1:0:0)$", fontsize = fsz) + +ax.plot(6.5-0.1,5.5-0.3+0.35-0.2+0.2,"o", color = "#CA323D", ms = 25, markeredgewidth = 0.1, alpha = 0.9) +ax.text(6.7-0.1,5.44-0.3+0.35-0.2+0.2,r"$(0:1:0)$", fontsize = fsz) + +ax.plot(6.5-0.1,5.5-0.3+0.35-0.3-0.4+0.2,"o", color = "#2D4676", ms = 25, markeredgewidth = 0.1, alpha = 0.9) +ax.text(6.7-0.1,5.44-0.3+0.35-0.3-0.4+0.2,r"$(0:0:1)$", fontsize = fsz) + +plt.savefig("./plots/"+filename+"."+output_format, dpi = 600, bbox_inches='tight') + +exit(1) + +##os.system("cd plots") +##os.system("gs triangle_plot_hist.eps") +##os.system("cd ..") + diff --git a/plot_llh/mcmc_scan.py b/plot_llh/mcmc_scan.py new file mode 100644 index 0000000..455d5dd --- /dev/null +++ b/plot_llh/mcmc_scan.py @@ -0,0 +1,174 @@ +#! /usr/bin/env python +""" +Sample points from a gaussian likelihood +""" + +from __future__ import absolute_import, division + +import sys + +import argparse +import multiprocessing + +import numpy as np +from scipy.stats import multivariate_normal + +import emcee +import tqdm + +MEASURED_FR = [1, 1, 1] +SIGMA = 0.001 + + +def angles_to_fr(angles): + sphi4, c2psi = angles + + psi = (0.5)*np.arccos(c2psi) + + sphi2 = np.sqrt(sphi4) + cphi2 = 1. - sphi2 + spsi2 = np.sin(psi)**2 + cspi2 = 1. - spsi2 + + x = sphi2*cspi2 + y = sphi2*spsi2 + z = cphi2 + return x, y, z + + +def triangle_llh(theta): + """-Log likelihood function for a given theta.""" + fr = angles_to_fr(theta) + fr_bf = MEASURED_FR + cov_fr = np.identity(3) * SIGMA + return np.log(multivariate_normal.pdf(fr, mean=fr_bf, cov=cov_fr)) + + +def lnprior(theta): + """Priors on theta.""" + sphi4, c2psi = theta + + # Flavour ratio bounds + if 0. <= sphi4 <= 1.0 and -1.0 <= c2psi <= 1.0: + pass + else: return -np.inf + + return 0. + + +def lnprob(theta): + """Prob function for mcmc.""" + lp = lnprior(theta) + if not np.isfinite(lp): + return -np.inf + return lp + triangle_llh(theta) + + +def parse_args(): + """Parse command line arguments""" + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + '--measured-ratio', type=int, nargs=3, default=[1, 1, 1], + help='Set the central value for the measured flavour ratio at IceCube' + ) + parser.add_argument( + '--sigma-ratio', type=float, default=0.01, + help='Set the 1 sigma for the measured flavour ratio' + ) + parser.add_argument( + '--burnin', type=int, default=100, + help='Amount to burnin' + ) + parser.add_argument( + '--nwalkers', type=int, default=100, + help='Number of walkers' + ) + parser.add_argument( + '--nsteps', type=int, default=2000, + help='Number of steps to run' + ) + parser.add_argument( + '--seed', type=int, default=99, + help='Set the random seed value' + ) + parser.add_argument( + '--outfile', type=str, default='./untitled', + help='Path to output chains' + ) + args = parser.parse_args() + return args + + +def main(): + args = parse_args() + + np.random.seed(args.seed) + + global MEASURED_FR + global SIGMA + + MEASURED_FR = np.array(args.measured_ratio) / float(np.sum(args.measured_ratio)) + SIGMA = args.sigma_ratio + + print 'MEASURED_FR = {0}'.format(MEASURED_FR) + print 'SIGMA = {0}'.format(SIGMA) + + ndim = 2 + nwalkers = args.nwalkers + ntemps = 1 + burnin = args.burnin + betas = np.array([1e0, 1e-1, 1e-2, 1e-3, 1e-4]) + p0_base = [0.5, 0.5] + p0_std = [.2, 0.2] + + print 'p0_base', p0_base + print 'p0_std', p0_std + p0 = np.random.normal(p0_base, p0_std, size=[ntemps, nwalkers, ndim]) + print map(lnprior, p0[0]) + + # threads = multiprocessing.cpu_count() + threads = 1 + sampler = emcee.PTSampler( + ntemps, nwalkers, ndim, triangle_llh, lnprior, threads=threads + ) + + print "Running burn-in" + for result in tqdm.tqdm(sampler.sample(p0, iterations=burnin), total=burnin): + pos, prob, state = result + sampler.reset() + print "Finished burn-in" + + nsteps = args.nsteps + + print "Running" + for _ in tqdm.tqdm(sampler.sample(pos, iterations=nsteps), total=nsteps): + pass + print "Finished" + + outfile = args.outfile+'_{0:03d}_{1:03d}_{2:03d}_{3:.1E}'.format( + int(MEASURED_FR[0]*100), int(MEASURED_FR[1]*100), int(MEASURED_FR[2]*100), SIGMA + ) + + samples = sampler.chain[0, :, :, :].reshape((-1, ndim)) + print 'acceptance fraction', sampler.acceptance_fraction + print 'sum of acceptance fraction', np.sum(sampler.acceptance_fraction) + print 'np.unique(samples[:,0]).shape', np.unique(samples[:,0]).shape + + try: + print 'autocorrelation', sampler.acor + except: + print 'WARNING : NEED TO RUN MORE SAMPLES FOR FILE ' + outfile + print 'outfile = ', outfile + + fr_samples = np.array(map(angles_to_fr, samples)) + np.save(outfile+'.npy', fr_samples) + + print "DONE!" + + +main.__doc__ = __doc__ + + +if __name__ == '__main__': + main() + |