added the cte error option

code_new/NR_dynesty_t0_loop.py

-#!/usr/bin/env python
-# coding: utf-8
-
-# In[48]:
-
-
 #Import relevant modules, import data and all that
 import numpy as np
 from scipy import interpolate
@@ -44,14 +38,10 @@ try:
     parser.sections()
 except SystemExit: 
     parser = ConfigParser()
-    parser.read('config.ini')
+    parser.read('config_n1.ini')
     parser.sections()
     pass
 
-
-# In[49]:
-
-
 # path
 rootpath=parser.get('nr-paths','rootpath')
 
@@ -66,18 +56,10 @@ overwrite = parser.get('setup','overwrite')
 downfactor = np.int(parser.get('setup','plot_down_factor'))
 sampler = parser.get('setup','sampler')
 
-
-# In[50]:
-
-
 if not os.path.exists(output_folder):
     os.mkdir(output_folder)
     print("Directory " , output_folder ,  " Created ")
 
-
-# In[51]:
-
-
 # time config
 
 tshift=parser.get('time-setup','tshift')
@@ -89,10 +71,6 @@ tend = np.float(tend)
 t_align=parser.get('time-setup','t_align')
 t_align = np.float(t_align)
 
-
-# In[52]:
-
-
 # n-tones & nlive
 
 nmax=parser.get('n-tones','nmax')
@@ -101,13 +79,17 @@ nmax = np.int(nmax)
 npoints=parser.get('n-live-points','npoints')
 npoints = np.int(npoints)
 
-
-# In[81]:
-
-
 # model
 model=parser.get('rd-model','model')
 error_str = eval(parser.get('rd-model','error_str'))
+if error_str:
+    error_val=np.float(parser.get('rd-model','error_val'))
+    if error_val==0:
+        error_type='NR_estimate'
+    else:
+        error_type=error_val
+
+
 if model == 'w-tau':
     tau_var_str='tau'
 elif model == 'w-q':
@@ -119,35 +101,26 @@ print('model:',model)
 print('nmax:',nmax)
 print('tshift:',tshift)
 print('error:', error_str)
+print('error value:',error_type)
 
-
-# In[83]:
-
-
+if error_str:
+    output_folder_1=output_folder+'/'+model+'-nmax'+str(nmax)+'_'+str(error_str)+'_'+str(error_type)
+else:
     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[84]:
-
-
 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[85]:
+samples_file=output_folder_1+'/posterior_samples-'+str(simulation_number)+'tshift_'+str(tshift)+'_'+model+'_nmax_'+str(nmax)+'.csv'
 
 
 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[86]:
-
-
 # loading priors
 w_mins=np.empty(nmax+1)
 w_maxs=np.empty(nmax+1)
@@ -193,10 +166,6 @@ if model == 'w-tau-fixed':
     priors_max = np.concatenate((a_maxs,ph_maxs))
     prior_dim = len(priors_min)
 
-
-# In[87]:
-
-
 vary_fund = True
 
 #sampler parameters
@@ -222,8 +191,8 @@ def EasyMatchT(t,h1,h2,tmin,tmax):
     #Computes the match for complex waveforms
     pos = np.argmax(t >= (tmin));
     
-    h1red=h1[pos:-1];
-    h2red=h2[pos:-1];
+    h1red=h1[pos:];
+    h2red=h2[pos:];
     
     norm1=np.sum(np.abs(h1red)**2)
     norm2=np.sum(np.abs(h2red)**2)
@@ -233,6 +202,18 @@ def EasyMatchT(t,h1,h2,tmin,tmax):
     
     return res
 
+def EasySNRT(t,h1,h2,tmin,tmax):
+    #Computes the match for complex waveforms
+    pos = np.argmax(t >= (tmin));
+    
+    h1red=h1[pos:];
+    h2red=h2[pos:];
+
+    myTable=h1red*np.conjugate(h2red)
+    res=2*np.sqrt((np.sum(myTable)).real)
+    
+    return res
+
 def wRD_to_f_Phys(f,M):
     c=2.99792458*10**8;G=6.67259*10**(-11);MS=1.9885*10**30;
     return (c**3/(M*MS*G*2*np.pi))*f
@@ -242,8 +223,16 @@ def tauRD_to_t_Phys(tau,M):
     return ((M*MS*G)/c**3)*tau
 
 
-# In[88]:
-
+def twopoint_autocovariance(t,n):
+    #Computes the match for complex waveforms
+    dt=t[1]-t[0]
+    res = np.zeros(len(n))
+    taus = np.zeros(len(n))
+    for tau in range(0,int(len(n)/2)):
+        ntau=np.roll(n, tau)
+        taus[tau] = t[tau]
+        res[tau]=np.sum(n*ntau).real
+    return (taus[:int(len(n)/2)],res[:int(len(n)/2)])
 
 #This loads the 22 mode data
 gw = {}
