Skip to content

GitLab

Commit d1ac5cd8 authored by Plaszczynski Stephane's avatar Plaszczynski Stephane
Browse files

Merge branch '11-add-python-packaging' into 'dev' 

Resolve "camel.py"

See merge request cosmotools/CAMEL!7
parents 1a9f4c6b baefbd21
Hide whitespace changes
Inline Side-by-side
Showing with 185 additions and 184 deletions
work/tools/python/camel.py → work/tools/python/camel/camel.py
View file @ d1ac5cd8
......@@ -46,7 +46,7 @@ def mergeMC( filename, burnin=0.9, num=None, ext="txt", par=None, nelts=None):
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:
......@@ -55,15 +55,15 @@ def mergeMC( filename, burnin=0.9, num=None, ext="txt", par=None, nelts=None):
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:
......@@ -72,8 +72,8 @@ def mergeMC( filename, burnin=0.9, num=None, ext="txt", par=None, nelts=None):
if nelts:
if nelts < len(chain):
chain = chain[np.array(np.random.random(nelts)*len(chain),int),:]
dchain = dict( zip(par,np.transpose(chain)))
dchain = dict( list(zip(par,np.transpose(chain))))
return(dchain)
......@@ -89,10 +89,10 @@ def MCparnames( filename, num=1, ext="txt"):
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"]
......@@ -124,7 +124,7 @@ def scan(filename):
index = index+1
print(l[index].replace("\n",""))
while not l[index].startswith("it "):
index = index+1
......@@ -141,38 +141,38 @@ def scan(filename):
#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))
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)
......@@ -197,12 +197,12 @@ def sig2ci( val, sig):
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)
......@@ -213,16 +213,16 @@ def cov2cor( mat):
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 = range(length_min,length_max,gap)[1:-1]
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)))
......@@ -239,21 +239,21 @@ def GelmanRubin( chains, gap=10000, length_min=1000, length_max=None, new=False)
#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
......@@ -261,22 +261,22 @@ def Geweke( chains, gap=10000, length_max=None):
if length_max > nsamples:
print( "Not enough samples...")
exit()
it = range(length_min,length_max,gap)[1:]
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)
########################################################################
......@@ -402,7 +402,7 @@ to_cosmomcName = {
'Aps143x217' : 'Aps143x217',
}
from_cosmomcName = dict(zip(to_cosmomcName.values(),to_cosmomcName.keys()))
from_cosmomcName = dict(list(zip(list(to_cosmomcName.values()),list(to_cosmomcName.keys()))))
......@@ -427,7 +427,7 @@ def posterior1d( chains, params, nbin=50, smooth=1,
if np.size(lw) == 1:
lw = [lw]*nchain
elif np.size(lw) < nchain:
lw = [1]*nchain
lw = [1]*nchain
if np.prod(subplot) < npar:
subplot[1] = np.ceil(np.double(npar)/subplot[0])
......@@ -442,16 +442,16 @@ def posterior1d( chains, params, nbin=50, smooth=1,
for c in range(nchain):
chain = chains[c]
if chain.has_key(par):
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:
......@@ -477,7 +477,7 @@ def posterior2d( chains, par1, par2, *args, **kwargs):
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)
......@@ -486,12 +486,12 @@ def posterior2d( chains, par1, par2, *args, **kwargs):
linewidths = kwargs.pop("linewidths", None)
fill = kwargs.get("fill", True)
datapoints = kwargs.get("datapoints", False)
contours = kwargs.get("contours", 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:
......@@ -508,22 +508,22 @@ def posterior2d( chains, par1, par2, *args, **kwargs):
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.)
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))
......@@ -539,12 +539,12 @@ def hist2d( chain, par1, par2, *args, **kwargs):
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:
......@@ -557,13 +557,13 @@ def hist2d( chain, par1, par2, *args, **kwargs):
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))
......@@ -572,7 +572,7 @@ def hist2d( chain, par1, par2, *args, **kwargs):
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)
......@@ -600,7 +600,7 @@ def triangle( chains, params, nbin=50, parnames=None, names=[], colors=[], lines
if len(linestyles) < nchain:
linestyles = ['-']*nchain
if weights == None:
weights = [np.ones(len(c['chi2'])) for c in chains]
......@@ -627,10 +627,10 @@ def triangle( chains, params, nbin=50, parnames=None, names=[], colors=[], lines
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))
......@@ -656,7 +656,7 @@ def triangle( chains, params, nbin=50, parnames=None, names=[], colors=[], lines
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])
......@@ -680,7 +680,7 @@ def triangle( chains, params, nbin=50, parnames=None, names=[], colors=[], lines
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)
......@@ -694,7 +694,7 @@ def triangle( chains, params, nbin=50, parnames=None, names=[], colors=[], lines
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
import matplotlib.ticker as mtick
#force tuple
if not( isinstance(chains, list) or isinstance(chains, tuple)):
......@@ -719,7 +719,7 @@ def rectangle( chains, par0, params, nbin=50, parnames=None, names=[], colors=[]
if len(linestyles) < nchain:
linestyles = ['-']*nchain
if weights == None:
weights = [np.ones(len(c['chi2'])) for c in chains]
......@@ -765,10 +765,10 @@ def rectangle( chains, par0, params, nbin=50, parnames=None, names=[], colors=[]
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):
......@@ -797,7 +797,7 @@ def MCcorrelation( chain, par, **kwargs):
data = [chain[p] for p in par]
matC = np.corrcoef( data)
plot_correlation( matC, par, **kwargs)
......@@ -810,8 +810,8 @@ def plot_correlation( matC, par, **kwargs):
covmat = np.tril( matC, k=-1)
covmat[covmat==0] = np.NaN
plt.imshow( covmat, vmin=vmin, vmax=vmax, **kwargs)
plt.yticks( range(len(covmat)), [parname.get(p,p) for p in par])
plt.xticks( range(len(covmat)), '')
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)
......@@ -854,10 +854,10 @@ def FeldmanCousins( xmin, sigma, CL=95):
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)
########################################################################
......@@ -876,11 +876,11 @@ def read_bestfit( dirname):
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:
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(zip(param, value))
bestfit = dict(list(zip(param, value)))
else :
bestfit = {'chi2':np.inf}
......@@ -889,15 +889,15 @@ def read_bestfit( dirname):
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]:
......@@ -916,17 +916,17 @@ def extract_bestfit( filename, skipone=False):
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:
......@@ -938,11 +938,11 @@ def extract_hess( filename, ext="LASymMatrix", reverse=False):
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
......@@ -950,7 +950,7 @@ def extract_hess( filename, ext="LASymMatrix", reverse=False):
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):
......@@ -958,7 +958,7 @@ def extract_hess( filename, ext="LASymMatrix", reverse=False):
if reverse:
matcov = flipud(matcov)
return(matcov,par)
......@@ -990,26 +990,26 @@ simple_lik_name = {
def writeExtendedBestFit( dbdir, bfwrite=False, outpath=""):
#Extract Bestfit from logfile
path, dirs, files = os.walk("%s" % (dbdir)).next()
path, dirs, files = next(os.walk("%s" % (dbdir)))
for idir in dirs:
p,d,f = os.walk("%s/%s" % (dbdir,idir)).next()
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)
......@@ -1027,7 +1027,7 @@ def writeExtendedBestFit( dbdir, bfwrite=False, outpath=""):
parlist = Parnames( camelparfile)
except ValueError:
return()
param_line = " ".join(parlist)
if CheckProfile:
......@@ -1045,14 +1045,14 @@ def writeExtendedBestFit( dbdir, bfwrite=False, outpath=""):
sf = open( name, 'r')
all_lines=sf.readlines()
# gets all lines in file
# 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
......@@ -1065,13 +1065,13 @@ def writeExtendedBestFit( dbdir, bfwrite=False, outpath=""):
#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
......@@ -1085,12 +1085,12 @@ def writeExtendedBestFit( dbdir, bfwrite=False, outpath=""):
if line_satisfied:
read=False
param_line = param_line.replace("\n"," ") + "\n"
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)
......@@ -1100,18 +1100,18 @@ def writeExtendedBestFit( dbdir, bfwrite=False, outpath=""):
sf = open( name, 'r')
# gets Profile line
goodline = False
values =""
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 += " "