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

updates

parent 9aa1e24e
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
code_new/NR_dynesty_t0_loop.ipynb
+ 153
− 89
View file @ ae6a0d9d
Source diff could not be displayed: it is too large. Options to address this: view the blob.
code_new/NR_dynesty_t0_loop.py
+ 17
− 174
View file @ ae6a0d9d
#!/usr/bin/env python
# coding: utf-8
# In[2]:
#Import relevant modules, import data and all that
import numpy as np
from scipy import interpolate
......@@ -48,10 +42,6 @@ except SystemExit:
parser.sections()
pass
# In[3]:
# path
rootpath=parser.get('nr-paths','rootpath')
......@@ -65,18 +55,10 @@ output_folder = parser.get('output-folder','output-folder')
overwrite = parser.get('setup','overwrite')
downfactor = np.int(parser.get('setup','plot_down_factor'))
# In[4]:
if not os.path.exists(output_folder):
os.mkdir(output_folder)
print("Directory " , output_folder , " Created ")
# In[5]:
# time config
tshift=parser.get('time-setup','tshift')
......@@ -88,10 +70,6 @@ tend = np.float(tend)
t_align=parser.get('time-setup','t_align')
t_align = np.float(t_align)
# In[6]:
# n-tones & nlive
nmax=parser.get('n-tones','nmax')
......@@ -100,10 +78,6 @@ nmax = np.int(nmax)
npoints=parser.get('n-live-points','npoints')
npoints = np.int(npoints)
# In[7]:
# model
model=parser.get('rd-model','model')
error_str = eval(parser.get('rd-model','error_str'))
......@@ -119,18 +93,17 @@ print('nmax:',nmax)
print('tshift:',tshift)
print('error:', error_str)
# In[8]:
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 ")
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[9]:
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'
# loading priors
w_mins=np.empty(nmax+1)
......@@ -177,10 +150,6 @@ if model == 'w-tau-fixed':
priors_max = np.concatenate((a_maxs,ph_maxs))
prior_dim = len(priors_min)
# In[10]:
vary_fund = True
#sampler parameters
......@@ -225,10 +194,6 @@ 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[11]:
#This loads the 22 mode data
gw = {}
gw[simulation_number] = h5py.File(simulation_path_1, 'r')
......@@ -257,10 +222,6 @@ times5 = gw5_sxs_bbh_0305[:,0]
tmax5=FindTmaximum(gw5_sxs_bbh_0305)
times5 = times5 - tmax5
# In[12]:
#Select the data from 0 onwards
position = np.argmax(times >= (t_align))
position5 = np.argmax(times5 >= (t_align))
......@@ -269,10 +230,6 @@ 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[13]:
#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')
......@@ -281,10 +238,6 @@ 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[14]:
#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$')
......@@ -293,29 +246,17 @@ 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[15]:
# 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[16]:
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[17]:
if timesrd5[-1]>= timesrd[-1]:
timesrd_final = timesrd
else:
......@@ -329,18 +270,10 @@ gwdatanew_im5 = gwnew_im5(timesrd_final)
gwdatanew = gwdatanew_re - 1j*gwdatanew_im
gwdatanew5 = gwdatanew_re5- 1j*gwdatanew_im5
# In[18]:
mismatch=1-EasyMatchT(timesrd_final,gwdatanew,gwdatanew5,0,0+90)
error=np.sqrt(2*mismatch)
print(mismatch)
# In[19]:
# Phase alignement
phas = np.angle(gwdatanew)
phas = np.unwrap(phas)
......@@ -349,10 +282,6 @@ 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[20]:
position = np.argmax(timesrd_final >= (t_align))
dphase = phas5[position]-phas[position]
print(dphase)
......@@ -368,10 +297,6 @@ 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[21]:
mismatch=1-EasyMatchT(timesrd_final,gwdatanew,gwdatanew5,0,+90)
print(mismatch)
if error_str:
......@@ -379,10 +304,6 @@ if error_str:
else :
error = 1
# In[22]:
#Test the new interpolated data
plt.figure(figsize = (12, 8))
plt.plot(timesrd_final, gwdatanew.real, "r", alpha=0.3, lw=2, label='Lev6')
......@@ -390,10 +311,6 @@ 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[23]:
#Test the error data
plt.figure(figsize = (12, 8))
plt.plot(timesrd_final, error.real, "b", alpha=0.3, lw=2, label='error real')
......@@ -401,10 +318,6 @@ 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[24]:
#Take the piece of waveform you want
position_in = np.argmax(timesrd_final >= tshift)
position_end = np.argmax(timesrd_final >= tend)
......@@ -413,10 +326,6 @@ 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[25]:
#Fitting
#RD model for nmax tones. Amplitudes are in (xn*Exp[i yn]) version. Used here.
def model_dv_q(theta):
......@@ -488,16 +397,8 @@ def log_probability(theta):
return -np.inf
return lp + log_likelihood(theta)
# In[26]:
dict = {'w-tau': model_dv_tau , 'w-q': model_dv_q, 'w-tau-fixed': model_dv}
# In[27]:
#I need to provid an initial guess for 4*(nmax+1) the parameters
np.random.seed(42)
nll = lambda *args: -log_likelihood(*args)
......@@ -508,17 +409,9 @@ print("Maximum likelihood estimates:")
vars_ml=soln.x
print(vars_ml)
# In[28]:
f2=dynesty.NestedSampler(log_likelihood, prior_transform, prior_dim, nlive=npoints,sample='rwalk')
f2.run_nested()
# In[29]:
wstr = r'$\omega_'
if model == 'w-tau':
......@@ -547,10 +440,6 @@ labels = np.concatenate((w_lab,tau_lab,amp_lab,pha_lab))
if model=='w-tau-fixed':
labels = np.concatenate((amp_lab,pha_lab))
# In[30]:
if model=='w-tau-fixed':
rg = (nmax+1)
else:
......@@ -570,56 +459,31 @@ else:
amps_aux = samps_tr[i][half_points:-1]
npamps[i] = np.quantile(amps_aux, 0.5)
# In[31]:
res = f2.results
res.samples_u.shape
res.summary()
samps=f2.results.samples
# In[32]:
evidence = res.logz[-1]
evidence_error = res.logzerr[-1]
# In[33]:
summary_titles=['n','id','t_shift','dlogz','dlogz_err']
f = output_folder_1+'/summary'+str(simulation_number)+'_'+model+'_nmax_'+str(nmax)+'.csv'
# In[35]:
if os.path.exists(f):
if os.path.exists(sumary_data):
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:
with open(sumary_data, 'a') as file:
writer = csv.writer(file)
if (outvalues.shape)[0]>1 :
writer.writerows(outvalues)
else:
writer.writerow(outvalues[0])
# In[36]:
samps=f2.results.samples
samps_tr=np.transpose(samps)
# In[37]:
sigma_vars_m = np.empty(prior_dim)
sigma_vars_p = np.empty(prior_dim)
sigma_vars = np.empty(prior_dim)
......@@ -629,29 +493,16 @@ 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[38]:
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[69]:
key =['max val','lower bound','higher bound']
file=output_folder_1+'/best_values_'+str(simulation_number)+'_'+model+'_nmax_'+str(nmax)+'.csv'
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)
if os.path.exists(file):
df2.to_csv(file, mode='a', header=False,index = True)
if os.path.exists(best_data):
df2.to_csv(best_data, mode='a', header=False,index = True)
else:
df2.to_csv(file,index = True)
# In[42]:
df2.to_csv(best_data,index = True)
if model == 'w-q':
tau_val = np.pi*w*tau
......@@ -662,10 +513,6 @@ elif model == 'w-tau':
elif model == 'w-tau-fixed':
truths = npamps
# In[43]:
fg, ax = dyplot.cornerplot(res, color='blue',
show_titles=True,
labels=labels,
......@@ -674,15 +521,15 @@ fg, ax = dyplot.cornerplot(res, color='blue',
truth_color='red',
)
fg.savefig(corner_plot, format = 'png', bbox_inches = 'tight')
# In[127]:
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')
from dynesty import plotting as dyplot
# In[148]:
lnz_truth = ndim * -np.log(2 * 10.) # analytic evidence solution
fig, axes = dyplot.runplot(res, lnz_truth=lnz_truth)
fig.tight_layout()
fig.savefig(diagnosis_plot, format = 'png', dpi = 384, 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}$')
......@@ -698,10 +545,6 @@ plt.xlabel('t')
plt.ylabel('h')
plt.show()
# In[149]:
fit_plot=output_folder_1+'/Fit_results_'+str(simulation_number)+'tshift_'+str(tshift)+'_'+model+'_nmax_'+str(nmax)+'.png'
figband.savefig(fit_plot)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment

AEI Privacy Policy Imprint