diff --git a/code_new/NR_dynesty_t0_loop.py b/code_new/NR_dynesty_t0_loop.py
index 83a5aab76b3ea3dd56ca13951abae244395462f8..4936606d14fd3c32cc0d91ea471029821906dddd 100644
--- a/code_new/NR_dynesty_t0_loop.py
+++ b/code_new/NR_dynesty_t0_loop.py
@@ -1,10 +1,42 @@
#!/usr/bin/env python
# coding: utf-8
-# In[286]:
+# In[1]:
+
+
+# Copyright (C) 2020 Xisco Jimenez Forteza
+#
+# This program is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 3 of the License, or (at your
+# option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
+# Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with this program; if not, write to the Free Software Foundation, Inc.,
+# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+
+
+#
+# =============================================================================
+#
+# Preamble
+#
+# =============================================================================
+#
+"""Generate ringdown templates in the time and perform parameter estimation on them.
+"""
+
+
+# In[2]:
#Import relevant modules, import data and all that
+import time
import numpy as np
from scipy import interpolate
import corner
@@ -17,11 +49,8 @@ import codecs
#rc('text', usetex=True)
plt.rcParams.update({'font.size': 16.5})
-import ptemcee
-#from pycbc.pool import choose_pool
from multiprocessing import Pool
import h5py
-import inspect
import pandas as pd
import json
import qnm
@@ -39,7 +68,10 @@ from scipy.optimize import fsolve
from scipy.optimize import least_squares
+# In[3]:
+
+# Cell that calls the arguments from your 'config.ini' file.
try:
parser = argparse.ArgumentParser(description="Simple argument parser")
parser.add_argument("-c", action="store", dest="config_file")
@@ -50,12 +82,12 @@ try:
parser.sections()
except SystemExit:
parser = ConfigParser()
- parser.read('config_fixed_n1_m_af.ini')
+ parser.read('config_n2_q10.ini')
parser.sections()
pass
-# In[287]:
+# In[4]:
# path
@@ -73,8 +105,13 @@ downfactor = np.int(parser.get('setup','plot_down_factor'))
sampler = parser.get('setup','sampler')
nr_code = parser.get('setup','nr_code')
+if parser.has_option('setup','export'):
+ export=eval(parser.get('setup','export'))
+else:
+ export=True
+
-# In[288]:
+# In[5]:
if parser.has_option('setup','nb_cores'):
@@ -83,7 +120,7 @@ else:
nbcores = 1
-# In[289]:
+# In[6]:
if not os.path.exists(output_folder):
@@ -91,7 +128,7 @@ if not os.path.exists(output_folder):
print("Directory " , output_folder , " Created ")
-# In[290]:
+# In[7]:
# time config
@@ -105,7 +142,7 @@ t_align=parser.get('time-setup','t_align')
t_align = np.float(t_align)
-# In[291]:
+# In[8]:
# n-tones & nlive
@@ -117,7 +154,7 @@ npoints=parser.get('n-live-points','npoints')
npoints = np.int(npoints)
-# In[292]:
+# In[9]:
# model
@@ -146,7 +183,7 @@ print('error:', error_str)
print('error value:',error_type)
-# In[293]:
+# In[10]:
if error_str:
@@ -159,7 +196,7 @@ if not os.path.exists(output_folder_1):
print("Directory " , output_folder_1 , " Created ")
-# In[294]:
+# In[11]:
corner_plot=output_folder_1+'/Dynesty_'+str(simulation_number)+'_'+model+'_nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'corner_plot.png'
@@ -169,14 +206,14 @@ fit_plot=output_folder_1+'/Fit_results_'+str(simulation_number)+'tshift_'+str(ts
samples_file=output_folder_1+'/posterior_samples-'+str(simulation_number)+'tshift_'+str(tshift)+'_'+model+'_nmax_'+str(nmax)+'.csv'
-# In[295]:
+# In[12]:
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[296]:
+# In[13]:
vary_fund = True
@@ -192,7 +229,8 @@ mediancolor = '#f7695c' #'#9b2814'
#Import data and necessary functions
#TimeOfMaximum
def FindTmaximum(y):
- #Determines the maximum absolute value of the complex waveform
+ """ It determines the maximum absolute value of the complex waveform.
+ """
absval = np.sqrt(y[:,1]*y[:,1]+y[:,2]*y[:,2])
vmax=np.max(absval)
index = np.argmax(absval == vmax)
@@ -201,7 +239,8 @@ def FindTmaximum(y):
def EasyMatchT(t,h1,h2,tmin,tmax):
- #Computes the match for complex waveforms
+ """ It computes the time-domain match for (h1|h2) complex waveforms.
+ """
pos = np.argmax(t >= (tmin));
h1red=h1[pos:];
@@ -216,7 +255,8 @@ def EasyMatchT(t,h1,h2,tmin,tmax):
return res
def EasySNRT(t,h1,h2,tmin,tmax):
- #Computes the match for complex waveforms
+ """ It computes the time-domain snr for (h1|h2) complex waveforms.
+ """
pos = np.argmax(t >= (tmin));
h1red=h1[pos:];
@@ -228,16 +268,32 @@ def EasySNRT(t,h1,h2,tmin,tmax):
return res
def wRD_to_f_Phys(f,M):
+ """ It converts NR frequencies to physical units in [Hz].
+ """
+ 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
+
+def wRD_to_f_NR(f,M):
+ """ It converts Physical frequencies to NR units in [1/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
def tauRD_to_t_Phys(tau,M):
+ """ It converts NR times to physical units in [s].
+ """
c=2.99792458*10**8;G=6.67259*10**(-11);MS=1.9885*10**30;
return ((M*MS*G)/c**3)*tau
+def tauRD_to_t_NR(tau,M):
+ """ It converts Physical times to NR units in [M].
+ """
+ c=2.99792458*10**8;G=6.67259*10**(-11);MS=1.9885*10**30;
+ return 1/((M*MS*G)/c**3)*tau
def twopoint_autocovariance(t,n):
- #Computes the match for complex waveforms
+ """ It computes the two-point autocovariance function.
+ """
dt=t[1]-t[0]
res = np.zeros(len(n))
taus = np.zeros(len(n))
@@ -247,38 +303,71 @@ def twopoint_autocovariance(t,n):
res[tau]=np.sum(n*ntau).real
return (taus[:int(len(n)/2)],res[:int(len(n)/2)])
-def QNM_spectrum(mf,af,l,m):
+grav_220 = [qnm.modes_cache(s=-2,l=2,m=2,n=i) for i in range (0,dim)]
+def QNM_spectrum_old(mf,af,l,m):
+ """ It computes the RD frequencies and damping times in NR units.
+ """
+ omegas_new_v2 = [qnm.modes_cache(s=-2,l=l,m=m,n=i)(a=af)[0] for i in range (0,dim)]
+ w_m_a = (np.real(omegas_new_v2))/mf
+ tau_m_a=-1/(np.imag(omegas_new_v2))*mf
- omegas_new = [qnm.modes_cache(s=-2,l=l,m=m,n=i)(a=af)[0] for i in range (0,dim)]
+ return (w_m_a, tau_m_a)
+
+def QNM_spectrum(mf,af,l,m):
+ """ It computes the RD frequencies and damping times in NR units.
+ """
+ omegas_new=np.asarray([grav_220[i](a=af)[0] for i in range (0,dim)])
w_m_a = (np.real(omegas_new))/mf
tau_m_a=-1/(np.imag(omegas_new))*mf
return (w_m_a, tau_m_a)
+def QNM_Berti(mf,af,l,m):
+ """ It computes the RD frequencies and damping times in NR units.
+ """
+ position=np.argmax(rdowndata[0,0] >= (af))
+ #w_m_a=f1+f2*(1-af)**f3
+ w_m_a=[None]*(nmax+1)
+ tau_ma_a=[None]*(nmax+1)
+
+ for i in range(nmax+1):
+ qnm=rdowndata[0,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]:
+
+
+np.sqrt(12/2*1/tauRD_to_t_NR(1.3*10**-47,70)*(5*10**(-21))**2)
-# In[297]:
+
+# In[15]:
+
+
+np.sqrt(0.004/2*1/(1.3*10**-47)*(5*10**(-21))**2)
+
+
+# In[16]:
gw = {}
-gw[simulation_number] = h5py.File(simulation_path_1, 'r')
-if nr_code=='SXS':
- gw_sxs_bbh_0305 = gw[simulation_number]["Extrapolated_N3.dir"]["Y_l2_m2.dat"]
+gw[simulation_number] = h5py.File(simulation_path_1, 'r')
+if nr_code=='SXS':
+ gw_sxs_bbh_0305 = gw[simulation_number]["Extrapolated_N3.dir"]["Y_l2_m2.dat"]
times = gw_sxs_bbh_0305[:,0]
-
gw5 = {}
gw5[simulation_number] = h5py.File(simulation_path_2, 'r')
gw5_sxs_bbh_0305 = gw5[simulation_number]["Extrapolated_N3.dir"]["Y_l2_m2.dat"]
-# Remember to download metadata.json from the simulation with number: 0305. Download Lev6/metadata.json
-# This postprocesses the metadata file to find the final mass and final spin
-
-elif nr_code=='Maya':
- gw_sxs_bbh_0305_amp = np.asarray(gw[simulation_number]['amp_l2_m2/Y'])[6:]
- gw_sxs_bbh_0305_amp_err = np.asarray(gw[simulation_number]['amp_l2_m2/errors'])
- times = np.asarray(gw[simulation_number]['amp_l2_m2/X'])[6:]
- gw_sxs_bbh_0305_amp_int = interp1d(times, -gw_sxs_bbh_0305_amp, kind='cubic')
- gw_sxs_bbh_0305_amp_errs_int = interp1d(times, gw_sxs_bbh_0305_amp_err, kind='cubic')
-
+elif nr_code=='Maya':
+ gw_sxs_bbh_0305_amp = np.asarray(gw[simulation_number]['amp_l2_m2/Y'])[6:]
+ gw_sxs_bbh_0305_amp_err = np.asarray(gw[simulation_number]['amp_l2_m2/errors'])
+ times = np.asarray(gw[simulation_number]['amp_l2_m2/X'])[6:]
+ gw_sxs_bbh_0305_amp_int = interp1d(times, -gw_sxs_bbh_0305_amp, kind='cubic')
+ gw_sxs_bbh_0305_amp_errs_int = interp1d(times, gw_sxs_bbh_0305_amp_err, kind='cubic')
gw_sxs_bbh_0305_pha = np.asarray(gw[simulation_number]['phase_l2_m2/Y'])[6:]
gw_sxs_bbh_0305_pha_err = np.asarray(gw[simulation_number]['phase_l2_m2/errors'])
times = np.asarray(gw[simulation_number]['phase_l2_m2/X'])[6:]
@@ -290,15 +379,14 @@ elif nr_code=='Maya':
phs=gw_sxs_bbh_0305_pha_int(times)
phs_err=-gw_sxs_bbh_0305_pha_errs_int(times)
-
+
gw_sxs_bbh_0305 = np.asarray([times,amps*np.cos(phs),amps*np.sin(phs)]).T
gw5_sxs_bbh_0305 = np.asarray([times,(amps+amps_err)*np.cos(phs+phs_err),amps*np.sin(phs+phs_err)]).T
-
-elif nr_code=='LaZeV':
- gw_sxs_bbh_0305_amp = np.asarray(gw[simulation_number]['amp_l2_m2/Y'])[6:]
- gw_sxs_bbh_0305_amp_err = np.asarray(gw[simulation_number]['amp_l2_m2/errors'])
- times_1 = np.asarray(gw[simulation_number]['amp_l2_m2/X'])[6:]
- gw_sxs_bbh_0305_amp_int = interp1d(times_1, -gw_sxs_bbh_0305_amp, kind='cubic')
+elif nr_code=='LaZeV':
+ gw_sxs_bbh_0305_amp = np.asarray(gw[simulation_number]['amp_l2_m2/Y'])[6:]
+ gw_sxs_bbh_0305_amp_err = np.asarray(gw[simulation_number]['amp_l2_m2/errors'])
+ times_1 = np.asarray(gw[simulation_number]['amp_l2_m2/X'])[6:]
+ gw_sxs_bbh_0305_amp_int = interp1d(times_1, -gw_sxs_bbh_0305_amp, kind='cubic')
gw_sxs_bbh_0305_amp_errs_int = interp1d(times_1, gw_sxs_bbh_0305_amp_err, kind='cubic')
gw_sxs_bbh_0305_pha = np.asarray(gw[simulation_number]['phase_l2_m2/Y'])[6:]
@@ -306,19 +394,19 @@ elif nr_code=='LaZeV':
times = np.asarray(gw[simulation_number]['phase_l2_m2/X'])[6:]
gw_sxs_bbh_0305_pha_int = interp1d(times, -gw_sxs_bbh_0305_pha, kind='cubic')
gw_sxs_bbh_0305_pha_errs_int = interp1d(times, gw_sxs_bbh_0305_pha_err, kind='cubic')
-
+
amps=gw_sxs_bbh_0305_amp_int(times_1)
amps_err=gw_sxs_bbh_0305_amp_errs_int(times_1)
phs=gw_sxs_bbh_0305_pha_int(times_1)
phs_err=-gw_sxs_bbh_0305_pha_errs_int(times_1)
-
+
gw_sxs_bbh_0305 = np.asarray([times_1,amps*np.cos(phs),amps*np.sin(phs)]).T
gw5_sxs_bbh_0305 = np.asarray([times_1,(amps+amps_err)*np.cos(phs+phs_err),amps*np.sin(phs+phs_err)]).T
times=times_1
-# In[298]:
+# In[17]:
if nr_code=='SXS':
@@ -342,20 +430,26 @@ times = gw_sxs_bbh_0305[:,0]
tmax=FindTmaximum(gw_sxs_bbh_0305[round(len(gw_sxs_bbh_0305)/2):])
times = times - tmax
-
#times 6--> x axis of your data
times5 = gw5_sxs_bbh_0305[:,0]
tmax5=FindTmaximum(gw5_sxs_bbh_0305[round(len(gw_sxs_bbh_0305)/2):])
times5 = times5 - tmax5
-# In[299]:
+# In[18]:
-w , tau = QNM_spectrum(mf,af,2,2)
+if parser.has_option('setup','qnm_model'):
+ qnm_model='berti'
+ rdownfolders=np.asarray([rootpath+'/RDmodes/l2/n'+str(i+1)+'l2m2.dat' for i in range(nmax+1)])
+ rdowndata = np.asarray([np.loadtxt(rdownfolders[i]).T for i in range(len(rdownfolders))])
+ w , tau = QNM_Berti(mf,af,2,2)
+else:
+ qnm_model='qnm'
+ w , tau = QNM_spectrum(mf,af,2,2)
-# In[300]:
+# In[19]:
# loading priors
@@ -415,7 +509,7 @@ elif model == 'w-tau-fixed-m-af':
priors=np.column_stack((priors_min,priors_max))
-# In[301]:
+# In[20]:
#Select the data from 0 onwards
@@ -427,7 +521,7 @@ timesrd=gw_sxs_bbh_0305[position:-1][:,0][:]-tmax
timesrd5=gw5_sxs_bbh_0305[position5:-1][:,0][:]-tmax5
-# In[302]:
+# In[21]:
#Test plot real part (data was picked in the last cell). Aligning in time
@@ -439,7 +533,13 @@ plt.plot(timesrd5, np.sqrt(gw_sxs_bbh_0305rd5[:,1]**2+gw_sxs_bbh_0305rd5[:,2]**2
plt.legend()
-# In[303]:
+# In[22]:
+
+
+#[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]:
gwnew_re = interpolate.interp1d(timesrd, gw_sxs_bbh_0305rd[:,1], kind = 'cubic')
@@ -449,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[304]:
+# In[24]:
if timesrd5[-1]>= timesrd[-1]:
@@ -466,7 +566,7 @@ gwdatanew = gwdatanew_re - 1j*gwdatanew_im
gwdatanew5 = gwdatanew_re5- 1j*gwdatanew_im5
-# In[305]:
+# In[25]:
#taus, corr= twopoint_autocovariance(timesrd,gwdatanew-gwdatanew5)
@@ -478,7 +578,7 @@ gwdatanew5 = gwdatanew_re5- 1j*gwdatanew_im5
#taus[index]
-# In[306]:
+# In[26]:
mismatch=1-EasyMatchT(timesrd_final,gwdatanew,gwdatanew5,0,0+90)
@@ -488,7 +588,7 @@ print('mismatch:', mismatch)
print('snr:', EasySNRT(timesrd_final,gwdatanew,gwdatanew,0,0+90)/error**2)
-# In[307]:
+# In[27]:
if error_str and error_val==0:
@@ -501,7 +601,7 @@ if error_str and error_val==0:
plt.legend()
-# In[308]:
+# In[28]:
if parser.has_option('rd-model','phase_alignment'):
@@ -510,7 +610,7 @@ else:
phase_alignment=False
-# In[309]:
+# In[29]:
# Phase alignement
@@ -522,7 +622,7 @@ if phase_alignment:
position = np.argmax(timesrd_final >= (t_align))
dphase = phas5[position]-phas[position]
gwdatanew = (gwdatanew_re - 1j*gwdatanew_im)*np.exp(1j*dphase)
- gw_sxs_bbh_0305rd6=gw6_sxs_bbh_0305[position6:-1]
+ gw_sxs_bbh_0305rd6=gw5_sxs_bbh_0305[position5:-1]
timesrd=gw_sxs_bbh_0305[position:-1][:,0][:920]
mismatch=1-EasyMatchT(timesrd_final,gwdatanew,gwdatanew5,0,+90)
error=np.sqrt(2*mismatch)
@@ -534,11 +634,10 @@ if phase_alignment:
error_est=np.sqrt(error.imag**2+error.real**2)
else :
error = 1
+ EasySNRT(timesrd_final,gwdatanew,gwdatanew5/error,0,0+90)
- EasySNRT(timesrd_final,gwdatanew,gwdatanew5/error,0,0+90)
-
-# In[310]:
+# In[30]:
#Test the new interpolated data
@@ -550,7 +649,7 @@ if error_str and error_val==0:
plt.legend()
-# In[311]:
+# In[31]:
#Test the error data
@@ -560,7 +659,7 @@ if error_str and error_val==0:
plt.legend()
-# In[312]:
+# In[32]:
#Test the error data
@@ -572,11 +671,10 @@ if error_str and error_val==0:
plt.legend()
-# In[313]:
+# In[33]:
#Take the piece of waveform you want
-tshift=0
position_in = np.argmax(timesrd_final >= tshift)
position_end = np.argmax(timesrd_final >= tend)
timesrd_final_tsh = timesrd_final[position_in:position_end]
@@ -592,14 +690,23 @@ else:
error_tsh=1
-# In[314]:
+# In[34]:
+
+
+
+plt.errorbar(timesrd_final_tsh[::10], gwdatanew_re_tsh[::10], yerr=10*np.sqrt((error_tsh.real**2+error_tsh.imag**2))[::10],fmt='o',color='red')
+plt.xlabel(r'$t/M$')
+plt.ylabel(r'$r \, h_+$')
+
+
+# In[35]:
#Fitting
#RD model for nmax tones. Amplitudes are in (xn*Exp[i yn]) version. Used here.
def model_dv_q(theta):
- #x0, y0= theta
- #Your nmax might not align with the dim of theta. Better check it here.
+ """RD model parametrized with the quality factor q.
+ """
assert int(len(theta)/4) == dim, 'Please recheck your n and parameters'
wvars = theta[ : (dim)]
@@ -614,8 +721,8 @@ def model_dv_q(theta):
return ansatz
def model_dv_tau(theta):
- #x0, y0= theta
- #Your nmax might not align with the dim of theta. Better check it here.
+ """RD model parametrized with the damping time tau.
+ """
assert int(len(theta)/4) == dim, 'Please recheck your n and parameters'
wvars = theta[ : (dim)]
@@ -630,8 +737,8 @@ def model_dv_tau(theta):
return ansatz
def model_dv(theta):
- #x0, y0= theta
- #Your nmax might not align with the dim of theta. Better check it here.
+ """RD model parametrized with the damping time tau and with the QNM spectrum fixd to GR.
+ """
xvars = theta[ : (dim)]
yvars = theta[(dim) : 2*(dim)]
@@ -642,14 +749,14 @@ def model_dv(theta):
return ansatz
def model_dv_ma(theta):
- #x0, y0= theta
- #Your nmax might not align with the dim of theta. Better check it here.
+ """RD model parametrized with the damping time tau and with the QNM spectrum fixd to GR. The QNM spectrum is given from the mass and spin.
+ """
xvars = theta[ : (dim)]
yvars = theta[(dim) : 2*(dim)]
- mass_vars = np.float(theta[-2])
- spin_vars = np.float(theta[-1])
+ mass_vars = theta[-2]
+ spin_vars = theta[-1]
- w_m_a , tau_m_a = QNM_spectrum(mass_vars,spin_vars,2,2)
+ w_m_a , tau_m_a = dict_omega[qnm_model](mass_vars,spin_vars,2,2)
ansatz = 0
for i in range (0,dim):
@@ -661,12 +768,16 @@ def model_dv_ma(theta):
#It works for the (xn*Exp[iyn]) version.
def prior_transform(cube):
+ """RD uniform priors. The values for priors_min and priors_max must be given out of this function.
+ """
for i in range(prior_dim):
cube[i] = priors_min[i]+ cube[i]*(priors_max[i]-priors_min[i])
return cube
# LogLikelihood function. It is just a Gaussian loglikelihood based on computing the residuals^2
def log_likelihood(theta):
+ """chi2 likelihood.
+ """
modelev = dict[model](theta)
result = -np.sum(((gwdatanew_re_tsh - modelev.real)**2+(gwdatanew_im_tsh - modelev.imag)**2)/(2*(error_tsh.real**2+error_tsh.imag**2)))
if np.isnan(result):
@@ -677,95 +788,18 @@ def log_likelihood(theta):
# Logposterior distribution for the residuals case.
# The evidence is just a normalization factor
def log_probability(theta):
+ """Posterior likelihood.
+ """
lp = log_prior(theta)
if not np.isfinite(lp):
return -np.inf
return lp + log_likelihood(theta)
-
-# In[315]:
-
-
dict = {'w-tau': model_dv_tau , 'w-q': model_dv_q, 'w-tau-fixed': model_dv,'w-tau-fixed-m-af': model_dv_ma}
+dict_omega = {'berti': QNM_Berti , 'qnm': QNM_spectrum}
-# In[316]:
-
-
-#nmax=2
-#dim = nmax+1
-#ndim = 4*dim
-
-# loading priors
-#w_mins=np.empty(nmax+1)
-#w_maxs=np.empty(nmax+1)
-#tau_mins=np.empty(nmax+1)
-#tau_maxs=np.empty(nmax+1)
-#a_mins=np.empty(nmax+1)
-#a_maxs=np.empty(nmax+1)
-#ph_mins=np.empty(nmax+1)
-#ph_maxs=np.empty(nmax+1)
-
-
-#for i in range(nmax+1):
-# wp_min=parser.get('prior-w'+str(i),'w'+str(i)+'_min')
-# w_mins[i] = np.float(wp_min)
-
-# wp_max=parser.get('prior-w'+str(i),'w'+str(i)+'_max')
-# w_maxs[i] = np.float(wp_max)
-
-# taup_min=parser.get('prior-'+tau_var_str+str(i),tau_var_str+str(i)+'_min')
-# tau_mins[i] = np.float(taup_min)
-
-# taup_max=parser.get('prior-'+tau_var_str+str(i),tau_var_str+str(i)+'_max')
-# tau_maxs[i] = np.float(taup_max)
-
-# amp0_min=parser.get('prior-amp'+str(i),'amp'+str(i)+'_min')
-# a_mins[i] = np.float(amp0_min)
-
-# amp1_max=parser.get('prior-amp'+str(i),'amp'+str(i)+'_max')
-# a_maxs[i] = np.float(amp1_max)
-
-# phase_min=parser.get('prior-phase'+str(i),'phase'+str(i)+'_min')
-# ph_mins[i] = np.float(phase_min)*2*np.pi
-
-# phase_max=parser.get('prior-phase'+str(i),'phase'+str(i)+'_max')
-# ph_maxs[i] = np.float(phase_max)*2*np.pi
-
-#priors_min = np.concatenate((w_mins,tau_mins,a_mins,ph_mins))
-#priors_max = np.concatenate((w_maxs,tau_maxs,a_maxs,ph_maxs))
-#prior_dim = len(priors_min)
-
-#nll = lambda *args: -log_likelihood(*args)
-#initial = [w[0],w[1],w[2],tau[0],tau[1],tau[2],a_maxs[0]/2,a_maxs[1]/2,a_maxs[2]/2,1,1,1]
-#bic2=np.ones(len(range(20)))
-
-#for tshift in range(20):
-# position_in = np.argmax(timesrd_final >= tshift)
-# 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]
-# error_tsh=error[position_in:position_end]
-
-# soln = minimize(nll, initial)
-# vars_ml=soln.x
-# initial = vars_ml
-# bic2[tshift]=2*log_likelihood(vars_ml)+len(timesrd_final_tsh)*np.log(len(timesrd_final_tsh))
-
-
-# In[317]:
-
-
-#plt.figure(figsize = (12, 8))
-#plt.plot(np.abs(bic4), "r", alpha=0.3, lw=3, label=r'$bic: 4$')
-#plt.plot(np.abs(bic3), "b", alpha=0.3, lw=3, label=r'$bic: 3$')
-#plt.plot(np.abs(bic2), "cyan", alpha=0.3, lw=3, label=r'$bic: 2$')
-#plt.yscale('log')
-#plt.legend()
-
-
-# In[322]:
+# In[36]:
nll = lambda *args: -log_likelihood(*args)
@@ -783,12 +817,13 @@ else:
vars_ml=soln.x
-# In[323]:
+# In[37]:
mypool = Pool(nbcores)
mypool.size = nbcores
+start = time.process_time()
f2=dynesty.NestedSampler(log_likelihood, prior_transform, prior_dim, nlive=npoints,sample=sampler,pool=mypool)
if parser.has_option('setup','dlogz'):
dlogz=np.float(parser.get('setup','dlogz'))
@@ -796,50 +831,21 @@ if parser.has_option('setup','dlogz'):
else:
f2.run_nested()
+print(time.process_time() - start)
-# In[325]:
-
-
-wstr = r'$\omega_'
-
-if model == 'w-tau':
- taustr = r'$\tau_'
-elif model == 'w-q':
- taustr = r'$q_'
-elif model == 'w-tau-fixed':
- taustr = r'$dumb_var}'
-elif model == 'w-tau-fixed-m-af':
- taustr = r'$\tau_'
-
-ampstr = r'$A_'
-phasestr = r'$\phi_'
-
-w_lab = [None] * dim
-tau_lab = [None] * dim
-amp_lab = [None] * dim
-pha_lab = [None] * dim
-mass_lab = ['mass']
-spin_lab = ['spin']
-
-for i in range(dim):
- w_lab[i] = wstr+str(i)+'$'
- tau_lab[i] = taustr+str(i)+'$'
- amp_lab[i] = ampstr+str(i)+'$'
- pha_lab[i] = phasestr+str(i)+'$'
+# In[39]:
-labels = np.concatenate((w_lab,tau_lab,amp_lab,pha_lab))
-if model=='w-tau-fixed':
- labels = np.concatenate((amp_lab,pha_lab))
-
-if model=='w-tau-fixed-m-af':
- pha_lab[i] = phasestr+str(i)+'$'
-
- labels = np.concatenate((amp_lab,pha_lab,mass_lab,spin_lab))
+res = f2.results
+res.samples_u.shape
+res.summary()
+samps=f2.results.samples
+samps_tr=np.transpose(samps)
+half_points=int(round((len(samps_tr[0])/1.25)))
-# In[337]:
+# In[40]:
if model=='w-tau-fixed':
@@ -849,9 +855,6 @@ elif model=='w-tau-fixed':
else:
rg = (nmax+1)*2
-samps=f2.results.samples
-samps_tr=np.transpose(samps)
-half_points=int(round((len(samps_tr[0])/1.25)))
if model=='w-tau-fixed-a-mf':
npamps = np.empty((nmax+1))
@@ -865,29 +868,20 @@ else :
npamps[i] = np.quantile(amps_aux, 0.5)
-# In[338]:
-
-
-res = f2.results
-res.samples_u.shape
-res.summary()
-samps=f2.results.samples
-
-
-# In[339]:
+# In[41]:
evidence = res.logz[-1]
evidence_error = res.logzerr[-1]
-# In[340]:
+# In[42]:
summary_titles=['n','id','t_shift','dlogz','dlogz_err']
-# In[341]:
+# In[43]:
if not eval(overwrite):
@@ -904,14 +898,14 @@ if not eval(overwrite):
writer.writerow(outvalues[0])
-# In[342]:
+# In[44]:
samps=f2.results.samples
samps_tr=np.transpose(samps)
-# In[343]:
+# In[45]:
sigma_vars_m = np.empty(prior_dim)
@@ -924,14 +918,56 @@ for i in range(prior_dim):
sigma_vars_p[i] = np.quantile(amps_aux, 0.9)
-# In[344]:
+# In[46]:
+
+
+wstr = r'$\omega_'
+
+if model == 'w-tau':
+ taustr = r'$\tau_'
+elif model == 'w-q':
+ taustr = r'$q_'
+elif model == 'w-tau-fixed':
+ taustr = r'$dumb_var}'
+elif model == 'w-tau-fixed-m-af':
+ taustr = r'$\tau_'
+
+ampstr = r'$A_'
+phasestr = r'$\phi_'
+
+w_lab = [None] * dim
+tau_lab = [None] * dim
+amp_lab = [None] * dim
+pha_lab = [None] * dim
+mass_lab = ['mass']
+spin_lab = ['spin']
+
+for i in range(dim):
+ w_lab[i] = wstr+str(i)+'$'
+ tau_lab[i] = taustr+str(i)+'$'
+ amp_lab[i] = ampstr+str(i)+'$'
+ pha_lab[i] = phasestr+str(i)+'$'
+
+
+labels = np.concatenate((w_lab,tau_lab,amp_lab,pha_lab))
+
+if model=='w-tau-fixed':
+ labels = np.concatenate((amp_lab,pha_lab))
+
+if model=='w-tau-fixed-m-af':
+ pha_lab[i] = phasestr+str(i)+'$'
+
+ labels = np.concatenate((amp_lab,pha_lab,mass_lab,spin_lab))
+
+
+# In[47]:
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[345]:
+# In[48]:
key =['max val','lower bound','higher bound']
@@ -944,7 +980,7 @@ if not eval(overwrite):
df2.to_csv(best_data,index = True)
-# In[346]:
+# In[49]:
if model == 'w-q':
@@ -959,7 +995,7 @@ elif model == 'w-tau-fixed-m-af':
truths = np.concatenate((npamps,[mf],[af]))
-# In[347]:
+# In[50]:
fg, ax = dyplot.cornerplot(res, color='blue',
@@ -971,34 +1007,7 @@ fg, ax = dyplot.cornerplot(res, color='blue',
)
-# In[364]:
-
-
-afdist=samps_tr[-1][::downfactor]
-mdist=samps_tr[-2][::downfactor]
-
-w_dist=[None] * len(afdist)
-tau_dist=[None] * len(afdist)
-
-for i in range(len(afdist)):
- w_dist[i],tau_dist[i] = QNM_spectrum(mdist[i],afdist[i],2,2)
-
-for i in range(dim):
- w_lab[i] = wstr+str(i)+'$'
- tau_lab[i] = taustr+str(i)+'$'
-
-labels = np.concatenate((w_lab,tau_lab))
-truths = np.concatenate((w,tau))
-
-w_tau_dist=np.column_stack((w_dist,tau_dist))
-figure = corner.corner(w_tau_dist, labels=labels,
- quantiles=[0.05, 0.5, 0.95],
- show_titles=True, title_kwargs={"fontsize": 16},
- truths=truths,
- truth_color ='r')
-
-
-# In[365]:
+# In[52]:
if not eval(overwrite):
@@ -1007,7 +1016,7 @@ if not eval(overwrite):
figure.savefig(corner_plot_extra, format = 'png', bbox_inches = 'tight')
-# In[366]:
+# In[64]:
from dynesty import plotting as dyplot
@@ -1017,14 +1026,14 @@ fig, axes = dyplot.runplot(res, lnz_truth=lnz_truth)
fig.tight_layout()
-# In[367]:
+# In[54]:
if not eval(overwrite):
fig.savefig(diagnosis_plot, format = 'png', dpi = 384, bbox_inches = 'tight')
-# In[368]:
+# In[55]:
figband = plt.figure(figsize = (12, 9))
@@ -1042,14 +1051,14 @@ plt.ylabel('h')
plt.show()
-# In[371]:
+# In[56]:
if not eval(overwrite):
figband.savefig(fit_plot)
-# In[372]:
+# In[57]:
if not eval(overwrite):