python notebook

code_new/NR_dynesty_t0_loop.py

@@ -82,7 +82,7 @@ try:
     parser.sections()
 except SystemExit: 
     parser = ConfigParser()
-    parser.read('config_n2_q10.ini')
+    parser.read('config_fixed_n1_m_af.ini')
     parser.sections()
     pass
 
@@ -331,26 +331,26 @@ def QNM_Berti(mf,af,l,m):
     tau_ma_a=[None]*(nmax+1)
     
     for i in range(nmax+1):
-        qnm=rdowndata[0,1:3,position]
+        qnm=rdowndata[i,1:3,position]
         w_m_a[i] = qnm[0]/mf
         tau_ma_a[i] = -1/(qnm[1])*mf
 
     return w_m_a, tau_ma_a
 
 
-# In[14]:
+# In[15]:
 
 
 np.sqrt(12/2*1/tauRD_to_t_NR(1.3*10**-47,70)*(5*10**(-21))**2)
 
 
-# In[15]:
+# In[16]:
 
 
 np.sqrt(0.004/2*1/(1.3*10**-47)*(5*10**(-21))**2)
 
 
-# In[16]:
+# In[17]:
 
 
 gw = {}
@@ -406,7 +406,7 @@ elif nr_code=='LaZeV':
     times=times_1
 
 
-# In[17]:
+# In[18]:
 
 
 if nr_code=='SXS':
@@ -436,7 +436,7 @@ tmax5=FindTmaximum(gw5_sxs_bbh_0305[round(len(gw_sxs_bbh_0305)/2):])
 times5 = times5 - tmax5
 
 
-# In[18]:
+# In[19]:
 
 
 if parser.has_option('setup','qnm_model'):
@@ -449,7 +449,7 @@ else:
     w , tau = QNM_spectrum(mf,af,2,2)
 
 
-# In[19]:
+# In[20]:
 
 
 # loading priors
@@ -509,7 +509,7 @@ elif model ==  'w-tau-fixed-m-af':
     priors=np.column_stack((priors_min,priors_max))
 
 
-# In[20]:
+# In[21]:
 
 
 #Select the data from 0 onwards
@@ -521,7 +521,7 @@ timesrd=gw_sxs_bbh_0305[position:-1][:,0][:]-tmax
 timesrd5=gw5_sxs_bbh_0305[position5:-1][:,0][:]-tmax5
 
 
-# In[21]:
+# In[22]:
 
 
 #Test plot real part (data was picked in the last cell). Aligning in time
