diff --git a/code_new/NR_dynesty_t0_loop.py b/code_new/NR_dynesty_t0_loop.py index 4936606d14fd3c32cc0d91ea471029821906dddd..da45b76800bc0a7e15884d1409c2545825544f02 100644 --- a/code_new/NR_dynesty_t0_loop.py +++ b/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):