Skip to content
Snippets Groups Projects
Commit e4815c0e authored by frcojimenez's avatar frcojimenez
Browse files

added config -c option and 1-sigma band plots

parent 60415167
Branches
No related tags found
No related merge requests found
This diff is collapsed.
#!/usr/bin/env python #!/usr/bin/env python
# coding: utf-8 # coding: utf-8
# In[99]: # In[63]:
#Import relevant modules, import data and all that #Import relevant modules, import data and all that
...@@ -30,6 +30,10 @@ import dynesty ...@@ -30,6 +30,10 @@ import dynesty
from dynesty import plotting as dyplot from dynesty import plotting as dyplot
import os import os
import csv import csv
import argparse
import scipy.optimize as optimization
from scipy.optimize import minimize
#Remember to change the following global variables #Remember to change the following global variables
#rootpath: root path to nr data #rootpath: root path to nr data
...@@ -38,16 +42,25 @@ import csv ...@@ -38,16 +42,25 @@ import csv
#tshift: time shift after the strain peak #tshift: time shift after the strain peak
#vary_fund: whether you vary the fundamental frequency. Works in the model_dv function. #vary_fund: whether you vary the fundamental frequency. Works in the model_dv function.
try:
parser = argparse.ArgumentParser(description="Simple argument parser")
parser.add_argument("-c", action="store", dest="config_file")
result = parser.parse_args()
config_file=result.config_file
parser = ConfigParser()
parser.read(config_file)
parser.sections()
except SystemExit:
parser = ConfigParser() parser = ConfigParser()
parser.read('config_q10n.ini') parser.read('config.ini')
parser.sections() parser.sections()
pass
# In[100]: # In[36]:
# path # path
rootpath=parser.get('nr-paths','rootpath') rootpath=parser.get('nr-paths','rootpath')
simulation_path_1 = parser.get('nr-paths','simulation_path_1') simulation_path_1 = parser.get('nr-paths','simulation_path_1')
...@@ -57,17 +70,18 @@ simulation_number = parser.get('nr-paths','simulation_number') ...@@ -57,17 +70,18 @@ simulation_number = parser.get('nr-paths','simulation_number')
simulation_number = np.int(simulation_number) simulation_number = np.int(simulation_number)
output_folder = parser.get('output-folder','output-folder') output_folder = parser.get('output-folder','output-folder')
overwrite = parser.get('overwrite','overwrite')
# In[101]: # In[37]:
if not os.path.exists(output_folder): if not os.path.exists(output_folder):
os.mkdir(dirName) os.mkdir(output_folder)
print("Directory " , output_folder , " Created ") print("Directory " , output_folder , " Created ")
# In[102]: # In[38]:
# time config # time config
...@@ -82,7 +96,7 @@ t_align=parser.get('time-setup','t_align') ...@@ -82,7 +96,7 @@ t_align=parser.get('time-setup','t_align')
t_align = np.float(t_align) t_align = np.float(t_align)
# In[103]: # In[39]:
# n-tones & nlive # n-tones & nlive
...@@ -94,7 +108,7 @@ npoints=parser.get('n-live-points','npoints') ...@@ -94,7 +108,7 @@ npoints=parser.get('n-live-points','npoints')
npoints = np.int(npoints) npoints = np.int(npoints)
# In[104]: # In[40]:
# model # model
...@@ -106,8 +120,26 @@ elif model == 'w-q': ...@@ -106,8 +120,26 @@ elif model == 'w-q':
elif model == 'w-tau-fixed': elif model == 'w-tau-fixed':
tau_var_str='q' tau_var_str='q'
print('model:',model)
print('nmax',nmax)
# In[41]:
tshift
# In[105]: # In[42]:
output_folder_1=output_folder+'/'+model+'-nmax'+str(nmax)
if not os.path.exists(output_folder_1):
os.mkdir(output_folder_1)
print("Directory " , output_folder_1 , " Created ")
# In[43]:
# loading priors # loading priors
...@@ -156,7 +188,7 @@ if model == 'w-tau-fixed': ...@@ -156,7 +188,7 @@ if model == 'w-tau-fixed':
prior_dim = len(priors_min) prior_dim = len(priors_min)
# In[106]: # In[44]:
vary_fund = True vary_fund = True
...@@ -195,8 +227,16 @@ def EasyMatchT(t,h1,h2,tmin,tmax): ...@@ -195,8 +227,16 @@ def EasyMatchT(t,h1,h2,tmin,tmax):
return res return res
def wRD_to_f_Phys(f,M):
c=2.99792458*10**8;G=6.67259*10**(-11);MS=1.9885*10**30;
return (c**3/(M*MS*G*2*np.pi))*f
# In[107]: def tauRD_to_t_Phys(tau,M):
c=2.99792458*10**8;G=6.67259*10**(-11);MS=1.9885*10**30;
return ((M*MS*G)/c**3)*tau
# In[45]:
#This loads the 22 mode data #This loads the 22 mode data
...@@ -228,7 +268,7 @@ tmax5=FindTmaximum(gw5_sxs_bbh_0305) ...@@ -228,7 +268,7 @@ tmax5=FindTmaximum(gw5_sxs_bbh_0305)
times5 = times5 - tmax5 times5 = times5 - tmax5
# In[108]: # In[46]:
#Select the data from 0 onwards #Select the data from 0 onwards
...@@ -240,7 +280,7 @@ timesrd=gw_sxs_bbh_0305[position:-1][:,0][:-1]-tmax ...@@ -240,7 +280,7 @@ timesrd=gw_sxs_bbh_0305[position:-1][:,0][:-1]-tmax
timesrd5=gw5_sxs_bbh_0305[position5:-1][:,0][:-1]-tmax5 timesrd5=gw5_sxs_bbh_0305[position5:-1][:,0][:-1]-tmax5
# In[109]: # In[47]:
#Test plot real part (data was picked in the last cell). Aligning in time #Test plot real part (data was picked in the last cell). Aligning in time
...@@ -252,7 +292,7 @@ plt.plot(timesrd5, np.sqrt(gw_sxs_bbh_0305rd5[:,1]**2+gw_sxs_bbh_0305rd5[:,2]**2 ...@@ -252,7 +292,7 @@ plt.plot(timesrd5, np.sqrt(gw_sxs_bbh_0305rd5[:,1]**2+gw_sxs_bbh_0305rd5[:,2]**2
plt.legend() plt.legend()
# In[110]: # In[48]:
#Test plot im part (data was picked in the last cell). Aligning in time #Test plot im part (data was picked in the last cell). Aligning in time
...@@ -264,7 +304,7 @@ plt.plot(timesrd5, np.sqrt(gw_sxs_bbh_0305rd5[:,1]**2+gw_sxs_bbh_0305rd5[:,2]**2 ...@@ -264,7 +304,7 @@ plt.plot(timesrd5, np.sqrt(gw_sxs_bbh_0305rd5[:,1]**2+gw_sxs_bbh_0305rd5[:,2]**2
plt.legend() plt.legend()
# In[111]: # In[49]:
# Depending on nmax, you load nmax number of freqs. and damping times from the qnm package # Depending on nmax, you load nmax number of freqs. and damping times from the qnm package
...@@ -273,7 +313,7 @@ w = (np.real(omegas))/mf ...@@ -273,7 +313,7 @@ w = (np.real(omegas))/mf
tau=-1/(np.imag(omegas))*mf tau=-1/(np.imag(omegas))*mf
# In[112]: # In[50]:
gwnew_re = interpolate.interp1d(timesrd, gw_sxs_bbh_0305rd[:,1], kind = 'cubic') gwnew_re = interpolate.interp1d(timesrd, gw_sxs_bbh_0305rd[:,1], kind = 'cubic')
...@@ -283,7 +323,7 @@ gwnew_re5 = interpolate.interp1d(timesrd5, gw_sxs_bbh_0305rd5[:,1], kind = 'cubi ...@@ -283,7 +323,7 @@ gwnew_re5 = interpolate.interp1d(timesrd5, gw_sxs_bbh_0305rd5[:,1], kind = 'cubi
gwnew_im5 = interpolate.interp1d(timesrd5, gw_sxs_bbh_0305rd5[:,2], kind = 'cubic') gwnew_im5 = interpolate.interp1d(timesrd5, gw_sxs_bbh_0305rd5[:,2], kind = 'cubic')
# In[113]: # In[51]:
if timesrd5[-1]>= timesrd[-1]: if timesrd5[-1]>= timesrd[-1]:
...@@ -300,7 +340,7 @@ gwdatanew = gwdatanew_re - 1j*gwdatanew_im ...@@ -300,7 +340,7 @@ gwdatanew = gwdatanew_re - 1j*gwdatanew_im
gwdatanew5 = gwdatanew_re5- 1j*gwdatanew_im5 gwdatanew5 = gwdatanew_re5- 1j*gwdatanew_im5
# In[114]: # In[52]:
mismatch=1-EasyMatchT(timesrd_final,gwdatanew,gwdatanew5,0,0+90) mismatch=1-EasyMatchT(timesrd_final,gwdatanew,gwdatanew5,0,0+90)
...@@ -308,7 +348,7 @@ error=np.sqrt(2*mismatch) ...@@ -308,7 +348,7 @@ error=np.sqrt(2*mismatch)
print(mismatch) print(mismatch)
# In[115]: # In[53]:
# Phase alignement # Phase alignement
...@@ -320,7 +360,7 @@ plt.plot(timesrd_final, phas, "r", alpha=0.3, lw=3, label=r'$phase$') ...@@ -320,7 +360,7 @@ plt.plot(timesrd_final, phas, "r", alpha=0.3, lw=3, label=r'$phase$')
plt.plot(timesrd_final, phas5, "blue", alpha=0.3, lw=3, label=r'$phase$') plt.plot(timesrd_final, phas5, "blue", alpha=0.3, lw=3, label=r'$phase$')
# In[116]: # In[54]:
position = np.argmax(timesrd_final >= (t_align)) position = np.argmax(timesrd_final >= (t_align))
...@@ -339,15 +379,15 @@ plt.plot(timesrd_final, phas, "r", alpha=0.3, lw=3, label=r'$phase$') ...@@ -339,15 +379,15 @@ plt.plot(timesrd_final, phas, "r", alpha=0.3, lw=3, label=r'$phase$')
plt.plot(timesrd_final, phas5, "blue", alpha=0.3, lw=3, label=r'$phase$') plt.plot(timesrd_final, phas5, "blue", alpha=0.3, lw=3, label=r'$phase$')
# In[117]: # In[55]:
mismatch=1-EasyMatchT(timesrd_final,gwdatanew,gwdatanew5,0,+90) mismatch=1-EasyMatchT(timesrd_final,gwdatanew,gwdatanew5,0,+90)
print(mismatch) print(mismatch)
error = np.sqrt(2*gwdatanew*gwdatanew-2*gwdatanew*gwdatanew5) error = np.sqrt(gwdatanew*gwdatanew-2*gwdatanew*gwdatanew5+gwdatanew5*gwdatanew5)
# In[118]: # In[56]:
#Test the new interpolated data #Test the new interpolated data
...@@ -358,7 +398,7 @@ plt.plot(timesrd_final, error.real, "b", alpha=0.3, lw=2, label='error') ...@@ -358,7 +398,7 @@ plt.plot(timesrd_final, error.real, "b", alpha=0.3, lw=2, label='error')
plt.legend() plt.legend()
# In[119]: # In[57]:
#Test the error data #Test the error data
...@@ -369,7 +409,7 @@ plt.plot(timesrd_final, np.sqrt(error.imag**2+error.real**2), "r", alpha=0.3, lw ...@@ -369,7 +409,7 @@ plt.plot(timesrd_final, np.sqrt(error.imag**2+error.real**2), "r", alpha=0.3, lw
plt.legend() plt.legend()
# In[120]: # In[58]:
#Take the piece of waveform you want #Take the piece of waveform you want
...@@ -381,7 +421,7 @@ gwdatanew_im_tsh = gwdatanew_im[position_in:position_end] ...@@ -381,7 +421,7 @@ gwdatanew_im_tsh = gwdatanew_im[position_in:position_end]
error_tsh=error[position_in:position_end] error_tsh=error[position_in:position_end]
# In[125]: # In[59]:
#Fitting #Fitting
...@@ -456,20 +496,34 @@ def log_probability(theta): ...@@ -456,20 +496,34 @@ def log_probability(theta):
return lp + log_likelihood(theta) return lp + log_likelihood(theta)
# In[126]: # In[60]:
dict = {'w-tau': model_dv_tau , 'w-q': model_dv_q, 'w-tau-fixed': model_dv} dict = {'w-tau': model_dv_tau , 'w-q': model_dv_q, 'w-tau-fixed': model_dv}
# In[128]: # In[64]:
#I need to provid an initial guess for 4*(nmax+1) the parameters
np.random.seed(42)
nll = lambda *args: -log_likelihood(*args)
initial = np.ones(prior_dim)
soln = minimize(nll, initial)
#x0_ml, y0_ml, a0_ml, b0_ml = soln.x
print("Maximum likelihood estimates:")
vars_ml=soln.x
print(vars_ml)
f2=dynesty.NestedSampler(log_likelihood, prior_transform, prior_dim, nlive=npoints) # In[75]:
f2=dynesty.NestedSampler(log_likelihood, prior_transform, prior_dim, nlive=npoints,sample='rwalk')
f2.run_nested() f2.run_nested()
# In[140]: # In[91]:
wstr = r'$\omega_' wstr = r'$\omega_'
...@@ -501,7 +555,7 @@ if model=='w-tau-fixed': ...@@ -501,7 +555,7 @@ if model=='w-tau-fixed':
labels = np.concatenate((amp_lab,pha_lab)) labels = np.concatenate((amp_lab,pha_lab))
# In[202]: # In[92]:
if model=='w-tau-fixed': if model=='w-tau-fixed':
...@@ -524,7 +578,7 @@ else: ...@@ -524,7 +578,7 @@ else:
npamps[i] = np.quantile(amps_aux, 0.5) npamps[i] = np.quantile(amps_aux, 0.5)
# In[203]: # In[93]:
res = f2.results res = f2.results
...@@ -533,43 +587,54 @@ res.summary() ...@@ -533,43 +587,54 @@ res.summary()
samps=f2.results.samples samps=f2.results.samples
# In[204]: # In[79]:
evidence = res.logz[-1] evidence = res.logz[-1]
evidence_error = res.logzerr[-1] evidence_error = res.logzerr[-1]
# In[205]: # In[80]:
outvalues = [nmax, simulation_number, tshift, evidence,evidence_error] summary_titles=['n','id','t_shift','dlogz','dlogz_err']
# In[206]: # In[81]:
f = open(output_folder+'/summary'+str(simulation_number)+'_'+model+'_nmax_'+str(nmax)+'.csv', 'w') f = output_folder_1+'/summary'+str(simulation_number)+'_'+model+'_nmax_'+str(nmax)+'.csv'
with f:
writer = csv.writer(f) # In[82]:
writer.writerows(map(lambda x: [x], outvalues))
if os.path.exists(f):
outvalues = np.array([[nmax, simulation_number, tshift, evidence,evidence_error]])
else:
outvalues = np.array([summary_titles,[nmax, simulation_number, tshift, evidence,evidence_error]])
with open(f, 'a') as file:
writer = csv.writer(file)
if (outvalues.shape)[0]>1 :
writer.writerows(outvalues)
else:
writer.writerow(outvalues[0])
# In[207]: # In[83]:
samps=f2.results.samples samps=f2.results.samples
samps_tr=np.transpose(samps) samps_tr=np.transpose(samps)
# In[208]: # In[84]:
sigma_vars_m = np.empty(ndim) sigma_vars_m = np.empty(prior_dim)
sigma_vars_p = np.empty(ndim) sigma_vars_p = np.empty(prior_dim)
sigma_vars = np.empty(ndim) sigma_vars = np.empty(prior_dim)
for i in range(prior_dim): for i in range(prior_dim):
amps_aux = samps_tr[i][half_points:-1] amps_aux = samps_tr[i][half_points:-1]
sigma_vars_m[i] = np.quantile(amps_aux, 0.1) sigma_vars_m[i] = np.quantile(amps_aux, 0.1)
...@@ -577,38 +642,39 @@ for i in range(prior_dim): ...@@ -577,38 +642,39 @@ for i in range(prior_dim):
sigma_vars_p[i] = np.quantile(amps_aux, 0.9) sigma_vars_p[i] = np.quantile(amps_aux, 0.9)
# In[209]: # In[85]:
sigma_vars_all = [sigma_vars,sigma_vars_m,sigma_vars_p] sigma_vars_all = [sigma_vars,sigma_vars_m,sigma_vars_p]
# In[210]: # In[86]:
sigma_vars_all=np.stack([sigma_vars,sigma_vars_m,sigma_vars_p], axis=0) sigma_vars_all=np.stack([sigma_vars,sigma_vars_m,sigma_vars_p], axis=0)
# In[211]: # In[87]:
dfslist = [pd.DataFrame(sigma_vars_all[i].reshape((-1,prior_dim)), columns=labels) for i in range(prior_dim-1)] key =['max val','lower bound','higher bound']
df2 = pd.concat(dfslist, keys=['max val','lower bound','higher bound']) file=output_folder_1+'/best_values_'+str(simulation_number)+'tshift_'+str(tshift)+'_'+model+'_nmax_'+str(nmax)+'.csv'
df2.to_csv(output_folder+'/best_values_'+str(simulation_number)+'_'+model+'_nmax_'+str(nmax)+'.csv', index = False) dfslist = [pd.DataFrame(sigma_vars_all[i].reshape((-1,prior_dim)), columns=labels, index = [key[i]]) for i in range(3)]
df2 = pd.concat(dfslist)
df2.to_csv(file, index = False)
# In[212]: # In[88]:
f = open(output_folder+'/best_sigmas_'+str(simulation_number)+'_'+model+'_'+str(nmax)+'.csv', 'w') f = open(output_folder_1+'/best_sigmas_'+str(simulation_number)+'tshift_'+str(tshift)+'_'+model+'_'+str(nmax)+'.csv', 'w')
with f: with f:
writer = csv.writer(f) writer = csv.writer(f)
writer.writerows(map(lambda x: [x], df2)) writer.writerows(map(lambda x: [x], df2))
# In[213]: # In[89]:
if model == 'w-q': if model == 'w-q':
...@@ -621,7 +687,7 @@ elif model == 'w-tau-fixed': ...@@ -621,7 +687,7 @@ elif model == 'w-tau-fixed':
truths = npamps truths = npamps
# In[216]: # In[90]:
fg, ax = dyplot.cornerplot(res, color='blue', fg, ax = dyplot.cornerplot(res, color='blue',
...@@ -633,8 +699,26 @@ fg, ax = dyplot.cornerplot(res, color='blue', ...@@ -633,8 +699,26 @@ fg, ax = dyplot.cornerplot(res, color='blue',
) )
# In[177]: # In[200]:
fg.savefig(output_folder_1+'/Dynesty_'+str(simulation_number)+'_'+model+'_nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'_chainplot.png', format = 'png', bbox_inches = 'tight')
# In[ ]:
fg.savefig(output_folder+'/Dynesty_'+str(simulation_number)+'_'+model+'_nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'_chainplot.png', format = 'png', bbox_inches = 'tight')
figband = plt.figure(figsize = (12, 9))
plt.plot(timesrd_final_tsh,gwdatanew_re_tsh, "green", alpha=0.9, lw=3, label=r'$res_{240}$')
plt.plot(timesrd_final_tsh,dict[model](vars_ml).real,'bo', alpha=0.9, lw=3, label=r'$fit$')
samples_res=samps[-20000:-1][0::100]
for i in samples_res:
plt.plot(tshift, dict[model](i).real, "r-", alpha=0.01, lw=3)
plt.title(r'Comparison of the MC fit data and the $1-\sigma$ error band')
plt.legend()
plt.xlabel('t')
plt.ylabel('h')
plt.show()
#!/usr/bin/env bash #!/usr/bin/env bash
times=(0 2.5 5 7.5 10 12 15 18 20) times=(0 0.2 0.4 0.8 1.2 2 2.5 5 7.5 10 12 15 18 20)
configfile=config.ini
pythonfile=/work/francisco.jimenez/sio/git/rdstackingproject/code_new/NR_dynesty_t0_loop.py
pythonscript=/work/francisco.jimenez/sio/pyenv/shims/python
awk -v a="10" '/^tshift/ && $3 != "supplied" { $3=a } { print }' $configfile > tmp && mv tmp $configfile
awk -v a="10" '/^tshift/ && $3 != "supplied" { $3=a } { print }' config.ini > tmp && mv tmp config.ini
for i in ${times[@]}; do for i in ${times[@]}; do
echo "submit a simulation with this parameter:" echo "submit a simulation with this parameter:"
echo "$i" echo "$i"
awk -v a="$i" '/^tshift/ && $3 != "supplied" { $3=a } { print }' config.ini > tmp && mv tmp config.ini awk -v a="$i" '/^tshift/ && $3 != "supplied" { $3=a } { print }' $configfile > tmp && mv tmp $configfile
/work/francisco.jimenez/sio/pyenv/shims/python /work/francisco.jimenez/sio/git/rdstackingproject/code_new/NR_dynesty_t0_loop.py $pythonscript $pythonfile
done done
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment