From c11167edb4a871b1719a300682c0365285ffa19e Mon Sep 17 00:00:00 2001
From: frcojimenez <frco.jimenez@gmail.com>
Date: Fri, 20 Nov 2020 16:54:30 +0000
Subject: [PATCH] fixing bug on plots
---
code_new/NR_dynesty_t0_loop.ipynb | 14 ++-
code_new/NR_dynesty_t0_loop.py | 193 +++++++++++++++++++++++++++++-
2 files changed, 205 insertions(+), 2 deletions(-)
diff --git a/code_new/NR_dynesty_t0_loop.ipynb b/code_new/NR_dynesty_t0_loop.ipynb
index 27bcbdd..4e22d9a 100644
--- a/code_new/NR_dynesty_t0_loop.ipynb
+++ b/code_new/NR_dynesty_t0_loop.ipynb
@@ -1121,7 +1121,7 @@
"samples_1sigma = filter(lambda sample: np.all(onesig_bounds[0] <= sample) and np.all(sample <= onesig_bounds[1]), samps)\n",
"samples_1sigma_down = list(samples_1sigma)[::downfactor]\n",
"for sample in samples_1sigma_down:\n",
- " plt.plot(timesrd_final_tsh, model_dv(sample).real, \"r-\", alpha=0.01, lw=3)\n",
+ " plt.plot(timesrd_final_tsh, dict[model](sample).real, \"r-\", alpha=0.01, lw=3)\n",
"plt.title(r'Comparison of the MC fit data and the $1-\\sigma$ error band')\n",
"plt.legend()\n",
"plt.xlabel('t')\n",
@@ -1151,6 +1151,18 @@
"display_name": "Python 3",
"language": "python",
"name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.7.7"
}
},
"nbformat": 4,
diff --git a/code_new/NR_dynesty_t0_loop.py b/code_new/NR_dynesty_t0_loop.py
index 8e1537b..d074f6d 100755
--- a/code_new/NR_dynesty_t0_loop.py
+++ b/code_new/NR_dynesty_t0_loop.py
@@ -1,3 +1,9 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# In[3]:
+
+
#Import relevant modules, import data and all that
import numpy as np
from scipy import interpolate
@@ -42,6 +48,10 @@ except SystemExit:
parser.sections()
pass
+
+# In[4]:
+
+
# path
rootpath=parser.get('nr-paths','rootpath')
@@ -55,10 +65,18 @@ output_folder = parser.get('output-folder','output-folder')
overwrite = parser.get('setup','overwrite')
downfactor = np.int(parser.get('setup','plot_down_factor'))
+
+# In[5]:
+
+
if not os.path.exists(output_folder):
os.mkdir(output_folder)
print("Directory " , output_folder , " Created ")
+
+# In[6]:
+
+
# time config
tshift=parser.get('time-setup','tshift')
@@ -70,6 +88,10 @@ tend = np.float(tend)
t_align=parser.get('time-setup','t_align')
t_align = np.float(t_align)
+
+# In[7]:
+
+
# n-tones & nlive
nmax=parser.get('n-tones','nmax')
@@ -78,6 +100,10 @@ nmax = np.int(nmax)
npoints=parser.get('n-live-points','npoints')
npoints = np.int(npoints)
+
+# In[8]:
+
+
# model
model=parser.get('rd-model','model')
error_str = eval(parser.get('rd-model','error_str'))
@@ -93,18 +119,34 @@ print('nmax:',nmax)
print('tshift:',tshift)
print('error:', error_str)
+
+# In[10]:
+
+
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[11]:
+
+
corner_plot=output_folder_1+'/Dynesty_'+str(simulation_number)+'_'+model+'_nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'corner_plot.png'
diagnosis_plot=output_folder_1+'/Dynesty_diagnosis'+str(simulation_number)+'_'+model+'_nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'.png'
fit_plot=output_folder_1+'/Fit_results_'+str(simulation_number)+'tshift_'+str(tshift)+'_'+model+'_nmax_'+str(nmax)+'.png'
+
+# In[12]:
+
+
sumary_data = output_folder_1+'/summary'+str(simulation_number)+'_'+model+'_nmax_'+str(nmax)+'.csv'
best_data=output_folder_1+'/best_values_'+str(simulation_number)+'_'+model+'_nmax_'+str(nmax)+'.csv'
+
+# In[66]:
+
+
# loading priors
w_mins=np.empty(nmax+1)
w_maxs=np.empty(nmax+1)
@@ -150,6 +192,10 @@ if model == 'w-tau-fixed':
priors_max = np.concatenate((a_maxs,ph_maxs))
prior_dim = len(priors_min)
+
+# In[67]:
+
+
vary_fund = True
#sampler parameters
@@ -194,6 +240,10 @@ 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[68]:
+
+
#This loads the 22 mode data
gw = {}
gw[simulation_number] = h5py.File(simulation_path_1, 'r')
@@ -222,6 +272,10 @@ times5 = gw5_sxs_bbh_0305[:,0]
tmax5=FindTmaximum(gw5_sxs_bbh_0305)
times5 = times5 - tmax5
+
+# In[69]:
+
+
#Select the data from 0 onwards
position = np.argmax(times >= (t_align))
position5 = np.argmax(times5 >= (t_align))
@@ -230,6 +284,10 @@ gw_sxs_bbh_0305rd5=gw5_sxs_bbh_0305[position5+1:-1]
timesrd=gw_sxs_bbh_0305[position:-1][:,0][:-1]-tmax
timesrd5=gw5_sxs_bbh_0305[position5:-1][:,0][:-1]-tmax5
+
+# In[70]:
+
+
#Test plot real part (data was picked in the last cell). Aligning in time
plt.figure(figsize = (12, 8))
plt.plot(timesrd, gw_sxs_bbh_0305rd[:,1], "r", alpha=0.3, lw=3, label=r'$Lev6$: real')
@@ -238,6 +296,10 @@ plt.plot(timesrd5, gw_sxs_bbh_0305rd5[:,1], "b", alpha=0.3, lw=3, label=r'$Lev5:
plt.plot(timesrd5, np.sqrt(gw_sxs_bbh_0305rd5[:,1]**2+gw_sxs_bbh_0305rd5[:,2]**2), "b", alpha=0.3, lw=3, label=r'$Lev5\,amp$')
plt.legend()
+
+# In[71]:
+
+
#Test plot im part (data was picked in the last cell). Aligning in time
plt.figure(figsize = (12, 8))
plt.plot(timesrd, gw_sxs_bbh_0305rd[:,2], "r", alpha=0.3, lw=3, label=r'$Lev6: imag$')
@@ -246,17 +308,29 @@ plt.plot(timesrd5, gw_sxs_bbh_0305rd5[:,2], "b", alpha=0.3, lw=3, label=r'$Lev5:
plt.plot(timesrd5, np.sqrt(gw_sxs_bbh_0305rd5[:,1]**2+gw_sxs_bbh_0305rd5[:,2]**2), "b", alpha=0.3, lw=3, label=r'$Lev5\,amp$')
plt.legend()
+
+# In[72]:
+
+
# Depending on nmax, you load nmax number of freqs. and damping times from the qnm package
omegas = [qnm.modes_cache(s=-2,l=2,m=2,n=i)(a=af)[0] for i in range (0,dim)]
w = (np.real(omegas))/mf
tau=-1/(np.imag(omegas))*mf
+
+# In[73]:
+
+
gwnew_re = interpolate.interp1d(timesrd, gw_sxs_bbh_0305rd[:,1], kind = 'cubic')
gwnew_im = interpolate.interp1d(timesrd, gw_sxs_bbh_0305rd[:,2], kind = 'cubic')
gwnew_re5 = interpolate.interp1d(timesrd5, gw_sxs_bbh_0305rd5[:,1], kind = 'cubic')
gwnew_im5 = interpolate.interp1d(timesrd5, gw_sxs_bbh_0305rd5[:,2], kind = 'cubic')
+
+# In[74]:
+
+
if timesrd5[-1]>= timesrd[-1]:
timesrd_final = timesrd
else:
@@ -270,10 +344,18 @@ gwdatanew_im5 = gwnew_im5(timesrd_final)
gwdatanew = gwdatanew_re - 1j*gwdatanew_im
gwdatanew5 = gwdatanew_re5- 1j*gwdatanew_im5
+
+# In[75]:
+
+
mismatch=1-EasyMatchT(timesrd_final,gwdatanew,gwdatanew5,0,0+90)
error=np.sqrt(2*mismatch)
print(mismatch)
+
+# In[76]:
+
+
# Phase alignement
phas = np.angle(gwdatanew)
phas = np.unwrap(phas)
@@ -282,6 +364,10 @@ phas5 = np.unwrap(phas5)
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$')
+
+# In[77]:
+
+
position = np.argmax(timesrd_final >= (t_align))
dphase = phas5[position]-phas[position]
print(dphase)
@@ -297,6 +383,10 @@ phas5 = np.unwrap(phas5)
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$')
+
+# In[78]:
+
+
mismatch=1-EasyMatchT(timesrd_final,gwdatanew,gwdatanew5,0,+90)
print(mismatch)
if error_str:
@@ -304,6 +394,10 @@ if error_str:
else :
error = 1
+
+# In[79]:
+
+
#Test the new interpolated data
plt.figure(figsize = (12, 8))
plt.plot(timesrd_final, gwdatanew.real, "r", alpha=0.3, lw=2, label='Lev6')
@@ -311,6 +405,10 @@ plt.plot(timesrd_final, gwdatanew5.real, "b", alpha=0.3, lw=2, label='Lev5')
plt.plot(timesrd_final, error.real, "b", alpha=0.3, lw=2, label='error')
plt.legend()
+
+# In[80]:
+
+
#Test the error data
plt.figure(figsize = (12, 8))
plt.plot(timesrd_final, error.real, "b", alpha=0.3, lw=2, label='error real')
@@ -318,6 +416,10 @@ plt.plot(timesrd_final, error.imag, "r", alpha=0.3, lw=2, label='error imag')
plt.plot(timesrd_final, np.sqrt(error.imag**2+error.real**2), "r", alpha=0.3, lw=2, label='all error')
plt.legend()
+
+# In[81]:
+
+
#Take the piece of waveform you want
position_in = np.argmax(timesrd_final >= tshift)
position_end = np.argmax(timesrd_final >= tend)
@@ -326,6 +428,10 @@ gwdatanew_re_tsh = gwdatanew_re[position_in:position_end]
gwdatanew_im_tsh = gwdatanew_im[position_in:position_end]
error_tsh=error[position_in:position_end]
+
+# In[82]:
+
+
#Fitting
#RD model for nmax tones. Amplitudes are in (xn*Exp[i yn]) version. Used here.
def model_dv_q(theta):
@@ -397,8 +503,16 @@ def log_probability(theta):
return -np.inf
return lp + log_likelihood(theta)
+
+# In[83]:
+
+
dict = {'w-tau': model_dv_tau , 'w-q': model_dv_q, 'w-tau-fixed': model_dv}
+
+# In[84]:
+
+
#I need to provid an initial guess for 4*(nmax+1) the parameters
np.random.seed(42)
nll = lambda *args: -log_likelihood(*args)
@@ -409,9 +523,17 @@ print("Maximum likelihood estimates:")
vars_ml=soln.x
print(vars_ml)
+
+# In[85]:
+
+
f2=dynesty.NestedSampler(log_likelihood, prior_transform, prior_dim, nlive=npoints,sample='rwalk')
f2.run_nested()
+
+# In[87]:
+
+
wstr = r'$\omega_'
if model == 'w-tau':
@@ -440,6 +562,10 @@ labels = np.concatenate((w_lab,tau_lab,amp_lab,pha_lab))
if model=='w-tau-fixed':
labels = np.concatenate((amp_lab,pha_lab))
+
+# In[88]:
+
+
if model=='w-tau-fixed':
rg = (nmax+1)
else:
@@ -459,16 +585,32 @@ else:
amps_aux = samps_tr[i][half_points:-1]
npamps[i] = np.quantile(amps_aux, 0.5)
+
+# In[89]:
+
+
res = f2.results
res.samples_u.shape
res.summary()
samps=f2.results.samples
+
+# In[90]:
+
+
evidence = res.logz[-1]
evidence_error = res.logzerr[-1]
+
+# In[91]:
+
+
summary_titles=['n','id','t_shift','dlogz','dlogz_err']
+
+# In[92]:
+
+
if os.path.exists(sumary_data):
outvalues = np.array([[nmax, simulation_number, tshift, evidence,evidence_error]])
else:
@@ -481,9 +623,17 @@ with open(sumary_data, 'a') as file:
else:
writer.writerow(outvalues[0])
+
+# In[93]:
+
+
samps=f2.results.samples
samps_tr=np.transpose(samps)
+
+# In[94]:
+
+
sigma_vars_m = np.empty(prior_dim)
sigma_vars_p = np.empty(prior_dim)
sigma_vars = np.empty(prior_dim)
@@ -493,9 +643,17 @@ for i in range(prior_dim):
sigma_vars[i] = np.quantile(amps_aux, 0.5)
sigma_vars_p[i] = np.quantile(amps_aux, 0.9)
+
+# In[95]:
+
+
sigma_vars_all = [sigma_vars,sigma_vars_m,sigma_vars_p]
sigma_vars_all=np.stack([sigma_vars,sigma_vars_m,sigma_vars_p], axis=0)
+
+# In[96]:
+
+
key =['max val','lower bound','higher bound']
dfslist = [pd.DataFrame(np.concatenate(([tshift],sigma_vars_all[1])).reshape((-1,prior_dim+1)), columns=np.concatenate((['tshift'],labels)), index = [key[i]]) for i in range(3)]
df2 = pd.concat(dfslist)
@@ -504,6 +662,10 @@ if os.path.exists(best_data):
else:
df2.to_csv(best_data,index = True)
+
+# In[97]:
+
+
if model == 'w-q':
tau_val = np.pi*w*tau
truths = np.concatenate((w,tau_val,npamps))
@@ -513,6 +675,10 @@ elif model == 'w-tau':
elif model == 'w-tau-fixed':
truths = npamps
+
+# In[98]:
+
+
fg, ax = dyplot.cornerplot(res, color='blue',
show_titles=True,
labels=labels,
@@ -521,16 +687,32 @@ fg, ax = dyplot.cornerplot(res, color='blue',
truth_color='red',
)
+
+# In[99]:
+
+
fg.savefig(corner_plot, format = 'png', bbox_inches = 'tight')
+
+# In[100]:
+
+
from dynesty import plotting as dyplot
lnz_truth = ndim * -np.log(2 * 10.) # analytic evidence solution
fig, axes = dyplot.runplot(res, lnz_truth=lnz_truth)
fig.tight_layout()
+
+# In[44]:
+
+
fig.savefig(diagnosis_plot, format = 'png', dpi = 384, bbox_inches = 'tight')
+
+# In[101]:
+
+
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$')
@@ -538,13 +720,22 @@ onesig_bounds = np.array([np.percentile(samps[:, i], [16, 84]) for i in range(le
samples_1sigma = filter(lambda sample: np.all(onesig_bounds[0] <= sample) and np.all(sample <= onesig_bounds[1]), samps)
samples_1sigma_down = list(samples_1sigma)[::downfactor]
for sample in samples_1sigma_down:
- plt.plot(timesrd_final_tsh, model_dv(sample).real, "r-", alpha=0.01, lw=3)
+ plt.plot(timesrd_final_tsh, dict[model](sample).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()
+
+# In[102]:
+
+
figband.savefig(fit_plot)
+# In[ ]:
+
+
+
+
--
GitLab