Commit 68461d04 authored by Matthieu Tristram's avatar Matthieu Tristram
Browse files

Merge branch 'master' of https://gitlab.in2p3.fr/plaszczy/CAMEL

parents 1acee541 a8a7c1d9
......@@ -96,7 +96,7 @@ gauss1 SPT_high_220_cal 1.01 .02
#-----------------------------
#CMB Very-high-l (SPT-low)
SPT_Low = HighEll/SPT_low_plik.par
SPT_Low = HighEll/SPT_low.lik
par SPT_low_Aps nui 20.32 0.01 1. 60
par SPT_low_cal nui 1.00 0.01 0 2
#-----------------------------
......
#!/bin/bash
#EXPECTED_ARGS=5
#parFile=$1
#var=$2
#xmin=$3
#xmax=$4
#N=$5
#if [ $# -ne $EXPECTED_ARGS ]
#then
# echo "Usage: `basename $0` parfile var xmin xmax N"
# exit 1
#fi
#val array
##############define the scanned variable###############
# define the variable that is scanned
var="100*theta_s"
#values
# 1. define automatically: N points between xmin and xmax
N=50
xmin=1.04
xmax=1.0422
N=50
#compute automatically the values
for ((i=0;i<N;i++)) ; do
val[$i]=$(echo | awk -v xmin=$xmin -v xmax=$xmax -v N=$N -v i=$i 'BEGIN{step=(xmax-xmin)/(N-1)} {printf("%f",xmin+(i)*step)}')
done
#val=(0.9 1 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.4)
#val=(0.95 1.03 1.07 1.12 1.17)
#val=(0.95 1 1.03 1.05 1.07 1.10 1.12 1.15 1.20 1.25 1.30)
#2. or add explcitely the list of values
#var="Alens"
#val=(0.95 1 1.03 1.05 1.07 1.10 1.12 1.15 1.20 1.25 1.30)
#var=tau_reio
#val=(0.05 0.06 0.07 0.08 0.09 0.10 0.11 0.12 0.13 0.14 0.15 0.15 0.16 0.17 0.18 0.19 0.20 0.21 0.22 0.23 0.24 0.25)
#val=(0.025 0.035 0.045 0.055 0.065 0.075 0.085 0.095)
#val=(0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 0.055 0.06 0.065 0.07 0.075 0.08 0.085 0.09 0.095 0.10 0.11 0.12)
#val=(0.16 0.17 0.18 0.10 0.20 0.21 0.22 0.23 0.24 0.25)
echo "About to run a scan on the ---$var--- variable with the following values"
echo ${val[*]}
echo "Is this OK? [y/n]"
read answer
if [ $answer != 'y' ] ; then
echo "exiting"
exit 1
fi
#var=N_eff
#TE
#val=(1.5 1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3. 3.2 3.5)
#val=(0.7 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3. 3.2 3.5 )
#EE
#val=(0.1 0.5 1. 1.5 2 2.5 3 3.5 4 4.5 5)
#;val=(0.2 0.3 0.4 0.6 0.7 0.8 0.9 1.1 1.2 1.3 1.4 1.6 1.7 1.8 1.9 2.1 2.2 2.3 2.4 2.6 2.7 2.8 2.9 3.1 3.2 3.3 3.4 3.6 3.7)
#TT
#val=(0.5 0.7 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3. 3.2 3.4 )
#val=(3.5 3.6 3.7 3.8 3.9 4)
######################################
#escape problematic characters
varesc=$(printf "%q" $var)
cd `dirname $1`
dirpar=$PWD
cd -
......@@ -55,7 +43,7 @@ parbase=`basename $1`
file=$dirpar/$parbase
dirout=${parbase%".par"}_mprof_"${var}"
dirout=${parbase%".par"}_mprof_"$varesc"
if [ -d "$dirout" ] ; then
echo "directory $dirout exists: still want to run? (o/n)"
......@@ -92,7 +80,7 @@ for i in $(seq 0 $((${#val[@]}-1))) ; do
zeval=${val[$i]}
echo $zeval
OUTDIR=$HERE/$dirout/"$var"_${zeval}
OUTDIR=$HERE/$dirout/"$varesc"_${zeval}
mkdir $OUTDIR
cd $OUTDIR
###################################################
......@@ -119,9 +107,9 @@ cp $CAMELROOT/batch/awk/genrand.awk .
#
#input file -> transforme en pico et enleve var
grep -v "$var" $file > parfile_in
grep -v "$varesc" $file > parfile_in
#compute var
#compute var: here no escape
echo "fix $var cosmo $zeval" >> parfile_in
echo "init par:"
......
""" Uitlity to extact the best-fit in each sub-directory given in argument
and put it in a single file """
import numpy as np
import sys,glob,os,string
if len(sys.argv)==1:
print("Usage: bestbestfit.py dir-list (you may use *)")
sys.exit()
header=False
for dir in sys.argv[1:]:
bestfits=glob.glob(dir+"/best_fit*")
chi2min=sys.float_info.max
data=np.zeros(1)
for bf in bestfits:
if not header:
with open(bf,'r') as f:
line = f.readline().rstrip()
labels=line.split()
ichi2=labels.index('chi2')
varname=labels[ichi2-1]
#open output file
fname="profile_"+varname+".txt"
if os.path.isfile(fname):
os.remove(fname)
fout=open(fname,"a")
fout.write(line)
fout.write('\n')
header=True
d=np.loadtxt(bf,skiprows=1)
if d[ichi2]<chi2min:
data=d
chi2min=d[ichi2]
fout.write(string.join(map(str,data),"\t"))
fout.write('\n')
fout.close()
print("output written to file: "+fname)
import numpy as np
import sys
import matplotlib.pyplot as plt
file=sys.argv[1]
#read labels
with open(file,'r') as f:
label = f.readline().rstrip()
l=label.split()
d=np.loadtxt(file,skiprows=1)
ichi2=l.index('chi2')
ivar=ichi2-1 #
##determine how many variables are not fixed
fixed=[sum(d[:,i]-d[0,i])==0 for i in range(ivar-1)]
nsub=sum(np.invert(fixed))
ncol=int(np.sqrt(nsub))
nrow=int(nsub/ncol)
if (nsub%nrow != 0):
nrow=nrow+1
##plot
fig=plt.figure(figsize=(20,15))
iplot=0
for i in range(ivar-1):
if not fixed[i]:
iplot=iplot+1
ax=fig.add_subplot(nrow,ncol,iplot)
ax.plot(d[:,ivar],d[:,i])
ax.set_xlabel(l[ivar])
ax.set_ylabel(l[i])
ax.tick_params(labelsize=8)
plt.tight_layout()
plt.show()
\ No newline at end of file
""" Utility to draw a profile leikelihood from a bestfit file using matplotlib"""
from __future__ import print_function
import numpy as np
import sys
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
file=sys.argv[1]
#read labels
with open(file,'r') as f:
label = f.readline().rstrip()
l=label.split()
d=np.loadtxt(file,skiprows=1)
ichi2=l.index('chi2')
ivar=ichi2-1 # the scanned variable index is always below the chi2 value
chi2=d[:,ichi2]
var=d[:,ivar]
f=interp1d(var,chi2,kind='quadratic')
x=np.linspace(var[0],var[-1],100)
y=f(x)
ymin=y.min()
xmin=x[y==y.min()][0]
w=np.where(y-ymin<1)
low=w[0][0]
high=w[0][-1]
xlow=(x[low]+x[low-1])/2
xhigh=(x[high]+x[high+1])/2
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(var,chi2-ymin,'ok')
ax.plot(x,y-ymin,'r')
ax.plot([var[0],var[-1]],[1,1],'b--')
ax.plot([xlow,xlow],[-0.1,1],'b')
ax.plot([xhigh,xhigh],[-0.1,1],'b')
#ax.annotate("",xy=(xlow,-0.1),xytext=(xlow,1),arrowprops=dict(arrowstyle='->'))
#ax.annotate("",xy=(xhigh,-0.1),xytext=(xhigh,1),arrowprops=dict(arrowstyle='->'))
ax.set_xlabel(l[ivar])
ax.set_ylabel("chi2-min(chi2)")
ax.set_ylim((-0.1,5))
result="{:f}-{:f}+{:f}".format(xmin,xmin-xlow,xhigh-xmin)
print("{}={}".format(l[ivar],result))
ax.text(xmin,4,result,horizontalalignment='center')
plt.show()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment