added error options

code_new/NR_dynesty_t0_loop.py

 #!/usr/bin/env python
 # coding: utf-8
 
-# In[3]:
+# In[48]:
 
 
 #Import relevant modules, import data and all that
@@ -49,7 +49,7 @@ except SystemExit:
     pass
 
 
-# In[4]:
+# In[49]:
 
 
 # path
@@ -64,9 +64,10 @@ simulation_number = np.int(simulation_number)
 output_folder = parser.get('output-folder','output-folder')
 overwrite = parser.get('setup','overwrite')
 downfactor = np.int(parser.get('setup','plot_down_factor'))
+sampler = parser.get('setup','sampler')
 
 
-# In[5]:
+# In[50]:
 
 
 if not os.path.exists(output_folder):
@@ -74,7 +75,7 @@ if not os.path.exists(output_folder):
     print("Directory " , output_folder ,  " Created ")
 
 
-# In[6]:
+# In[51]:
 
 
 # time config
@@ -89,7 +90,7 @@ t_align=parser.get('time-setup','t_align')
 t_align = np.float(t_align)
 
 
-# In[7]:
+# In[52]:
 
 
 # n-tones & nlive
@@ -101,7 +102,7 @@ npoints=parser.get('n-live-points','npoints')
 npoints = np.int(npoints)
 
 
-# In[8]:
+# In[81]:
 
 
 # model
@@ -120,16 +121,16 @@ print('tshift:',tshift)
 print('error:', error_str)
 
 
-# In[10]:
+# In[83]:
 
 
-output_folder_1=output_folder+'/'+model+'-nmax'+str(nmax)
+output_folder_1=output_folder+'/'+model+'-nmax'+str(nmax)+'_'+str(error_str)
 if not os.path.exists(output_folder_1):
     os.mkdir(output_folder_1)
     print("Directory " , output_folder_1 ,  " Created ")
 
 
-# In[11]:
+# In[84]:
 
 
 corner_plot=output_folder_1+'/Dynesty_'+str(simulation_number)+'_'+model+'_nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'corner_plot.png'
@@ -137,14 +138,14 @@ diagnosis_plot=output_folder_1+'/Dynesty_diagnosis'+str(simulation_number)+'_'+m
 fit_plot=output_folder_1+'/Fit_results_'+str(simulation_number)+'tshift_'+str(tshift)+'_'+model+'_nmax_'+str(nmax)+'.png'
 
 
-# In[12]:
+# In[85]:
 
 
 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]:
+# In[86]:
 
 
 # loading priors
@@ -193,7 +194,7 @@ if model == 'w-tau-fixed':
     prior_dim = len(priors_min)
 
 
-# In[67]:
+# In[87]:
 
 
 vary_fund = True
@@ -241,7 +242,7 @@ def tauRD_to_t_Phys(tau,M):
     return ((M*MS*G)/c**3)*tau
 
 
-# In[68]:
+# In[88]:
 
 
 #This loads the 22 mode data
@@ -273,7 +274,7 @@ tmax5=FindTmaximum(gw5_sxs_bbh_0305)
 times5 = times5 - tmax5
 
 
-# In[69]:
+# In[89]:
 
 
 #Select the data from 0 onwards
@@ -285,7 +286,7 @@ timesrd=gw_sxs_bbh_0305[position:-1][:,0][:-1]-tmax
 timesrd5=gw5_sxs_bbh_0305[position5:-1][:,0][:-1]-tmax5
 
 
-# In[70]:
+# In[90]:
 
 
 #Test plot real part (data was picked in the last cell). Aligning in time
