Commit 73f5831b authored by Plaszczynski Stephane's avatar Plaszczynski Stephane
Browse files

first import

parent 7fcb7673
BEGIN{
nerr=0
nmaxit=0
fn=FILENAME
}
/function Value/ {
chi2min=$4
}
/calls/{
niter=$NF
}
/TIMER/ {
tt=$0
split($3,a,"=")
tsec=a[2]
}
/Error running/{
nerr++
}
/MAX_IT/{
nmaxit++
}
# changement de fichier: dump
FNR==0{
print FILENAME
fn=FILENAME
}
END{
if (niter>0)
print file,chi2min,niter,tsec/niter,nerr,nmaxit
}
#!/bin/env awk
BEGIN {
t=0
N=0
it=0
}
$1==it && NF>10{
#print $0
t+=$NF
N++
it++
}
/Gauss/{
print $0
}
/WARNING/ {
print $0
}
/function Value/ {
chi2min=$4
print "chi2min=",chi2min
}
/calls/{
print $0
niter=$NF
print "niter=",niter
}
/Chi2 from/{
chi2=$4
ndof=$6
}
/TIMER/ {
tt=$0
print tt
split($3,a,"=")
tsec=a[2]
}
/===/{
print $0
}
/converge/ { print $0 }
/Name/ {cosmo=0}
$2 ~ /\|\|/{
split($0,a,"|")
num=a[1]
varname=a[3]
val=a[7]
err=a[9]
type=a[5]
if (cosmo < 7)
printf("||%d|| %s\t|| %10.5g\t|| %10.5g\t||%s\t||\n",num,varname,val,err,type)
# printf("||%d|| %s\t|| %10.5g+-%10.5g\t||\n",num,varname,val,err)
else
printf("||%d|| %s\t|| %5.3g\t|| %5.3g\t||%s\t||\n",num,varname,val,err,type)
# printf("||%d|| %s\t|| %5.3g+- %5.3g\t||\n",num,varname,val,err)
cosmo++
}
END{
print "chi2min=",chi2min
print "#iter=",niter,"\ttot time=",tsec/3600.," h\ttime/iter=",tsec/niter,"s"
# print "Hillipop: chi2=",chi2," ndof=",ndof,"(",(chi2-ndof)/sqrt(2*ndof), "sigma)"
}
BEGIN {
srand(seed ? seed : 1)
}
$1 =="par"{
r=2*rand()-1
newval=$4+r*$5
printf("%s\t%s\t%s\t%f\t%f\t%f\t%f\n",$1,$2,$3,newval,$5,$6,$7)
}
$1 !="par"{
print $0
}
NR == 40000{
N=0
for (i=1;i<=NF;i++) {
s[i]=0
s2[i]=0
}
}
NR > 40000 {
N++
for (i=1;i<=NF;i++) {
s[i]+=$i
s2[i]+=($i)**2
}
}
END{
for (i=1;i<=NF;i++) {
mu=s[i]/N
a=sqrt(s2[i]/N)
var=(a-mu)*(a+mu)
printf("%3.12f\t%3.12f\t",mu,sqrt(var))
}
printf("\n")
}
BEGIN{
print "engine=class"
}
!/engine/ && !/precisionFile/ && $2 != "m_ncdm" && $2 != "N_nu" && $2 != "YHe"{
print $0
}
END{
print "class N_ncdm 1"
print "class m_ncdm 0.06"
print "class N_eff 2.046"
print "class k_pivot 0.05"
print "class lensing yes"
print "class sBBN\\ file bbn/sBBN.dat"
print "precisionFile=hpjul2.pre"
}
BEGIN{
print "engine=pico"
Alens=0
N_nu=0
M_ncdm=0
}
#############
NF==3 && $2 == "Alens"{
Alens=$3
print "fix Alens cosmo ",Alens
}
##############
NF==3 && $2 == "m_ncdm"{
m_ncdm=$3
print "fix m_ncdm cosmo ",m_ncdm
}
#############
$2 == "N_eff"{
if (NF==3) {
N_nu=$3+1
print "fix N_nu cosmo ",N_nu
}
else{
N_nu=$4+1
print $1," N_nu cosmo ",N_nu,$5,$6,$7
}
}
#############
!/class/ && $2 != "N_eff"{
print $0
}
############
END{
#fill if missing
if (m_ncdm==0)
print "fix m_ncdm cosm 0.06 "
if (Alens==0)
print "fix Alens cosm 1"
if (N_nu==0)
print "fix N_nu cosm 3.046"
print "fix YHe cosm 0.248"
}
n_chain=4L
t1=mrdfits(dir+"/samples.fits",1,hdr)
name=strtrim(hdr[5:n_elements(hdr)-4],2)
t2=mrdfits(dir+"/samples.fits",2)
t3=mrdfits(dir+"/samples.fits",3)
t4=mrdfits(dir+"/samples.fits",4)
sz1=size(t1,/dim)
sz2=size(t2,/dim)
sz3=size(t3,/dim)
sz4=size(t4,/dim)
dim=sz1[1]
stop=min([sz1[0]-1,sz2[0]-1,sz3[0]-1,sz4[0]-1])
rec=dblarr(stop+1, dim, n_chain)
flength=dim*stop
print, flength
rec[*,*,0]=t1[0:stop,*]
rec[*,*,1]=t2[0:stop,*]
rec[*,*,2]=t3[0:stop,*]
rec[*,*,3]=t4[0:stop,*]
lenght_min=long(50)
lenght_max=stop
gap=1000L
event=(lenght_max-lenght_min)/gap+1
R=dblarr(dim,event) ;vectore GR
;vettore di medie delle singole catene
;vector of means of single chain
mean_chain=fltarr(dim,n_chain)
B_n=fltarr(dim) ; variance between chains
W=fltarr(dim) ; variance of the chain
W_c=fltarr(dim,n_chain) ; W for fixed n and fixed burn_in=stop/2
mean_dist=fltarr(dim) ; for all dim mean of the chain
index_event=0
for index=lenght_min, lenght_max-1 do begin &$
for d=1,dim-1 do begin &$
if stop ne 0 then begin &$
for m=0,n_chain-1 do begin &$
burn_in=index/2 &$
mean_chain[d,m]=total(rec[burn_in :index,d,m])/(index-burn_in) &$
W_c[d,m]=total((rec[burn_in :index,d,m]-mean_chain[d,m])^2)/(index-burn_in) &$
endfor &$
endif &$
mean_dist[d]=total(mean_chain[d,*])/n_chain &$ ;mean of the distribution
if n_chain gt 1 then begin &$
B_n[d]=total((mean_chain[d,*]-mean_dist[d])^2)/(n_chain-1) &$
endif &$
W[d]=total(W_c[d,*])/n_chain &$
if W[d] ne 0 then R[d,index_event]=(float((index-1.)/index)*W[d] + B_n[d]*(1.+1./n_chain))/W[d] &$
endfor &$
index_event=index_event+1 &$
index=index+gap-1 &$
endfor
;print, R
;-------------plot R
if (n_elements(write) eq 0) then write=0
if (write ne 0) then openplot, "GR_"+dir+'.ps', xsize=20, ysize=12
loadct, 39
x=indgen(event)*gap+lenght_min
!x.style=1
y=fltarr(event)
y[*]=0.03
plot, x ,y, title='Gelman Rubin R vs length of the chain', xtitle='length of the chain', ytitle='R-1', yrange=[1e-3,1],ylog=ylog
linestyle=[0,0,0,0,0,0,2]
color=[30,90,150,110,200,220,70]
vec_one=fltarr(event)
vec_one=replicate(1,event)
for i=1,6 do oplot, x, R[i,*]-vec_one, color=color[i-1]
for i=7,dim-1 do oplot, x, R[i,*]-vec_one, color=color[6], linestyle=linestyle[6]
;items=[ 'omega_b' , 'omega_cdm','100*theta_s', 'tau_reio', 'n_s', 'log(10^10A_s)','nuisances']
items=name[1:6]
items=[items,"nuisances"]
legend, items,colors=color,linestyle=linestyle, /top, /right ; textcolors=color
;add more nuisances here
if (write ne 0) then closeplot, /png
;Alens
sig2=0.07^2
ipos=6
cov=readmatrix("covmat.dat")
sz=size(cov,/dim)
cov_new=make_array(sz+1,/double,val=0)
cov_new[0:sz[0]-1,0:sz[1]-1]=cov
cov_new[sz[0],sz[1]]=sig2
w=lindgen(sz[0])
w_new=[w[0:ipos-1],sz[0]+1,w[ipos:*]]
c1=cov_new[*,w_new]
c2=c1[w_new,*]
cov_new=c2
name_new=[name,"Alens"]
cor_new=cov2cor(cov_new)
plotcor,cor_new,name_new[w_new]
writematrix,cov_new,f="covmat_new.dat"
pro autozone,n
ncol=long(sqrt(n))
nrow=N/ncol
if (N mod nrow ne 0) then nrow++
!p.multi=[0,ncol,nrow]
return
end
file=dir+"/bestfit.det"
if (not file_test(file)) then stop,'error file'+file
t=read_ascii(file)
t=transpose(t.(0))
readcol,file,files,f='A'
n=n_elements(files)
sz=size(t,/dim)
n=sz[0]
;look for names
camels=file_search(dir+"/*/camel*.par")
names=parnames(camels[0],fix=fix)
namevar=fix[n_elements(fix)-1]
print,"profile on: ",namevar
print,+"# free params="+strn(n_elements(names))
print,names
; prof2
;ichi=sz[1]-12
;Mnu
ichi=sz[1]-11
;
imu=ichi-1
;imu=ichi-2 ; Alens class
m=t[*,imu]
chi2=t[*,ichi]
dchi2=chi2-min(chi2)
zone
window,/free,tit=file
plot,m,chi2-min(chi2),/ys,/psym,xtit=namevar,ytit=textoidl('\Delta \chi^2'),tit=file,yr=[0,4]
;best fit uniq
u=uniq(m)
mu=m[u]
mu=mu[sort(mu)]
Nu=n_elements(u)
bf=make_array(Nu,sz[1],/double)
for iu=0,Nu-1 do begin &$
w=where(m eq mu[iu]) &$
sub=t[w,*] &$
c=min(sub[*,ichi],imin) &$
print,mu[iu],double(c) &$
bf[iu,*]=transpose(sub[imin,*]) &$
endfor
oplot,mu,bf[*,ichi]-min(bf[*,ichi],imin),col=red
;imin=2
!x.margin=[15,5]
;window,/free,xs=1200,ys=700,tit=file
;cgdisplay,/free,xs=1200,ys=700,tit=file
cgdisplay,900,700,/free,tit=file
autozone,n_elements(names)
ii=0
for i=0,sz[1]-1 do begin &$
if (stddev(bf[*,i]) eq 0) then continue &$
if (ii ge n_elements(names)) then break &$
cgplot,mu,bf[*,i],/ys,xtit=namevar,ytit=names[ii],psym=-16,symsize=0.1 &$
print,names[ii],bf[imin,i] &$
ii++ &$
endfor
print,namevar,mu[imin]
;fichier avec mnu rajoute chi2 rajoute au debut
print,"writing",dir+"/bestfit.all"
openw,lun,dir+"/bestfit.all",/get_lun
for i=0,n_elements(mu)-1 do begin &$
line=[bf[i,imu],bf[i,ichi],transpose(bf[i,*])] &$
printf,lun,strjoin(strtrim(string(line),2)," ") &$
endfor
free_lun,lun
print,bf[imin,*]
!x.margin=[6,5]
xfit=0.2
yup=18
xr=[-.5,.5]
yr=[0,yup]
;dirs=["hlpTT_bflike_vhl_LCDMClass_nuscan","hlpTT_vhl_lensing_LCDMClass_nuscan"]
;dirs=["hlpTT_bflike_vhl_LCDMClass_nuscan","hlpTT_bflike_vhl_LCDMClass_nuscan0"]
dirs=["hlpTT_bflike_vhl_LCDMClass_nuscan","hlpTT_bflike_vhl_JLA_LCDMClass_nuscan"]
names=dirs
cols=['black','red','blue','green','orange','magenta','violet','purple','navy']
;paires
;cols=['red','red5','blue','navy','green','forest green','violet','violet red']
varname=textoidl('\Sigma m_\nu [eV]')
cgplot,findgen(10),xr=xr,yr=yr,xtit=varname,ytit=textoidl('\Delta \chi^2'),ys=1,/xs,/nodata
xp=[xr[0],0,0,xr[0]]
yp=[yr[0],yr[0],yup,yup]
polyfill,xp,yp,/data,/line_fill,ORIENTATION=45, THICK=1.5, SPACING=0.1,col=grey
in=''
N=n_elements(dirs)
for i=0,N-1 do begin &$
if (i eq 0) then over=0 else over=1 &$
t=read_ascii(dirs[i]+'/bestfit.all') &$
t=transpose(t.(0)) &$
sz=size(t,/dim) &$
ichi=1 &$
itau=0 &$
chi2=t[*,ichi] &$
dchi2=chi2-min(chi2,imin) &$
tau=t[*,itau] &$
iord=sort(tau) &$
tau=tau[iord] &$
dchi2=dchi2[iord] &$
diff=dchi2-shift(dchi2,-1) &$
w=where(dchi2 lt 5 and abs(diff) lt 1) &$ ; points used
;w=where(dchi2 gt 0) &$
x=tau &$
chi2min1=min(chi2) &$
chi2-=chi2min1 &$
cut=where(x lt xfit) &$
par=svdfit(x[cut],chi2[cut],3,/double) &$
xmle=-par[1]/(2*par[2]) &$
chi2min=poly(xmle,par) &$
sig=1/sqrt(par[2]) &$
print,dirs[i]+' xmle=',+strn(xmle),' sigma='+strn(sig)+' nsig='+strn(xmle/sig) &$
cgplot,x,chi2-chi2min,psym=16,col=cols[i],/over &$
xx=range(xr[0],xr[1],N=100) &$
cgplot,xx,poly(xx,par)-chi2min,linestyle=2,col=cols[i],/over &$
xx=range(0,xfit,N=100) &$
cgplot,xx,poly(xx,par)-chi2min,col=cols[i],/over &$
;feldman
xf=[-.5,-.4,-.3,-.2,-.1,0.,.1,.2,.3,.4,.5] &$
yf95=[1.49,1.58,1.67,1.77,1.86,1.96,2.06,2.16,2.26,2.36,2.46] &$
mu2=interpol(yf95,xf,xmle/sig) &$
x0=mu2*sig &$
print, 'FC=',x0 &$
y0=poly(x0,par)-chi2min &$
cgplot,[0,x0],[y0,y0],col=cols[i],thick=2,/over &$
cgplot,[x0,x0],[0,y0],col=cols[i],thick=2,/over &$
endfor
;hlp
;plik
;cgtext,0.6,0.8,names[0],/normal,col=cols[0]
;cgtext,0.6,0.75,names[1],/normal,col=cols[1]
al_legend,names,textcol=cgcolor(cols),/top,/right,box=0
.r parnames
if (n_elements(write) eq 0) then write=0
fileout=dir+'.ps'
if (n_elements(t) eq 0) then begin &$
file=dir+'/samples_cut.fits' &$
print,"reading : ",file &$
t=mrdfits(file,1,hdr)&$
sz=size(t,/dim) &$
n=sz[0] &$
;name=parnames(dir+'/camel.par') &$
;name=['chi2',name] &$
name=strtrim(hdr[5:n_elements(hdr)-4],2) &$
derived,t,name &$
endif
;noms/val
;if (write eq 0) then
if (write ne 0) then openplot,fileout,xs=30,ys=30,charsize=2 else cgdisplay,900,700,tit=dir,/free
N=n_elements(name)
autozone,N
!p.charsize=1.2
delvar,name2
if (file_test(dir+"/bf")) then begin &$
print, "best fit file found",dir+"/bf" &$
readcol,dir+'/bf',delim='||',f='x,a,f',name2,bf &$
name2=strtrim(name2,2) &$
endif
line=['']
!x.margin=[5,5]
for i=1,N-1 do begin &$
if (total(t[*,i]-t[0,i]) eq 0) then continue &$
m=mean(t[*,i]) &$
s=stddev(t[*,i]) &$
hist_plot,t[*,i],/noplot,x=x,y=y,nbin=101 &$
;hist_plot,t[*,i],/noplot,x=x,y=y,nbin=50 &$
yy=smooth(y,15,/edge_trunc)&$
yy/=max(yy) &$
cgplot,x,yy,xtit=name[i],yr=[0,1],/xs,xticks=3,yticks=3,ytickv=[0.,0.5,1] &$
ibf=-1 &$
if (n_elements(name2) ne 0) then begin &$
ibf=strindex(name2,name[i],/silent) &$
cgplot,[bf[ibf],bf[ibf]],!y.crange,/over,line=2 &$
endif &$
if (i le 6 ) then dec=5 else dec=3 &$
ss=cgNumber_Formatter(m,dec=dec)+'+-'+cgNumber_Formatter(s,dec=dec) &$
cgtext,m-2*s,0.6,ss,/data,charsize=1,col='blue' &$
if (ibf ge 0) then begin&$
print,'||',name[i],'||',cgNumber_Formatter(bf[ibf],dec=3),'||',ss,'||' &$
endif else begin &$
print,'||',name[i],'|| ||',ss,'||' &$
endelse &$
line=[line,strn(m)] &$
endfor
if (write ne 0) then closeplot,/pdf,/view
;if (write ne 0) then closeplot,/png
.r strindex
.r derived
;if (n_elements(mcr) eq 0) then mcr=[80000L,150000L]
;print,'range=',mcr
!X.OMargin = [5, 2]
!Y.OMargin = [5, 5]
!y.margin=[3,2]
!x.margin=[5,2]
if (n_elements(write) eq 0) then write=0
fileout=dir+'_vs_'+dir2+'.ps'
delvar,t
file=dir+'/samples_cut.fits' &$
print,"reading : ",file &$
t=mrdfits(file,1,hdr)&$
name=strtrim(hdr[5:n_elements(hdr)-4],2) &$
derived,t,name
sz=size(t,/dim) &$
;le plus grand
delvar,t2
file=dir2+'/samples_cut.fits' &$
print,"reading : ",file&$
t2=mrdfits(file,1,hdr) &$
name2=strtrim(hdr[5:n_elements(hdr)-4],2) &$
derived,t2,name2
if (n_elements(write) eq 0) then write=0
if (write ne 0) then openplot,fileout,xs=30,ys=30,charsize=1.5 else cgdisplay,900,700,tit=dir,/free
if (write eq 0) then !p.thick=3 else !p.thick=8
if (write eq 0) then !p.charsize=1
N=n_elements(name)
autozone,N
!x.margin=[5,5]
for i=1,N-1 do begin &$
if (total(t[*,i]-t[0,i]) eq 0) then continue &$
;dans liste2
j=strindex(name2,name[i],/silent) &$
print,name[i],j &$
if (j lt 0) then continue &$
print,i,j,name[i],name2[j] &$
hist_plot,t[*,i],/noplot,x=x,y=y,nbin=101 &$
yy=smooth(y,15,/edge_trunc)&$
yy/=max(yy) &$
hist_plot,t2[*,j],/noplot,x=x2,y=y2,nbin=101 &$
yy2=smooth(y2,15,/edge_trunc)&$
yy2/=max(yy2) &$
xx=[x,x2] &$
cgplot,x,yy,xtit=name[i],xr=minmax(xx),color='dark green',yticks=4,yr=[0,1] &$