@@ -273,21 +262,13 @@ times5 = gw5_sxs_bbh_0305[:,0]
 tmax5=FindTmaximum(gw5_sxs_bbh_0305)
 times5 = times5 - tmax5
 
-
-# In[89]:
-
-
 #Select the data from 0 onwards
 position = np.argmax(times >= (t_align))
 position5 = np.argmax(times5 >= (t_align))
-gw_sxs_bbh_0305rd=gw_sxs_bbh_0305[position+1:-1]
-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[90]:
-
+gw_sxs_bbh_0305rd=gw_sxs_bbh_0305[position+1:]
+gw_sxs_bbh_0305rd5=gw5_sxs_bbh_0305[position5+1:]
+timesrd=gw_sxs_bbh_0305[position:-1][:,0][:]-tmax
+timesrd5=gw5_sxs_bbh_0305[position5:-1][:,0][:]-tmax5
 
 #Test plot real part (data was picked in the last cell). Aligning in time
 plt.figure(figsize = (12, 8))
@@ -297,10 +278,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[91]:
-
-
 #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$')
@@ -309,29 +286,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[92]:
-
-
 # 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[93]:
-
-
 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[94]:
-
-
 if timesrd5[-1]>= timesrd[-1]: 
     timesrd_final = timesrd
 else:
@@ -345,17 +310,36 @@ gwdatanew_im5 = gwnew_im5(timesrd_final)
 gwdatanew = gwdatanew_re - 1j*gwdatanew_im
 gwdatanew5 = gwdatanew_re5- 1j*gwdatanew_im5
 
-
-# In[95]:
-
+taus, corr= twopoint_autocovariance(timesrd,gwdatanew-gwdatanew5)
+plt.figure(figsize = (12, 8))
+plt.plot(taus, corr,'ro')
+plt.show()
+vmax=np.max(corr)
+index = np.argmax(corr == vmax)
+taus[index]
 
 mismatch=1-EasyMatchT(timesrd_final,gwdatanew,gwdatanew5,0,0+90)
 error=np.sqrt(2*mismatch)
-print(mismatch)
-
+print('error estimate:',error)
+print('mismatch:', mismatch)
+print('snr:', EasySNRT(timesrd_final,gwdatanew,gwdatanew,0,0+90)/error**2)
 
-# In[96]:
+if error_str and error_val==0:
+    error = np.sqrt(gwdatanew*gwdatanew-2*gwdatanew*gwdatanew5+gwdatanew5*gwdatanew5)
+    error_est=np.sqrt(error.imag**2+error.real**2)
+    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')
+    plt.plot(timesrd_final, error_est, "b", alpha=0.3, lw=2, label='error')
+    plt.legend()
 
+if error_str and error_val==0:
+    error = np.sqrt(gwdatanew*gwdatanew-2*gwdatanew*gwdatanew5+gwdatanew5*gwdatanew5)
+    error_est=np.sqrt(error.imag**2+error.real**2)
+    plt.figure(figsize = (12, 8))
+    plt.xlim(0,80)
+    plt.plot(timesrd_final, np.abs(gwdatanew)/error_est, "r", alpha=0.3, lw=2, label='Lev6')
+    plt.legend()
 
 # Phase alignement
 phas = np.angle(gwdatanew)
@@ -365,10 +349,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[97]:
-
-
 position = np.argmax(timesrd_final >= (t_align))
 dphase = phas5[position]-phas[position]
 print(dphase)
@@ -384,44 +364,40 @@ 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[98]:
-
-
 mismatch=1-EasyMatchT(timesrd_final,gwdatanew,gwdatanew5,0,+90)
-print(mismatch)
+error=np.sqrt(2*mismatch)
+print('error estimate:',error)
+print('mismatch:', mismatch)
+print('snr:', EasySNRT(timesrd_final,gwdatanew,gwdatanew5,0,0+90)/error)
 if error_str:
     error = np.sqrt(gwdatanew*gwdatanew-2*gwdatanew*gwdatanew5+gwdatanew5*gwdatanew5)
+    error_est=np.sqrt(error.imag**2+error.real**2)
 else :
     error = 1 
 
-
-# In[100]:
-
+EasySNRT(timesrd_final,gwdatanew,gwdatanew5/error,0,0+90)
 
 #Test the new interpolated data
-if error_str:
+if error_str and error_val==0:
     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')
     plt.plot(timesrd_final, error.real, "b", alpha=0.3, lw=2, label='error')
     plt.legend()
 
-
-# In[101]:
-
-
 #Test the error data
-if error_str:
+if error_str and error_val==0:
     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')
     plt.plot(timesrd_final, np.sqrt(error.imag**2+error.real**2), "r", alpha=0.3, lw=2, label='all error')
     plt.legend()
 
-
-# In[103]:
-
+#Test the error data
+if error_str and error_val==0:
+    plt.figure(figsize = (12, 8))
+    plt.xlim(1,40)
+    plt.ylim(-300,300)
+    plt.plot(timesrd_final,gwdatanew.real/np.sqrt(error.imag**2+error.real**2), "r", alpha=0.3, lw=2, label='all error')
+    plt.legend()
 
 #Take the piece of waveform you want
 position_in = np.argmax(timesrd_final >= tshift)
