# #CAMEL tools library # #Sep 2015 - M. Tristram - import numpy as np import matplotlib.pyplot as plt import scipy.ndimage as nd import sys import re import glob,os import gzip import mmap import pdb def readCosmoMC( filename, nchain=4): cosmopars = list(np.loadtxt( "%s.paramnames" % (filename), dtype='str', unpack=True, usecols=[0])) parlist = ['iter','chi2']+[from_cosmomcName.get(p,p) for p in cosmopars] data = extract_chains( filename+"_", nchain) chain = {} for p in range(len(parlist)): chain[parlist[p]] = np.array([s for tmp in data for s in tmp[p]]) # chain[parlist[p]] = np.concatenate( [tmp[p] for tmp in data]) return( chain) def mergeMC( filename, burnin=0.9, num=None, ext="txt", par=None, nelts=None): import subprocess import os import sys if num == None: num=[n for n in range(100) if os.path.isfile("%s%d.%s" % (filename,n,ext))] if par == None: par = MCparnames( filename, num[0], ext) print(num) chain = [] for name in ["%s%d.%s" % (filename,n,ext) for n in num]: f = open( name, 'rb') f.readline() #labels line_size=np.mean([len(f.readline()) for i in range(100)]) # print( line_size) nsamp = os.stat(name).st_size / line_size # print(nsamp) if 0 < burnin < 1: ist = int((1.-burnin)*nsamp) elif burnin < 0: ist = int(-burnin) else: ist = int(nsamp-burnin) if int(line_size*ist*1.2) < os.path.getsize(name): f.seek(-int(line_size*ist*1.2),os.SEEK_END) lines_found=f.readlines() f.close() # data=[np.array(x.split(),float) for x in lines_found[-ist:]] data=[[float(c) for c in x.split()] for x in lines_found[-ist:]] if len(chain)==0: chain = data else: chain = np.concatenate((chain,data),axis=0) if nelts: if nelts < len(chain): chain = chain[np.array(np.random.random(nelts)*len(chain),int),:] dchain = dict( list(zip(par,np.transpose(chain)))) return(dchain) #Read parameter names in MCMC sample file def MCparnames( filename, num=1, ext="txt"): f = open( "%s%d.%s" % (filename,num,ext)) par = [v for v in f.readline().replace("\n","").split("\t") if v != ""] f.close() return(par) def Parnames( parfile, par_type=None): param_list = [] if not os.path.isfile(parfile): raise ValueError( "Filename do not exists: %s" % parfile) sf = open(parfile) par_type=[par_type] if par_type else ["par","fix"] for l in sf.readlines(): fval = re.findall(r'\S+',l) # print(fval) if len(fval) >1 : if fval[0] in par_type: if not fval[1] in param_list: param_list.append(fval[1]) sf.close() return( param_list) def Interval( chain, par, CL=68, symmetrical=False): return( getci( np.percentile( chain[par], [(100-CL)/2,50,100-(100-CL)/2]), symmetrical=symmetrical)) def Upper( chain, par, CL=95): return( np.percentile( chain[par], CL)) def scan(filename): f = open( filename, 'r') l=f.readlines() f.close() index=0 while not l[index].startswith("SCANNING PAR"): index = index+1 print(l[index].replace("\n","")) while not l[index].startswith("it "): index = index+1 par = [v.strip() for v in l[index].replace("\n","").split("|") if v != ""] index = index + 1 scan = [] while not l[index].startswith("===="): scan.append( [v for v in l[index].replace("\n","").split("\t") if v != ""]) index = index + 1 return( scan,par) #expand given mat with the given parameters on the diag def expand_cov( mat, par, newpar, newvar): match = [newpar.index(i) for i in newpar if i in par] notmatch = [newpar.index(i) for i in newpar if i not in par] sz = len(newpar) newmat = np.zeros( (sz, sz)) for i in match: for j in match: newmat[i,j] = mat[par.index(newpar[i]),par.index(newpar[j])] for i in notmatch: newmat[i,i] = newvar[i] return( newmat) def extract_chains( chainame, nchain): chain = [] for c in range(0,nchain): print( "%s%d.txt" % (chainame,c+1)) data = [] with open("%s%d.txt" % (chainame,c+1)) as f: f.readline() for line in f: data.append( [float(a) for a in line.split()]) chain.append( np.array(np.transpose(data))) return(chain) def getci( c, symmetrical=False): if symmetrical: return(c[1], np.mean([abs(c[0]-c[1]), abs(c[2]-c[1])])) else: return(c[1], c[0]-c[1], c[2]-c[1]) def ci2sig( val, ci, symmetrical=False): if symmetrical: return(val, np.mean([abs(ci[0]-val), abs(ci[1]-val)])) else: return(val, ci[0]-val, ci[1]-val) def sig2ci( val, sig): if len(sig) == 1: return( val, (val-sig,val+sig)) else: return( val, (val-np.abs(sig[0]),val+sig[1])) def cov2cor( mat): cor = mat.copy() for i in range( len(mat)): for j in range( len(mat)): cor[i,j] = mat[i,j]/np.sqrt(mat[i,i]*mat[j,j]) return(cor) ######################################################################## #CONVERGENCE ######################################################################## def GelmanRubin( chains, gap=10000, length_min=1000, length_max=None, new=False): #plot GelmanRubin test skipping first parameter=chi2 nchain = len(chains) nsamples = min( [len(c[0]) for c in chains]) if length_max==None: length_max = nsamples if length_max > nsamples: print( "Not enough samples...") exit() it = list(range(length_min,length_max,gap))[1:-1] R = [] for isamp in it: print( "%d%%" % int(it.index(isamp)*100./len(it))) if new: #Gelman-Rubin stat on the last "gap" samples for each iteration n = gap tmp = [data[1:,isamp:isamp+gap] for data in chains] else: #Gelman-Runbin stat on the last half samples for each iteration #SP fix for python3: add int n = int((isamp - length_min)/2) tmp = [data[1:,length_min+n:isamp] for data in chains] #within Chain mchain = np.mean( tmp,2) schain = np.var( tmp,2, ddof=1) #between chains W = np.mean( schain,0) B = np.var( mchain,0, ddof=1)*n #sqrt ? R.append( ( (n-1.)/n*W + B/n )/W ) return(it,np.array(R)) def Geweke( chains, gap=10000, length_max=None): length_min = 1000 nchain = len(chains) nsamples = min( [len(c[0]) for c in chains]) if length_max==None: length_max = nsamples if length_max > nsamples: print( "Not enough samples...") exit() it = list(range(length_min,length_max,gap))[1:] T = [] for isamp in it: print( "%d%%" % int(it.index(isamp)*100./len(it))) n = length_max-isamp tmp1 = [data[:,isamp:isamp+n*0.1] for data in chains] tmp2 = [data[:,isamp+n*0.5:length_max] for data in chains] #Tscore Z1 = np.mean( tmp1,2) Z2 = np.mean( tmp2,2) s = np.sqrt( np.var(tmp1,2)/(n*0.1) + np.var(tmp2,2)/(n*0.5) ) T.append( (Z1-Z2)/s) return(it,T) ######################################################################## ######################################################################## #PARAMETER NAMES ######################################################################## parname = { 'omega_b' : '$\\Omega_\mathrm{b}h^2$', 'omega_cdm' : '$\\Omega_\mathrm{c}h^2$', 'Omega_b' : '$\\Omega_\mathrm{b}$', 'Omega_cdm' : '$\\Omega_\mathrm{c}$', '100*theta_s' : '$100\\theta_\mathrm{s}$', 'tau_reio' : '$\\tau$', 'z_reio' : '$\mathrm{z}_\mathrm{reio}$', 'n_s' : '$\mathrm{n}_\mathrm{s}$', 'log(10^10A_s)' : '$log(10^{10}A_\mathrm{s})$', 'H0' : '$H_0$', 'sigma8' : '$\\sigma_8$', 'Alens' : '$\mathrm{A}_\mathrm{L}$', 'm_ncdm' : '$\\Sigma m_\\nu$', 'N_eff' : '$N_\mathrm{eff}$', 'r' : '$r$', 'A_planck' : '$A_\mathrm{planck}$', 'A_cib_217' : '$A_\mathrm{CIB}^{217}$', 'xi_sz_cib' : '$\\xi_\mathrm{SZxCIB}$', 'A_sz' : '$A_\mathrm{SZ}$', 'ps_A_100_100' : '$A_\mathrm{PS}^\mathrm{Planck}\mathrm{(100)}$', 'ps_A_143_143' : '$A_\mathrm{PS}^\mathrm{Planck}\mathrm{(143)}$', 'ps_A_217_217' : '$A_\mathrm{PS}^\mathrm{Planck}\mathrm{(217)}$', 'ps_A_143_217' : '$A_\mathrm{PS}^\mathrm{Planck}\mathrm{(143x217)}$', 'ksz_norm' : '$A_\mathrm{kSZ}$', 'gal545_A_100' : '$A_\mathrm{gal545}^\mathrm{100}$', 'gal545_A_143' : '$A_\mathrm{gal545}^\mathrm{143}$', 'gal545_A_143_217' : '$A_\mathrm{gal545}^\mathrm{143x217}$', 'gal545_A_217' : '$A_\mathrm{gal545}^\mathrm{217}$', 'calib_100T' : '$c^\mathrm{Planck}_{100}$', 'calib_217T' : '$c^\mathrm{Planck}_{217}$', 'c0' : '$c^\mathrm{Planck}_\mathrm{100hm1}$', 'c1' : '$c^\mathrm{Planck}_\mathrm{100hm2}$', 'c2' : '$c^\mathrm{Planck}_\mathrm{143hm1}$', 'c3' : '$c^\mathrm{Planck}_\mathrm{143hm2}$', 'c4' : '$c^\mathrm{Planck}_\mathrm{217hm1}$', 'c5' : '$c^\mathrm{Planck}_\mathrm{217hm2}$', 'Aps100x100' : '$A_\mathrm{PS}^\mathrm{Planck}\mathrm{(100x100)}$', 'Aps100x143' : '$A_\mathrm{PS}^\mathrm{Planck}\mathrm{(100x143)}$', 'Aps100x217' : '$A_\mathrm{PS}^\mathrm{Planck}\mathrm{(100x217)}$', 'Aps143x143' : '$A_\mathrm{PS}^\mathrm{Planck}\mathrm{(143x143)}$', 'Aps217x217' : '$A_\mathrm{PS}^\mathrm{Planck}\mathrm{(217x217)}$', 'Aps143x217' : '$A_\mathrm{PS}^\mathrm{Planck}\mathrm{(143x217)}$', 'Asz' : '$A_\mathrm{SZ}$', 'Acib' : '$A_\mathrm{CIB}$', 'Aksz' : '$A_\mathrm{kSZ}$', 'Aszxcib' : '$A_\mathrm{SZxCIB}$', 'AdustTT' : '$A_\mathrm{dust}^\mathrm{TT}$', 'AdustTP' : '$A_\mathrm{dust}^\mathrm{TE}$', 'AdustPP' : '$A_\mathrm{dust}^\mathrm{EE}$', 'SPT_ADust' : '$A^\mathrm{dust}_\mathrm{SPT}$', 'SPT_high_95_cal' : '$c^\mathrm{SPT\ high}_{95}$', 'SPT_high_150_cal' : '$c^\mathrm{SPT\ high}_{150}$', 'SPT_high_220_cal' : '$c^\mathrm{SPT\ high}_{220}$', 'SPT_high_Aps_95x95' : '$A_{PS}^\mathrm{SPT\ high}\mathrm{(95x95)}$', 'SPT_high_Aps_95x150' : '$A_{PS}^\mathrm{SPT\ high}\mathrm{(95x150)}$', 'SPT_high_Aps_95x220' : '$A_{PS}^\mathrm{SPT\ high}\mathrm{(95x220)}$', 'SPT_high_Aps_150x150' : '$A_{PS}^\mathrm{SPT\ high}\mathrm{(150x150)}$', 'SPT_high_Aps_150x220' : '$A_{PS}^\mathrm{SPT\ high}\mathrm{(150x220)}$', 'SPT_high_Aps_220x220' : '$A_{PS}^\mathrm{SPT\ high}\mathrm{(220x220)}$', 'SPT_low_Aps' : '$A_{PS}^\mathrm{SPT\ low}$', 'SPT_low_cal' : '$c^\mathrm{SPT\ low}$', 'ACT_equat_148_cal' : '$c^\mathrm{ACT\ equat}_{148}$', 'ACT_equat_220_cal' : '$c^\mathrm{ACT\ equat}_{220}$', 'ACT_equat_ADust' : '$A_\mathrm{dust}^\mathrm{ACT\ equat}$', 'ACT_equat_Aps_148x148' : '$A_{PS}^\mathrm{ACT\ equat}\mathrm{(148x148)}$', 'ACT_equat_Aps_148x220' : '$A_{PS}^\mathrm{ACT\ equat}\mathrm{(148x220)}$', 'ACT_equat_Aps_220x220' : '$A_{PS}^\mathrm{ACT\ equat}\mathrm{(220x220)}$', 'ACT_south_148_cal' : '$c^\mathrm{ACT\ south}_{148}$', 'ACT_south_220_cal' : '$c^\mathrm{ACT\ south}_{220}$', 'ACT_south_ADust' : '$A_\mathrm{dust}^\mathrm{ACT\ south}$', 'ACT_south_Aps_148x148' : '$A_{PS}^\mathrm{ACT\ south}\mathrm{(148x148)}$', 'ACT_south_Aps_148x220' : '$A_{PS}^\mathrm{ACT\ south}\mathrm{(148x220)}$', 'ACT_south_Aps_220x220' : '$A_{PS}^\mathrm{ACT\ south}\mathrm{(220x220)}$', 'chi2' : '$\\chi^2$' } to_cosmomcName = { 'omega_b' : 'omegabh2', 'omega_cdm' : 'omegach2', '100*theta_s' : 'theta', 'tau_reio' : 'tau', 'z_reio' : 'zrei', 'n_s' : 'ns', 'log(10^10A_s)' : 'logA', 'A_planck' : 'calPlanck', 'A_cib_217' : 'acib217', 'xi_sz_cib' : 'xi', 'A_sz' : 'asz143', 'ps_A_100_100' : 'aps100', 'ps_A_143_143' : 'aps143', 'ps_A_217_217' : 'aps217', 'ps_A_143_217' : 'aps143217', 'ksz_norm' : 'aksz', 'gal545_A_100' : 'kgal100', 'gal545_A_143' : 'kgal143', 'gal545_A_143_217' : 'kgal143217', 'gal545_A_217' : 'kgal217', 'calib_100T' : 'cal0', 'calib_217T' : 'cal2', 'c0' : 'c0', 'c1' : 'c1', 'c2' : 'c2', 'c3' : 'c3', 'c4' : 'c4', 'c5' : 'c5', 'Aps100x100' : 'Aps100x100', 'Aps100x143' : 'Aps100x143', 'Aps100x217' : 'Aps100x217', 'Aps143x143' : 'Aps143x143', 'Aps217x217' : 'Aps217x217', 'Aps143x217' : 'Aps143x217', } from_cosmomcName = dict(list(zip(list(to_cosmomcName.values()),list(to_cosmomcName.keys())))) ######################################################################## #PLOTS ######################################################################## def posterior1d( chains, params, nbin=50, smooth=1, names=[], colors=[], linestyles=[], lw=[], extent=[], subplot=[1,1], figsize=None): nchain = len(chains) npar = len(params) fig=plt.figure(figsize=figsize) if len(colors) < nchain: # ax = plt.gca() # color_cycle = ax._get_lines.prop_cycler #ax._get_lines.color_cycle # colors = [next(color_cycle) for i in range(nchain)] colors = plt.rcParams['axes.color_cycle'] if len(linestyles) < nchain: linestyles = ['-']*nchain if np.size(lw) == 1: lw = [lw]*nchain elif np.size(lw) < nchain: lw = [1]*nchain if np.prod(subplot) < npar: subplot[1] = np.ceil(np.double(npar)/subplot[0]) if type(extent) == tuple: extent = [extent]*npar if len(extent) != npar: extent=None for par in params: ax=plt.subplot( subplot[0], subplot[1], params.index(par)+1) for c in range(nchain): chain = chains[c] if par in chain: hist1d( chain[par], smooth=smooth, color=colors[c], linestyle=linestyles[c], lw=lw[c], bins=nbin) ax.tick_params(axis='both', which='major', labelsize=16) ax.locator_params(tight=True, nbins=6) plt.xlabel( parname.get(par,par)) if extent: ax.set_xlim( extent[params.index(par)]) #legend if names != []: # if len(params) >= 2: if len(params) ==1: ax.legend(names, loc='upper center') else: if len(params) == subplot[0]*subplot[1]: ax.legend(names, loc='upper center', frameon=False, ncol=len(names), bbox_to_anchor=(0.,1.2)) plt.subplots_adjust( top=0.85) else: ax.legend(names, bbox_to_anchor=(1.05, 0), loc='lower left', borderaxespad=0.) plt.subplots_adjust( left=0.1, right=0.85) # ax.legend(names, loc='upper right') def posterior2d( chains, par1, par2, *args, **kwargs): """ Plot a 2-D histogram of samples. """ import scipy.stats import matplotlib.cm as cm x = chain[par1] y = chain[par2] # ax = kwargs.pop("ax", plt.gca()) extent = kwargs.pop("extent", [[x.min(), x.max()], [y.min(), y.max()]]) bins = kwargs.pop("bins", 50) colors = kwargs.pop("colors", None) alpha = kwargs.pop("alpha", 0.8) linewidths = kwargs.pop("linewidths", None) fill = kwargs.get("fill", True) datapoints = kwargs.get("datapoints", False) contours = kwargs.get("contours", False) levels = kwargs.get("levels", [0.68,0.95]) cmap = kwargs.get("cmap", None) # cmap = cm.get_cmap(kwargs.get("cmap", None)) X = np.linspace(extent[0][0], extent[0][1], bins + 1) Y = np.linspace(extent[1][0], extent[1][1], bins + 1) try: H, X, Y = np.histogram2d(x.flatten(), y.flatten(), bins=(X, Y), weights=kwargs.get('weights', None)) except ValueError: raise ValueError("It looks like at least one of your sample columns " "have no dynamic range. You could try using the " "`extent` argument.") if datapoints: color = 'k' if colors==None else colors[0] plt.plot(x, y, "o", color=color, ms=1.5, zorder=-1, alpha=0.1, rasterized=True) gkde = scipy.stats.gaussian_kde([x, y]) dx = (extent[0][1]-extent[0][0])/(bins+1.) dy = (extent[1][1]-extent[1][0])/(bins+1.) x2,y2 = np.mgrid[extent[0][0]:extent[0][1]:dx, extent[1][0]:extent[1][1]:dy] z = np.array(gkde.evaluate([x2.flatten(),y2.flatten()])).reshape(x2.shape) W = np.append(ctr_level( z, levels),z.max()) if fill: plt.contourf(x2, y2, z, W, alpha=alpha, cmap=cm.get_cmap(cmap), antialiased=True, colors=colors) if contours: plt.contour(x2, y2, z, W,linewidths=linewidths, colors=colors, antialiased=False) plt.xlim(extent[0]) plt.ylim(extent[1]) plt.xlabel(parname.get(par1,par1)) plt.ylabel(parname.get(par2,par2)) #backward compatibility def hist2d( chain, par1, par2, *args, **kwargs): """ Plot a 2-D histogram of samples. """ import scipy.stats import matplotlib.cm as cm x = chain[par1] y = chain[par2] extent = kwargs.pop("extent", [[x.min(), x.max()], [y.min(), y.max()]]) bins = kwargs.pop("bins", 50) nsmooth = kwargs.pop("nsmooth", 0.) cmap = kwargs.get("cmap", None) X = np.linspace(extent[0][0], extent[0][1], bins + 1) Y = np.linspace(extent[1][0], extent[1][1], bins + 1) try: H, X, Y = np.histogram2d(x.flatten(), y.flatten(), bins=(X, Y), weights=kwargs.get('weights', None)) except ValueError: raise ValueError("It looks like at least one of your sample columns " "have no dynamic range. You could try using the " "`extent` argument.") if nsmooth != 0.: H = nd.gaussian_filter( H, nsmooth) npts = 5. # plt.imshow( H.T, origin='bottom', **kwargs) # plt.xticks(np.arange(npts+1)/npts*bins,["%4.2f" % i for i in np.linspace(extent[0][0],extent[0][1],npts+1)]) # plt.yticks(np.arange(npts+1)/npts*bins,["%4.2f" % i for i in np.linspace(extent[1][0],extent[1][1],npts+1)]) plt.pcolormesh( X, Y, H.T, **kwargs) plt.xlabel(parname.get(par1,par1)) plt.ylabel(parname.get(par2,par2)) def triangle( chains, params, nbin=50, parnames=None, names=[], colors=[], linestyles=[], smooth=1., fontsize=12, contour=True, fill=False, cmaps=[], bestfit=None, Norm=True, weights=None): import matplotlib.ticker as mtick fig=plt.figure(figsize=(14,12)) plt.subplots_adjust(hspace=0.001,wspace=0.001) #force tuple if not( isinstance(chains, list) or isinstance(chains, tuple)): chains = (chains,) nchain = len(chains) npar = len(params) if len(colors) < nchain: # ax = plt.gca() # color_cycle = ax._get_lines.prop_cycler #ax._get_lines.color_cycle # colors = [next(color_cycle) for i in range(nchain)] colors = plt.rcParams['axes.color_cycle'] if len(cmaps) < nchain: cmaps = ['Blues','Greens','Greys','Oranges','Purples'] # cmaps = [m for m in cm.datad if not m.endswith("_r")] if type(fill) == bool: fill = [fill]*nchain if type(contour) == bool: contour = [contour]*nchain if len(linestyles) < nchain: linestyles = ['-']*nchain if weights == None: weights = [np.ones(len(c['chi2'])) for c in chains] if parnames == None: parnames = [parname.get(p,p) for p in params] rge = {} for par in params: tmp = [] for chain in chains: tmp = tmp+list(chain[par]) rge[par] = list([min(tmp),max(tmp)]) # rge[par] = list(np.mean(tmp) + np.dot(np.std(tmp), [-4,4])) for par in params: xpar = params.index(par) ax=plt.subplot( npar, npar, xpar*npar+xpar+1) for c in range(nchain): chain = chains[c] weight = weights[c] h,X=np.histogram( chain[par], bins=nbin, range=rge[par], normed=True, weights=weight) Y = nd.gaussian_filter1d( h, smooth,mode='nearest') Y = Y/np.max(Y) if Norm else Y plt.plot( (X[1:]+X[:-1])/2., Y, color=colors[c], linestyle=linestyles[c]) if bestfit: plt.axvline( (bestfit[c])[par], ls='--', color=colors[c]) ax.set_xlim( rge[par]) ax.yaxis.set_visible(False) ax.xaxis.set_visible(xpar==(npar-1)) ax.locator_params(tight=True, nbins=5) if 'PS' in parname.get(parnames[xpar],parnames[xpar]): ax.xaxis.set_major_formatter(mtick.FormatStrFormatter('%.0e')) #ax.xaxis.set_major_formatter(mtick.FormatStrFormatter('%.2e')) for ypar in range(xpar)[::-1]: par2 = params[ypar] ax=plt.subplot( npar, npar, xpar*npar+ypar+1) for c in range(nchain): chain = chains[c] weight = weights[c] hmap,binX,binY = np.histogram2d( chain[par], chain[par2], bins=nbin, range=[rge[par],rge[par2]], weights=weight) Y,X = np.meshgrid((binX[1:]+binX[:-1])/2.,(binY[1:]+binY[:-1])/2.) hmap = nd.gaussian_filter( hmap.T, 2) if fill[c]: plt.contourf( X,Y, hmap, levels=ctr_level( hmap, [1e-2,0.68,0.95]), colors=None, extent=tuple(rge[par]+rge[par2]), cmap=plt.cm.get_cmap(cmaps[c])) if contour[c]: plt.contour( X,Y, hmap, levels=ctr_level( hmap, [0.68,0.95]), colors=colors[c], linestyles=linestyles[c], extent=tuple(rge[par]+rge[par2])) if bestfit: plt.plot( (bestfit[c])[par2], (bestfit[c])[par], '+', color=colors[c]) # ax.locator_params( nbins=npar-1) ax.ticklabel_format(useOffset=False) # ax.tick_params(axis='both', which='major', labelsize=10) # ax.tick_params(axis='both', which='minor', labelsize=5) ax.yaxis.set_visible(ypar==0) ax.xaxis.set_visible(xpar==(npar-1)) ax.locator_params(tight=True, nbins=5) if 'PS' in parname.get(parnames[ypar],parnames[ypar]): ax.xaxis.set_major_formatter(mtick.FormatStrFormatter('%.0e')) #plot parameter names for xpar in range(npar): ax=plt.subplot( npar, npar, xpar*npar+1) ax.set_ylabel(parname.get(parnames[xpar],parnames[xpar])) ax.tick_params(axis='both', which='major', labelsize=fontsize) ax=plt.subplot( npar, npar, (npar-1)*npar+xpar+1) ax.set_xlabel(parname.get(parnames[xpar],parnames[xpar])) ax.tick_params(axis='both', which='major', labelsize=fontsize) # plt.locator_params(tight=True, nbins=5) #write legend if len(names) == nchain: ax=plt.subplot( npar, npar, int(0.25*npar)*npar+0.85*npar-1) for c in range(nchain): plt.text(0.5,1.-0.2*(c+1), names[c], color=colors[c]) ax.axis('off') return(fig) def rectangle( chains, par0, params, nbin=50, parnames=None, names=[], colors=[], linestyles=[], smooth=1., fontsize=12, contour=True, fill=False, cmaps=[], Norm=True, weights=None, ylim=None): import matplotlib.ticker as mtick #force tuple if not( isinstance(chains, list) or isinstance(chains, tuple)): chains = (chains,) nchain = len(chains) npar = len(params) fig=plt.figure( figsize=(5*npar,5)) plt.subplots_adjust(hspace=0.001,wspace=0.001) if len(colors) < nchain: colors = plt.rcParams['axes.color_cycle'] if len(cmaps) < nchain: cmaps = ['Blues','Greens','Greys','Oranges','Purples'] if type(fill) == bool: fill = [fill]*nchain if type(contour) == bool: contour = [contour]*nchain if len(linestyles) < nchain: linestyles = ['-']*nchain if weights == None: weights = [np.ones(len(c['chi2'])) for c in chains] if parnames == None: parnames = [parname.get(p,p) for p in params] if ylim == None: tmp = [] for chain in chains: tmp = tmp+list(chain[par0]) ylim = list(np.mean(tmp) + np.dot(np.std(tmp), [-4,4])) rge = {} for par in params: tmp = [] for chain in chains: tmp = tmp+list(chain[par]) rge[par] = list(np.mean(tmp) + np.dot(np.std(tmp), [-4,4])) for par in params: xpar = params.index(par) ax=plt.subplot( 1, npar, xpar+1) for c in range(nchain): chain = chains[c] weight = weights[c] hmap,binX,binY = np.histogram2d( chain[par], chain[par0], bins=nbin, range=[rge[par],ylim], weights=weight) X,Y = np.meshgrid((binX[1:]+binX[:-1])/2.,(binY[1:]+binY[:-1])/2.) hmap = nd.gaussian_filter( hmap.T, 2) if fill[c]: plt.contourf( X,Y, hmap, levels=np.append(list(ctr_level( hmap, [0.68,0.95])[::-1]),hmap.max()), colors=None, extent=tuple(rge[par]+ylim), cmap=plt.cm.get_cmap(cmaps[c])) if contour[c]: plt.contour( X,Y, hmap, levels=ctr_level( hmap, [0.68,0.95]), colors=colors[c], linestyles=linestyles[c], extent=tuple(rge[par]+ylim)) ax.ticklabel_format(useOffset=False) ax.locator_params(tight=True, nbins=5) ax.set_xlabel(parname.get(parnames[xpar],parnames[xpar])) ax.tick_params(axis='both', which='major', labelsize=fontsize) ax.yaxis.set_visible(xpar==0) if xpar == 0: ax.set_ylabel(parname.get(par0,par0)) # if 'PS' in parname.get(parnames[par],parnames[par]): # ax.xaxis.set_major_formatter(mtick.FormatStrFormatter('%.0e')) #write legend if len(names) == nchain: for c in range(nchain): ax.text( 0.9, 1.-0.1*(c+1), names[c], color=colors[c], transform=ax.transAxes, horizontalalignment='right') return(fig) def hist1d(x, bins=50, weights=None, smooth=0., normmax=True, **kwargs): extents = kwargs.get("extents",[x.min(),x.max()]) n, b = np.histogram(x, weights=weights, bins=bins,normed=True) ns = nd.gaussian_filter1d(n,smooth,mode='nearest') if normmax: ns /= ns.max() xx = (b[:-1]+b[1:])/2. yy = ns plt.plot( xx, yy, **kwargs) plt.xlim( extents) def MCcorrelation( chain, par, **kwargs): data = [chain[p] for p in par] matC = np.corrcoef( data) plot_correlation( matC, par, **kwargs) def plot_correlation( matC, par, **kwargs): vmin = kwargs.pop("vmin", -1) vmax = kwargs.pop("vmax", 1) plt.figure( figsize=(14,10)) #covmat = np.triu(np.corrcoef( chain.values()[1:]), k=1) covmat = np.tril( matC, k=-1) covmat[covmat==0] = np.NaN plt.imshow( covmat, vmin=vmin, vmax=vmax, **kwargs) plt.yticks( list(range(len(covmat))), [parname.get(p,p) for p in par]) plt.xticks( list(range(len(covmat))), '') plt.colorbar() for i in range( len(par[:-1])): # text( i, i+0.25, parname.get(par[i+1],par[i+1]),rotation=-90) plt.text( i-0.25, i+0.25, parname.get(par[i],par[i]),rotation='vertical',verticalalignment='bottom') ######################################################################## ######################################################################## def ctr_level(histo2d, lvl): """ Extract the contours for the 2d plots """ h = histo2d.flatten()*1. h.sort() cum_h = np.cumsum(h[::-1]) cum_h /= cum_h[-1] alvl = np.searchsorted(cum_h, lvl) clist = h[-alvl] return clist[::-1] def FeldmanCousins( xmin, sigma, CL=95): import os from scipy.interpolate import interp1d fcdata = np.loadtxt( os.getenv("CAMELROOT")+"/work/tools/FeldmanCousins.data").T listCL = [68,90,95,99] if not(CL in listCL): raise ValueError('CL not in the list') clmin = interp1d( fcdata[0], fcdata[listCL.index(CL)*2+1]) clmax = interp1d( fcdata[0], fcdata[listCL.index(CL)*2+2]) x0=xmin/sigma uplim = clmax(x0) return( uplim*sigma) ######################################################################## ######################################################################## def read_bestfit( dirname): if '.gz' in dirname : f= gzip.open(dirname , 'rb') else : f = open( dirname, 'r') param = np.array( f.readline().split(), dtype=str) newline = f.readline() if not "nan" in newline and not "inf" in newline and not "Abort" in newline: value = np.array( newline.split(), dtype=np.double) bestfit = dict(list(zip(param, value))) else : bestfit = {'chi2':np.inf} f.close() return( bestfit) def extract_bestfit( filename, skipone=False): f = open( filename, 'r') lines=f.readlines() f.close() index = 0 while not "# ext." in lines[index]: index = index+1 if skipone: index = index+1 while not "# ext." in lines[index]: index = index+1 print(index) index = index+2 param = [] value = [] error = [] while "||" in lines[index]: i,p,t,v,e = lines[index].split("||") if t.strip() != "const": param.append(p.strip()) value.append(float(v.replace("*","").strip())) error.append(float(e.replace("*","").strip())) index = index + 1 return(param,value,error) def extract_hess( filename, ext="LASymMatrix", reverse=False): f = open( filename, 'r') lines=f.readlines() f.close() npar = 0 par = [] for l in lines: if l.startswith('par\t'): npar = npar + 1 par.append( l.split("\t")[1]) print( "number of parameter : %d/%d" % (npar,len(par))) if reverse: lines.reverse() index = 0 while not ext in lines[index]: index = index+1 print(index) # while lines[index] == "\n" or lines[index] == "": # index = index+1 mat=[] for n in range(npar): mat.append([v for v in lines[index+1+n].replace("\n","").split(" ") if v != ""]) n = len(mat[0]) matcov = np.zeros( (n,n)) for i in range(n): matcov[i,:] = np.array(mat[i]) if reverse: matcov = flipud(matcov) return(matcov,par) def gauss_prior(x,m,s): L = np.exp( -0.5*((x-m)/s)**2) return( L) ######################################################################## simple_lik_name = { 'smica_g30_ftl_full_pttptt' : 'lensing', 'lowl_SMW_70_dx11d_2014_10_03_v5c_Ap' : 'bflike', 'HiLLiPOP:DX11dHM_superExt_CO_TT' : 'hlpTT', 'Hillipop:DX11dHM_superExt_CO_TT' : 'hlpPS', 'SPT_low' : 'SPT_low', 'SPT_high' : 'SPT_high', 'ACT_equat' : 'ACT_equat', 'ACT_south' : 'ACT_south', 'commander_rc2_v1.1_l2_29_B' : 'cmder', 'Lollipop' : 'lollipop', 'lollipop' : 'lollipop', 'lollipopn' : 'lollipop', 'BAO1D_dr12' : 'BAO(1D)', 'boss_dr12' : 'BAO(3D)', 'JLA' : 'SNIa(JLA)' } def writeExtendedBestFit( dbdir, bfwrite=False, outpath=""): #Extract Bestfit from logfile path, dirs, files = next(os.walk("%s" % (dbdir))) for idir in dirs: p,d,f = next(os.walk("%s/%s" % (dbdir,idir))) # this will find the last 'camelZZ.par' file in folder camelparfile="" for s in f: if ("camel" in s and ".par" in s): camelparfile="%s/%s/%s" % (dbdir,idir,s) # on cherche si on a lance le code avec Profile ou Minimize (camelrun) CheckProfile=False name ="" for file in glob.glob("%s/%s/camelrun*" % (dbdir,idir)): name = file print(file) try: #name="%s/%s/%s" % (dbdir,idir,"camelrun") cr=open(name) all_lines=cr.readlines() for c in all_lines: if "Profile" in c and not "#" in c: CheckProfile=True cr.close() except ValueError: return() print(CheckProfile) #build a line with all param names avoiding double counts try: parlist = Parnames( camelparfile) except ValueError: return() param_line = " ".join(parlist) if CheckProfile: param_line += " profiled " param_line+=" chi2 valid edm maxlim " #print param_line read=True line_satisfied=False for s in [f1 for f1 in f if ".o" in f1]: if read: name="%s/%s/%s" % (dbdir,idir,s) if 'gz' in name : sf= gzip.open( name , 'rb') else : sf = open( name, 'r') all_lines=sf.readlines() # gets all lines in file for c in all_lines: if "=====" in c : #print c chnam = c.split(" ")[1] #print chnam if ".lik" in c : param_line += simple_lik_name[chnam.replace(".lik:","")] param_line += " " line_satisfied=True #print chnam #print "|"+param_line+"|" if "ollipop" in c : param_line += simple_lik_name[chnam.replace(":","")] param_line += " " line_satisfied=True #print chnam #print "|"+param_line+"|" if ".clik" in c : param_line += simple_lik_name[chnam.replace(".clik:","")] param_line += " " line_satisfied=True #print "|"+param_line+"|" if ".lal" in c : param_line += simple_lik_name[chnam.replace(".lal:","")] param_line += " " line_satisfied=True #print "|"+param_line+"|" if ".par" in c : param_line += simple_lik_name[chnam.replace(".par:","")] param_line += " " line_satisfied=True #print 'dans .par',chnam,chnam.replace(".par:","") if line_satisfied: read=False param_line = param_line.replace("\n"," ") + "\n" #print param_line sf.close() for s in [f1 for f1 in f if ".o" in f1]: name="%s/%s/%s" % (dbdir,idir,s) if 'gz' in name : sf= gzip.open( name , 'rb') else : sf = open( name, 'r') # gets Profile line goodline = False values ="" for c in sf.readlines(): if goodline: values += " " if (".lik" in c or ".clik" in c or ".lal" in c or ".par" in c or "llipop" in c) and "=====" in c : chval = c.split(" ")[2] values += (chval.replace("chi2=","")).replace("\n","") values += " " if "Profile" in c and not "bestfit" in c and not "nan" in c and not "inf" in c and not "Abort" in c: values = c goodline = True values = values.replace("\n","") sf.close() # gets extensions #write file if '.gz' in name : count = int(name.split(".")[-2]) else: count = int(name.split(".")[-1]) bestfitfile = "%s/%s/extbest_fit%d" % (dbdir,idir,count) if bfwrite : bestfitfile = "%s/%s/best_fit%d" % (dbdir,idir,count) if outpath != "" : print("outpath gene -> %s/" %(outpath)) #print os.path.isdir( "%s/" %(outpath)) if os.path.isdir( "%s/" %(outpath))==False : #print 'mkdir ?' #return(); os.mkdir( "%s" %(outpath) ) print("outpath -> %s/%s/" %(outpath,idir)) #print os.path.isdir( "%s/%s/" %(outpath,idir)) if os.path.isdir( "%s/%s/" %(outpath,idir))==False : #print 'mkdir ?' #return(); os.mkdir( "%s/%s" %(outpath,idir) ) bestfitfile = "%s/%s/best_fit%d" %(outpath,idir,count) # print(values,bestfitfile) # if (not os.path.isfile(bestfitfile)) : if goodline: newf=open( bestfitfile,'w') newf.write(param_line) newf.write(values.replace("Profile","")) newf.close() return() def writeDerivedBestFit( dbdir, bfwrite=False, outpath=""): #Extract Bestfit from logfile path, dirs, files = next(os.walk("%s" % (dbdir))) for idir in dirs: p,d,f = next(os.walk("%s/%s" % (dbdir,idir))) # this will find the last 'camelZZ.par' file in folder camelparfile="" for s in f: if ("camel" in s and ".par" in s): camelparfile="%s/%s/%s" % (dbdir,idir,s) # on cherche si on a lance le code avec Profile ou Minimize (camelrun) CheckProfile=False try: name="%s/%s/%s" % (dbdir,idir,"camelrun") cr=open(name) all_lines=cr.readlines() for c in all_lines: if "Profile" in c: CheckProfile=True cr.close() except ValueError: return() #build a line with all param names avoiding double counts try: parlist = Parnames( camelparfile) except ValueError: return() param_line = " ".join(parlist) if CheckProfile: param_line += " profiled " param_line+=" chi2 valid edm maxlim " # add derived params param_line += " H0 sigma8 rs_drag " read=True line_satisfied=False # add likelihoods for s in [f1 for f1 in f if ".o" in f1]: if read: name="%s/%s/%s" % (dbdir,idir,s) if 'gz' in name : sf= gzip.open( name , 'rb') else : sf = open( name, 'r') all_lines=sf.readlines() # gets all lines in file for c in all_lines: if "=====" in c : chnam = c.split(" ")[1] if ".lik" in c : param_line += simple_lik_name[chnam.replace(".lik:","")] param_line += " " line_satisfied=True #print "|"+param_line+"|" if ".clik" in c : param_line += simple_lik_name[chnam.replace(".clik:","")] param_line += " " line_satisfied=True #print "|"+param_line+"|" if ".lal" in c : param_line += simple_lik_name[chnam.replace(".lal:","")] param_line += " " line_satisfied=True #print "|"+param_line+"|" if ".par" in c : param_line += simple_lik_name[chnam.replace(".par:","")] param_line += " " line_satisfied=True #print 'dans .par',chnam,chnam.replace(".par:","") if line_satisfied: read=False param_line = param_line.replace("\n"," ") + "\n" #print param_line sf.close() for s in [f1 for f1 in f if ".o" in f1]: # selects an ouput file name="%s/%s/%s" % (dbdir,idir,s) if 'gz' in name : sf= gzip.open( name , 'rb') else : sf = open( name, 'r') print(name) if '.gz' in name : count = int(name.split(".")[-2]) else: count = int(name.split(".")[-1]) # gets camelX.par cparfile = "camel%d.par" %(count) camelparfile="%s/%s/%s" % (dbdir,idir,cparfile) cpf = open(camelparfile,"r") # CHECKS AND IF NEEDED COMPLETE THE PAR FILE xxl = mmap.mmap(cpf.fileno(), 0, access=mmap.ACCESS_READ) if not xxl.find('derived') or not xxl.find('H0') : adst ="\n derived H0" wpf = open(camelparfile,"a") wpf.write(adst) wpf.close() if not xxl.find('derived') or not xxl.find('sigma8') : adst ="\n derived sigma8" wpf = open(camelparfile,"a") wpf.write(adst) wpf.close() if not xxl.find('derived') or not xxl.find('rs_drag') : adst ="\n derived rs_drag" wpf = open(camelparfile,"a") wpf.write(adst) wpf.close() cpf.close() goodline = False gderived = False values ="" vder = "" for c in sf.readlines(): if goodline: values += " " if (".lik" in c or ".clik" in c or ".lal" in c or ".par" in c ) and "=====" in c : chval = c.split(" ")[2] values += (chval.replace("chi2=","")).replace("\n","") values += " " if not gderived: print("before os.system") os.system("$CAMELROOT/Linux-x86_64/writeSpectra %s 1000 toto.fits > tmp.out" %(camelparfile) ) addin = open("tmp.out","r") adl = addin.readlines() print("after run externe") gderived = True # gets Profile line for d in adl : if "derived" in d and "asked" in d and "H0" in d : chval = d.split(" ")[3] vder += (chval.replace("H0=","")).replace("\n","") vder += " " if "derived" in d and "asked" in d and "sigma8" in d and not "(" in d : chval = d.split(" ")[3] vder += (chval.replace("sigma8=","")).replace("\n","") vder += " " if "derived" in d and "asked" in d and "rs_drag" in d : chval = d.split(" ")[3] vder += (chval.replace("rs_drag=","")).replace("\n","") vder += " " addin.close() print("after get derived %s" %(vder)) if "Profile" in c and not "bestfit" in c and not "nan" in c and not "inf" in c and not "Abort" in c: values = c goodline = True values = values.replace("\n","") print("values ok") sf.close() # gets extensions #write file if '.gz' in name : count = int(name.split(".")[-2]) else: count = int(name.split(".")[-1]) bestfitfile = "%s/%s/extbest_fit%d" % (dbdir,idir,count) if bfwrite : bestfitfile = "%s/%s/best_fit%d" % (dbdir,idir,count) if outpath != "" : if os.path.isdir( "%s/%s/" %(outpath,idir))==False : os.mkdir( "%s/%s/" %(outpath,idir) ) bestfitfile = "%s/%s/best_fit%d" %(outpath,idir,count) print(bestfitfile) print(goodline) # print(values,bestfitfile) # if (not os.path.isfile(bestfitfile)) : if goodline: newf=open( bestfitfile,'w') newf.write(param_line) values= values + vder print(vder) print(values) print(bestfitfile) newf.write(values.replace("Profile","")) newf.close() return() ######################################################################## #PROFILE class prof: def __init__( self, dbdir, parname, nestimate=3, nameforclass=None, offset=0., writeBestfit=False, remove=[]): """ Class profile for CAMEL Parameters ---------- dbdir: directory where the minimisations runs for each value of the param profiled are stored parname: parameter name (prefix of the subfolders in dbdir) Options ------- nestimate: minimal number of minimisation per profile point nameforclass: name of parameter for CLASS if different from parname offset: offset to be applied to the likelihood values writeBestfit: extract bestfit from log file (usefull if "best-fit" files do not exist) """ import os self.dbdir = dbdir self.parname = parname self._nameforclass = nameforclass if nameforclass else self.parname self._result = False self._str_result = "" self._is_fitted = False self._offset = offset self._CL = 68 path, dirs, files = next(os.walk("%s" % (self.dbdir))) if sum([self.parname+"_" in d for d in dirs])==0: raise ValueError( "Wrong directory") self._s = [] self._v = [] for idir in dirs: if not self.parname in idir: continue p,d,f = next(os.walk("%s/%s" % (self.dbdir,idir))) if writeBestfit: self._writeBestFit(idir,f) #check there are more than nestimate bestfit in directory idir if sum([s.startswith("best_fit") for s in f]) >= nestimate: # gets the value indicated in the subfolder's name keep=True for ii in remove: if ii == idir.replace(self.parname+'_',''): keep=False if keep: # values in string self._s.append( idir.replace(self.parname+'_','')) # values in double + offset self._v.append( np.double(idir.replace(self.parname+'_',''))+self._offset) if len(self._s) <= 2: raise ValueError( "not enough profiled-value available:\n %s/%s" % (self.dbdir,idir)) sortidx = np.argsort(self._v) self._s = [self._s[s] for s in sortidx] self._v = [self._v[s] for s in sortidx] #get profile v,p = self.get() self._p = p self.ymin = min(p) def write( self, filename): """ Write the profile in ASCII file """ params = [self.parname]+list(np.extract( self._params(self._s[0])!=self.parname, self._params(self._s[0]))) param_line = " ".join( params)+"\n" f=open( filename,'w') f.write(param_line) # print(param_line) for val in self._s: data = self._point_minimum( val) # print( " ".join([str(data[p]) for p in params])) print(params) f.write( " ".join([str(data[p]) for p in params])+"\n") f.close() def _writeBestFit(self, dir, f): #Extract Bestfit from logfile # this will find the last 'camelZZ.par' file in folder camelparfile="" for s in f: if ("camel" in s and ".par" in s): camelparfile="%s/%s/%s" % (self.dbdir,dir,s) #build a line with all param names avoiding double counts try: parlist = Parnames( camelparfile) except ValueError: return() param_line = " ".join(parlist)+" chi2 valid edm maxlim \n" # param_line = " ".join(parlist)+" chi2 valid maxlim edm\n" count=0 for s in [f1 for f1 in f if ".o" in f1]: name="%s/%s/%s" % (self.dbdir,dir,s) # gets Profile line sf = open(name) for c in sf.readlines(): if "Profile" in c and not "bestfit" in c and not "nan" in c and not "inf" in c and not "Abort" in c: values = c sf.close() #write file count = int(name.split(".")[-1]) bestfitfile = "%s/%s/best_fit%d" % (self.dbdir,dir,count) # print(values,bestfitfile) # if not os.path.isfile(bestfitfile): if values: newf=open( bestfitfile,'w') newf.write(param_line) newf.write(values.replace("Profile","")) newf.close() def point_allfits( self, val, par='chi2'): """ Get all fit chi2 at a given value from multi-Minimize runs """ import os path, dirs, files = next(os.walk("%s/%s_%s" % (self.dbdir,self.parname,val))) files = [f for f in files if f.startswith("best")] nmin = len(files) allfits = [] for f in np.sort(files): #read all parameters fname= "%s/%s_%s/%s" % (self.dbdir,self.parname,val,f) if os.stat(fname).st_size !=0: bftmp = read_bestfit( "%s/%s_%s/%s" % (self.dbdir,self.parname,val,f)) allfits.append( bftmp[par]) return( allfits) def _point_minimum( self, val, verbose=False): """ Get minimum at a given value from multi-Minimize runs """ import os path, dirs, files = next(os.walk("%s/%s_%s" % (self.dbdir,self.parname,val))) files = [f for f in files if f.startswith("best")] nmin = len(files) bf = {'chi2':np.inf} for f in np.sort(files): #read all parameters fname= "%s/%s_%s/%s" % (self.dbdir,self.parname,val,f) if verbose: print(fname) print(os.stat(fname).st_size) if os.stat(fname).st_size !=0: bftmp = read_bestfit("%s/%s_%s/%s" % (self.dbdir, self.parname, val, f)) if verbose: prtval = 'chi2' if verbose == True else verbose print( f, bftmp[prtval]) if bftmp['chi2'] < bf['chi2']: bf = bftmp return( bf) def _params( self, val): #Returns list of parameters from the fit file for a given value of the parameter import os path, dirs, files = next(os.walk("%s/%s_%s" % (self.dbdir,self.parname,val))) files = [f for f in files if f.startswith("best")] #read all parameters fname= "%s/%s_%s/%s" % (self.dbdir,self.parname,val,files[0]) f = open( fname, 'r') params = np.array( f.readline().split(), dtype=str) return( params) def get( self, par='chi2'): """ Return profile likelihood Option ------ par: str [default='chi2'] parameter to return Return ------ v, p: list list of profiled parameter value and profile value """ p = [] for val in self._s: p.append( self._point_minimum( val)[par]) return( self._v, p) def fit( self, method='cubic', order=3, deltaMax=np.inf, upper=False, CL=None, bounds=None): """ Fit the profile with different models Options ------- method: str [default='cubic'] 'cubic': 1D cubic interpolation 'spline': One-dimensional smoothing spline 'polyfit': Least squares polynomial fit order: int [default=3] order of polynomial fit (if demanded) deltaMax: float fit only points for which deltaChi2 < deltaMax upper: bool [default=False] derives upper limit with feldman cousins CL: int confidence level of the limit Returns ------- (xmin, (xlow,xhigh)): by default X: X% CL confidence level if "upper=True" """ from scipy.interpolate import interp1d,UnivariateSpline,InterpolatedUnivariateSpline if CL: self._CL = CL else: self._CL = 95 if upper else 68 datav = np.array(self._v)[(self._p-min(self._p)) < deltaMax] datap = np.array(self._p)[(self._p-min(self._p)) < deltaMax] if method=='spline': # self.f = UnivariateSpline(datav,datap) self.f = InterpolatedUnivariateSpline( datav, datap, k=order) elif method=='cubic': self.f = interp1d(datav,datap, kind='cubic') elif method=='polyfit': param = np.polyfit(datav,datap,order) self.f = np.poly1d(param) else: print( "method not supported...") print( "(methods available: spline, cubic, polyfit)") return() # raise ValueError( "method not supported...") if not bounds: bounds = (min(self._v),max(self._v)) xmin = self._minimum(bounds) if upper: xlow,xhigh = self._deltaChi2( CL=68) self._res = FeldmanCousins( xmin, xhigh-xmin, self._CL) self._str_res="%s < %f (%d%% C.L.)" % (parname.get(self._nameforclass,self._nameforclass),self._res, self._CL) print(self._str_res) else: xlow,xhigh = self._deltaChi2( CL=self._CL) self._res = (xmin,(xlow,xhigh)) self._str_res="$%f_{%+f}^{%+f}$" % (xmin,xlow-xmin,xhigh-xmin) print(self._str_res) self._is_fitted = True return( self._res) def get_fit(self): """ Return the fitted function """ if self._is_fitted: return(self.f) else: print( "Please fit before...") return() def get_result(self): """ Return the result from the fit """ print( self._str_res) return(self._res) def _minimum( self, bounds=None): #set and return minimum chi2 (xmin and ymin) from scipy.optimize import minimize_scalar if bounds: res=minimize_scalar(self.f,bounds=bounds,method='bounded') else: res=minimize_scalar(self.f) self._xmin=res.x self.ymin=self.f(self._xmin) return( self._xmin) def _deltaChi2( self, delta=1, CL=None): #set and return param values [xlow,xhigh] for deltaChi = delta #based on the fitted profile from scipy.optimize import brentq #-> returns zeros of a function if CL: CL2delta = {68:1,95:3.84,99:6.63} if CL in list(CL2delta.keys()): delta = CL2delta[CL] else: raise ValueError('CL not in the list') # function chi2-min(chi2)-delta fm1 = np.vectorize(lambda x: self.f(x)-self.ymin-delta) #position of chi2min+delta below minimum self._xlow = -np.inf if( self._xmin > self._v[0]): try: self._xlow = brentq(fm1,self._v[0],self._xmin) except ValueError: self._xlow = -np.inf #position of chi2min+delta above minimum try: self._xhigh = brentq(fm1,self._xmin,self._v[-1]) except ValueError: self._xhigh = np.inf # returns 68% confidence interval return( self._xlow, self._xhigh) def plot( self, draw_CL=True, result=True, extent=None,fontsize=18, **kwargs): """ Plot the profile Options ------- draw_CL: bool draw confidence limit as vertical lines result: bool print result on the plot extent: [] min,max of the x-axis kwargs: all plot arguments """ ax = plt.gca() from matplotlib import __version__ if __version__ >= '1.5': colordef = next(ax._get_lines.prop_cycler)['color'] else: colordef = next(ax._get_lines.color_cycle) color = kwargs.pop( "color", colordef) linestyle = kwargs.pop( "linestyle", '-') label = kwargs.pop( "label", None) marker = kwargs.pop( "marker", 'o') if extent==None: extent = (self._v[0],self._v[-1]) #plot profile refmin = self.ymin if self._is_fitted: #use min(ftopot) to set the y shift vartoplot=np.linspace(*extent, num=1000) ftoplot=self.f(vartoplot) refmin = min(ftoplot) baseplt,= ax.plot(self._v,self._p-refmin,color=color,marker=marker,linestyle='') ax.set_xlabel(parname.get(self._nameforclass,self._nameforclass),fontsize=fontsize) ax.set_ylabel("$\chi^2-\chi^2_{min}$",fontsize=fontsize) ax.set_ylim((-0.1,5)) color = baseplt.get_color() if draw_CL and self._CL==68: plt.axhline(y=1,color='k',linestyle='--') #plot fit if self._is_fitted: ax.plot(vartoplot,ftoplot-min(ftoplot), color=color, linestyle=linestyle, label=label, **kwargs) CI = (self._xlow,self._xhigh) if draw_CL: if np.isinf(CI[0]): res = FeldmanCousins( self._xmin, self._xhigh-self._xmin, self._CL) ax.plot([res,res],[-0.1,1], color=color, **kwargs) else: ax.plot([CI[0],CI[0]],[-0.1,1], color=color, **kwargs) ax.plot([CI[1],CI[1]],[-0.1,1], color=color, **kwargs) if result: if np.isinf(CI[0]): ax.text(max([self._xmin,0.]), 4, ' '+self._str_res, horizontalalignment='left' ) else: ax.text(self._xmin, 4, self._str_res, horizontalalignment='center') ########################################################################