@@ -533,13 +533,13 @@ plt.plot(timesrd5, np.sqrt(gw_sxs_bbh_0305rd5[:,1]**2+gw_sxs_bbh_0305rd5[:,2]**2
 plt.legend()
 
 
-# In[22]:
+# In[23]:
 
 
 #[plt.errorbar(csv_data_fixed[i]['t_shift'], -csv_data_fixed[i]['dlogz'], yerr=csv_data_fixed[i]['dlogz_err'], fmt='o',color=colors[i],label =tags_fixed[i]) for i in range(len(csv_data_fixed))] 
 
 
-# In[23]:
+# In[24]:
 
 
 gwnew_re = interpolate.interp1d(timesrd, gw_sxs_bbh_0305rd[:,1], kind = 'cubic')
@@ -549,7 +549,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[24]:
+# In[25]:
 
 
 if timesrd5[-1]>= timesrd[-1]: 
@@ -566,7 +566,7 @@ gwdatanew = gwdatanew_re - 1j*gwdatanew_im
 gwdatanew5 = gwdatanew_re5- 1j*gwdatanew_im5
 
 
-# In[25]:
+# In[26]:
 
 
 #taus, corr= twopoint_autocovariance(timesrd,gwdatanew-gwdatanew5)
@@ -578,7 +578,7 @@ gwdatanew5 = gwdatanew_re5- 1j*gwdatanew_im5
 #taus[index]
 
 
-# In[26]:
+# In[27]:
 
 
 mismatch=1-EasyMatchT(timesrd_final,gwdatanew,gwdatanew5,0,0+90)
@@ -588,7 +588,7 @@ print('mismatch:', mismatch)
 print('snr:', EasySNRT(timesrd_final,gwdatanew,gwdatanew,0,0+90)/error**2)
 
 
-# In[27]:
+# In[28]:
 
 
 if error_str and error_val==0:
@@ -601,7 +601,7 @@ if error_str and error_val==0:
     plt.legend()
 
 
-# In[28]:
+# In[29]:
 
 
 if parser.has_option('rd-model','phase_alignment'):
@@ -610,7 +610,7 @@ else:
     phase_alignment=False
 
 
-# In[29]:
+# In[30]:
 
 
 # Phase alignement
@@ -637,7 +637,7 @@ if phase_alignment:
         EasySNRT(timesrd_final,gwdatanew,gwdatanew5/error,0,0+90)
 
 
-# In[30]:
+# In[31]:
 
 
 #Test the new interpolated data
@@ -649,7 +649,7 @@ if error_str and error_val==0:
     plt.legend()
 
 
-# In[31]:
+# In[32]:
 
 
 #Test the error data
@@ -659,7 +659,7 @@ if error_str and error_val==0:
     plt.legend()
 
 
-# In[32]:
+# In[33]:
 
 
 #Test the error data
@@ -671,7 +671,7 @@ if error_str and error_val==0:
     plt.legend()
 
 
-# In[33]:
+# In[34]:
 
 
 #Take the piece of waveform you want
@@ -690,7 +690,7 @@ else:
     error_tsh=1
 
 
-# In[34]:
+# In[35]:
 
 
 
@@ -699,7 +699,7 @@ plt.xlabel(r'$t/M$')
 plt.ylabel(r'$r \, h_+$')
 
 
-# In[35]:
+# In[36]:
 
 
 #Fitting
@@ -799,7 +799,7 @@ dict = {'w-tau': model_dv_tau , 'w-q': model_dv_q, 'w-tau-fixed': model_dv,'w-ta
 dict_omega = {'berti': QNM_Berti , 'qnm': QNM_spectrum}
 
 
-# In[36]:
+# In[37]:
 
 
 nll = lambda *args: -log_likelihood(*args)
@@ -817,7 +817,7 @@ else:
     vars_ml=soln.x
 
 
-# In[37]:
+# In[38]:
 
 
 mypool = Pool(nbcores)
@@ -1007,39 +1007,45 @@ fg, ax = dyplot.cornerplot(res, color='blue',
 )
 
 
-# In[52]:
+# In[93]:
 
 
 if not eval(overwrite):
     fg.savefig(corner_plot, format = 'png', bbox_inches = 'tight')
     if model == 'w-tau-fixed-m-af':
+        truths=np.concatenate((w,tau))
+        fmass_spin=(samps.T)[-2:].T
+        fmass_spin_dist=[None]*len(fmass_spin)
+        labels_mf = np.concatenate((w_lab,tau_lab))
+        for i in range(len(fmass_spin)):
+            fmass_spin_dist[i]=np.concatenate(dict_omega[qnm_model](fmass_spin[i,0],fmass_spin[i,1],2,2))
+        
+        figure = corner.corner(fmass_spin_dist,truths=truths,quantiles=[0.05,0.95],labels=labels_mf,smooth=True,color='b',truth_color='r')
         figure.savefig(corner_plot_extra, format = 'png', bbox_inches = 'tight')
 
 
-# In[64]:
-
+# In[94]:
 
-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[54]:
+# In[95]:
 
 
 if not eval(overwrite):
     fig.savefig(diagnosis_plot, format = 'png', dpi = 384, bbox_inches = 'tight')
 
 
-# In[55]:
+# 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$')
-onesig_bounds = np.array([np.percentile(samps[:, i], [16, 84]) for i in range(len(samps[0]))]).T
+onesig_bounds = np.array([np.percentile(samps[:, i], [5, 95]) for i in range(len(samps[0]))]).T
 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:
@@ -1051,14 +1057,14 @@ plt.ylabel('h')
 plt.show()
 
 
-# In[56]:
+# In[99]:
 
 
 if not eval(overwrite):
     figband.savefig(fit_plot)
 
 
-# In[57]:
+# In[100]:
 
 
 if not eval(overwrite):