@@ -297,7 +298,7 @@ plt.plot(timesrd5, np.sqrt(gw_sxs_bbh_0305rd5[:,1]**2+gw_sxs_bbh_0305rd5[:,2]**2
 plt.legend()
 
 
-# In[71]:
+# In[91]:
 
 
 #Test plot im part (data was picked in the last cell). Aligning in time
@@ -309,7 +310,7 @@ plt.plot(timesrd5, np.sqrt(gw_sxs_bbh_0305rd5[:,1]**2+gw_sxs_bbh_0305rd5[:,2]**2
 plt.legend()
 
 
-# In[72]:
+# In[92]:
 
 
 # Depending on nmax, you load nmax number of freqs. and damping times from the qnm package
@@ -318,7 +319,7 @@ w = (np.real(omegas))/mf
 tau=-1/(np.imag(omegas))*mf
 
 
-# In[73]:
+# In[93]:
 
 
 gwnew_re = interpolate.interp1d(timesrd, gw_sxs_bbh_0305rd[:,1], kind = 'cubic')
@@ -328,7 +329,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')
 
 
-# In[74]:
+# In[94]:
 
 
 if timesrd5[-1]>= timesrd[-1]: 
@@ -345,7 +346,7 @@ gwdatanew = gwdatanew_re - 1j*gwdatanew_im
 gwdatanew5 = gwdatanew_re5- 1j*gwdatanew_im5
 
 
-# In[75]:
+# In[95]:
 
 
 mismatch=1-EasyMatchT(timesrd_final,gwdatanew,gwdatanew5,0,0+90)
@@ -353,7 +354,7 @@ error=np.sqrt(2*mismatch)
 print(mismatch)
 
 
-# In[76]:
+# In[96]:
 
 
 # Phase alignement
@@ -365,7 +366,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$')
 
 
-# In[77]:
+# In[97]:
 
 
 position = np.argmax(timesrd_final >= (t_align))
@@ -384,7 +385,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$')
 
 
-# In[78]:
+# In[98]:
 
 
 mismatch=1-EasyMatchT(timesrd_final,gwdatanew,gwdatanew5,0,+90)
@@ -395,10 +396,11 @@ else :
     error = 1 
 
 
-# In[79]:
+# In[100]:
 
 
 #Test the new interpolated data
+if error_str:
     plt.figure(figsize = (12, 8))
     plt.plot(timesrd_final, gwdatanew.real, "r", alpha=0.3, lw=2, label='Lev6')
     plt.plot(timesrd_final, gwdatanew5.real, "b", alpha=0.3, lw=2, label='Lev5')
@@ -406,10 +408,11 @@ plt.plot(timesrd_final, error.real, "b", alpha=0.3, lw=2, label='error')
     plt.legend()
 
 
-# In[80]:
+# In[101]:
 
 
 #Test the error data
+if error_str:
     plt.figure(figsize = (12, 8))
     plt.plot(timesrd_final, error.real, "b", alpha=0.3, lw=2, label='error real')
     plt.plot(timesrd_final, error.imag, "r", alpha=0.3, lw=2, label='error imag')
@@ -417,7 +420,7 @@ plt.plot(timesrd_final, np.sqrt(error.imag**2+error.real**2), "r", alpha=0.3, lw
     plt.legend()
 
 
-# In[81]:
+# In[103]:
 
 
 #Take the piece of waveform you want
@@ -426,10 +429,13 @@ position_end = np.argmax(timesrd_final >= tend)
 timesrd_final_tsh = timesrd_final[position_in:position_end]
 gwdatanew_re_tsh = gwdatanew_re[position_in:position_end]
 gwdatanew_im_tsh = gwdatanew_im[position_in:position_end]
+if error_str:
     error_tsh=error[position_in:position_end]
+else:
+    error_tsh=1
 
 
-# In[82]:
+# In[104]:
 
 
 #Fitting
@@ -504,13 +510,13 @@ def log_probability(theta):
     return lp + log_likelihood(theta)
 
 
-# In[83]:
+# In[105]:
 
 
 dict = {'w-tau': model_dv_tau , 'w-q': model_dv_q, 'w-tau-fixed': model_dv}
 
 
-# In[84]:
+# In[106]:
 
 
 #I need to provid an initial guess for 4*(nmax+1) the parameters
@@ -524,14 +530,14 @@ vars_ml=soln.x
 print(vars_ml)
 
 
-# In[85]:
+# In[107]:
 
 
-f2=dynesty.NestedSampler(log_likelihood, prior_transform, prior_dim, nlive=npoints,sample='rwalk')
+f2=dynesty.NestedSampler(log_likelihood, prior_transform, prior_dim, nlive=npoints,sample=sampler)
 f2.run_nested()
 
 
-# In[87]:
+# In[108]:
 
 
 wstr = r'$\omega_'
@@ -563,7 +569,7 @@ if model=='w-tau-fixed':
     labels = np.concatenate((amp_lab,pha_lab))
 
 
-# In[88]:
+# In[109]:
 
 
 if model=='w-tau-fixed':
@@ -586,7 +592,7 @@ else:
         npamps[i] = np.quantile(amps_aux, 0.5)
 
 
-# In[89]:
+# In[110]:
 
 
 res = f2.results
@@ -595,20 +601,20 @@ res.summary()
 samps=f2.results.samples
 
 
-# In[90]:
+# In[111]:
 
 
 evidence = res.logz[-1]
 evidence_error = res.logzerr[-1]
 
 
-# In[91]:
+# In[112]:
 
 
 summary_titles=['n','id','t_shift','dlogz','dlogz_err']
 
 
-# In[92]:
+# In[113]:
 
 
 if os.path.exists(sumary_data):
@@ -624,14 +630,14 @@ with open(sumary_data, 'a') as file:
         writer.writerow(outvalues[0])
 
 
-# In[93]:
+# In[114]:
 
 
 samps=f2.results.samples
 samps_tr=np.transpose(samps)
 
 
-# In[94]:
+# In[115]:
 
 
 sigma_vars_m = np.empty(prior_dim)
@@ -644,18 +650,18 @@ for i in range(prior_dim):
     sigma_vars_p[i] = np.quantile(amps_aux, 0.9)
 
 
-# In[95]:
+# In[116]:
 
 
 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]:
+# In[117]:
 
 
 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)]
+dfslist = [pd.DataFrame(np.concatenate(([tshift],sigma_vars_all[i])).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(best_data):
     df2.to_csv(best_data, mode='a', header=False,index = True)
@@ -663,7 +669,7 @@ else:
     df2.to_csv(best_data,index = True)
 
 
-# In[97]:
+# In[118]:
 
 
 if model == 'w-q':
@@ -676,7 +682,7 @@ elif model == 'w-tau-fixed':
     truths = npamps
 
 
-# In[98]:
+# In[119]:
 
 
 fg, ax = dyplot.cornerplot(res, color='blue', 
@@ -688,13 +694,13 @@ fg, ax = dyplot.cornerplot(res, color='blue',
 )
 
 
-# In[99]:
+# In[121]:
 
 
 fg.savefig(corner_plot, format = 'png', bbox_inches = 'tight')
 
 
-# In[100]:
+# In[122]:
 
 
 from dynesty import plotting as dyplot
@@ -704,13 +710,13 @@ fig, axes = dyplot.runplot(res, lnz_truth=lnz_truth)
 fig.tight_layout()
 
 
-# In[44]:
+# In[123]:
 
 
 fig.savefig(diagnosis_plot, format = 'png', dpi = 384, bbox_inches = 'tight')
 
 
-# In[101]:
+# In[124]:
 
 
 figband = plt.figure(figsize = (12, 9))
@@ -720,7 +726,7 @@ 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, dict[model](sample).real, "r-", alpha=0.01, lw=3)
+    plt.plot(timesrd_final_tsh, dict[model](sample).real, "r-", alpha=0.04, lw=3)
 plt.title(r'Comparison of the MC fit data and the $1-\sigma$ error band')
 plt.legend()
 plt.xlabel('t')
@@ -728,7 +734,7 @@ plt.ylabel('h')
 plt.show()
 
 
-# In[102]:
+# In[125]:
 
 
 figband.savefig(fit_plot)