Skip to content
Snippets Groups Projects
Commit be18f078 authored by frcojimenez's avatar frcojimenez
Browse files

updates on the m-spin sampling

parent cb4e509e
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
code_new/NR_dynesty_t0_loop.py
+ 262
− 253
View file @ be18f078
#!/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,16 +303,54 @@ 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
return (w_m_a, tau_m_a)
omegas_new = [qnm.modes_cache(s=-2,l=l,m=m,n=i)(a=af)[0] for i in range (0,dim)]
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]:
# In[297]:
np.sqrt(12/2*1/tauRD_to_t_NR(1.3*10**-47,70)*(5*10**(-21))**2)
# In[15]:
np.sqrt(0.004/2*1/(1.3*10**-47)*(5*10**(-21))**2)
# In[16]:
gw = {}
......@@ -264,12 +358,9 @@ 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:]
......@@ -277,8 +368,6 @@ elif nr_code=='Maya':
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:]
......@@ -293,7 +382,6 @@ elif nr_code=='Maya':
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'])
......@@ -318,7 +406,7 @@ elif nr_code=='LaZeV':
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]:
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)
# 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()
# In[325]:
print(time.process_time() - start)
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_'
# In[39]:
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))
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):
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment

AEI Privacy Policy Imprint