@@ -429,14 +405,10 @@ 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:
+if error_str and error_val==0:
     error_tsh=error[position_in:position_end]
 else:
-    error_tsh=1
-
-
-# In[104]:
-
+    error_tsh=error_val
 
 #Fitting
 #RD model for nmax tones. Amplitudes are in (xn*Exp[i yn]) version. Used here.
@@ -509,16 +481,8 @@ def log_probability(theta):
         return -np.inf
     return lp + log_likelihood(theta)
 
-
-# In[105]:
-
-
 dict = {'w-tau': model_dv_tau , 'w-q': model_dv_q, 'w-tau-fixed': model_dv}
 
-
-# In[106]:
-
-
 #I need to provid an initial guess for 4*(nmax+1) the parameters
 np.random.seed(42)
 nll = lambda *args: -log_likelihood(*args)
@@ -529,16 +493,8 @@ print("Maximum likelihood estimates:")
 vars_ml=soln.x
 print(vars_ml)
 
-
-# In[107]:
-
-
 f2=dynesty.NestedSampler(log_likelihood, prior_transform, prior_dim, nlive=npoints,sample=sampler)
-f2.run_nested()
-
-
-# In[108]:
-
+f2.run_nested(dlogz=0.01)
 
 wstr = r'$\omega_'
 
@@ -568,10 +524,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[109]:
-
-
 if model=='w-tau-fixed':
     rg = (nmax+1)
 else:
@@ -591,32 +543,17 @@ else:
         amps_aux = samps_tr[i][half_points:-1]
         npamps[i] = np.quantile(amps_aux, 0.5)
 
-
-# In[110]:
-
-
 res = f2.results
 res.samples_u.shape
 res.summary()
 samps=f2.results.samples
 
-
-# In[111]:
-
-
 evidence = res.logz[-1]
 evidence_error = res.logzerr[-1]
 
-
-# In[112]:
-
-
 summary_titles=['n','id','t_shift','dlogz','dlogz_err']
 
-
-# In[113]:
-
-
+if not eval(overwrite):
     if os.path.exists(sumary_data):
         outvalues = np.array([[nmax, simulation_number, tshift, evidence,evidence_error]])
     else:
@@ -629,49 +566,30 @@ with open(sumary_data, 'a') as file:
         else:
             writer.writerow(outvalues[0])
 
-
-# In[114]:
-
-
 samps=f2.results.samples
 samps_tr=np.transpose(samps)
 
-
-# In[115]:
-
-
 sigma_vars_m = np.empty(prior_dim)
 sigma_vars_p = np.empty(prior_dim)
 sigma_vars = np.empty(prior_dim)
 for i in range(prior_dim): 
-    amps_aux = samps_tr[i][half_points:-1]
+    amps_aux = samps_tr[i][half_points:]
     sigma_vars_m[i] = np.quantile(amps_aux, 0.1)
     sigma_vars[i] = np.quantile(amps_aux, 0.5)
     sigma_vars_p[i] = np.quantile(amps_aux, 0.9)
 
-
-# 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[117]:
-
-
 key =['max val','lower bound','higher bound']
 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 not eval(overwrite):
     if os.path.exists(best_data):
         df2.to_csv(best_data, mode='a', header=False,index = True)
     else:
         df2.to_csv(best_data,index = True)
 
-
-# In[118]:
-
-
 if model == 'w-q':
     tau_val = np.pi*w*tau
     truths = np.concatenate((w,tau_val,npamps))
@@ -681,10 +599,6 @@ elif model == 'w-tau':
 elif model == 'w-tau-fixed':
     truths = npamps
 
-
-# In[119]:
-
-
 fg, ax = dyplot.cornerplot(res, color='blue', 
                            show_titles=True,
                            labels=labels,
@@ -693,32 +607,18 @@ fg, ax = dyplot.cornerplot(res, color='blue',
                            truth_color='red',
 )
 
-
-# In[121]:
-
-
+if not eval(overwrite):
     fg.savefig(corner_plot, format = 'png', bbox_inches = 'tight')
 
-
-# In[122]:
-
-
 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[123]:
-
-
+if not eval(overwrite):
     fig.savefig(diagnosis_plot, format = 'png', dpi = 384, bbox_inches = 'tight')
 
-
-# In[124]:
-
-
 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$')
@@ -726,22 +626,18 @@ 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.04, lw=3)
+    plt.plot(timesrd_final_tsh, dict[model](sample).real, "r-", alpha=0.1, lw=1)
 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[125]:
-
-
+if not eval(overwrite):
     figband.savefig(fit_plot)
 
-
-# In[ ]:
-
-
-
-
+if not eval(overwrite):
+    with open(samples_file,'w') as file:
+        writer = csv.writer(file)
+        writer.writerow(labels)
+        writer.writerows(samps[::downfactor])