Skip to content
Snippets Groups Projects
Commit 68665660 authored by Francisco Jimenez Forteza's avatar Francisco Jimenez Forteza
Browse files

added the option of overwriting

parent 575cbeb1
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
code_new/source_code/RD_Fits.ipynb
+ 23
− 10
View file @ 68665660
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
"""Generate ringdown templates in the time and perform parameter estimation on them. """Generate ringdown templates in the time and perform parameter estimation on them.
""" """
``` ```
%% Output %% Output
'Generate ringdown templates in the time and perform parameter estimation on them.\n' 'Generate ringdown templates in the time and perform parameter estimation on them.\n'
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
#Import relevant modules, import data and all that #Import relevant modules, import data and all that
import time import time
import numpy as np import numpy as np
import corner import corner
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator from matplotlib.ticker import MaxNLocator
from matplotlib import rc from matplotlib import rc
from configparser import ConfigParser from configparser import ConfigParser
plt.rcParams.update({'font.size': 16.5}) plt.rcParams.update({'font.size': 16.5})
from multiprocessing import Pool from multiprocessing import Pool
import random import random
import dynesty import dynesty
from dynesty import plotting as dyplot from dynesty import plotting as dyplot
from dynesty.utils import resample_equal from dynesty.utils import resample_equal
from dynesty import utils as dyfunc from dynesty import utils as dyfunc
import os import os
import argparse import argparse
import scipy.optimize as optimization import scipy.optimize as optimization
from scipy.optimize import minimize from scipy.optimize import minimize
import rdown as rd import rdown as rd
import rdown_pe as rd_pe import rdown_pe as rd_pe
import rdown_utilities as rd_ut import rdown_utilities as rd_ut
import read_data as rdata import read_data as rdata
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
## Loading and running data tested with NR data ## Loading and running data tested with NR data
## Loading and running data tested with Mock data ## Loading and running data tested with Mock data
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Cell that calls the arguments from your 'config.ini' file. # Cell that calls the arguments from your 'config.ini' file.
try: try:
parser = argparse.ArgumentParser(description="Simple argument parser") parser = argparse.ArgumentParser(description="Simple argument parser")
parser.add_argument("-c", action="store", dest="config_file") parser.add_argument("-c", action="store", dest="config_file")
result = parser.parse_args() result = parser.parse_args()
config_file=result.config_file config_file=result.config_file
parser = ConfigParser() parser = ConfigParser()
parser.read(config_file) parser.read(config_file)
parser.sections() parser.sections()
except SystemExit: except SystemExit:
parser = ConfigParser() parser = ConfigParser()
parser.read('config_n0_to_1_mock.ini') parser.read('./run_0_mock/config_n0_to_1_mock.ini')
parser.sections() parser.sections()
pass pass
``` ```
%% Output %% Output
usage: ipykernel_launcher.py [-h] [-c CONFIG_FILE] usage: ipykernel_launcher.py [-h] [-c CONFIG_FILE]
ipykernel_launcher.py: error: unrecognized arguments: -f /work/francisco.jimenez/.local/share/jupyter/runtime/kernel-eb21fd9f-ea57-481c-b51d-230c5ef5a4a2.json ipykernel_launcher.py: error: unrecognized arguments: -f /work/francisco.jimenez/.local/share/jupyter/runtime/kernel-eb21fd9f-ea57-481c-b51d-230c5ef5a4a2.json
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Load variables from config file # Load variables from config file
(simulation_path_1,simulation_path_2, metadata_file , simulation_number, output_folder, (simulation_path_1,simulation_path_2, metadata_file , simulation_number, output_folder,
export, overwrite, sampler,nr_code, nbcores,tshift,tend,t_align, export, overwrite, sampler,nr_code, nbcores,tshift,tend,t_align,
nmax , npoints, model, error_str, fitnoise, l_int, index_mass,index_spin, nmax , npoints, model, error_str, fitnoise, l_int, index_mass,index_spin,
error_type, error_val, af, mf,tau_var_str,nm_mock)=rdata.read_config_file(parser) error_type, error_val, af, mf,tau_var_str,nm_mock)=rdata.read_config_file(parser)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Show configuration options # Show configuration options
dim = nmax+1 dim = nmax+1
ndim = 4*dim ndim = 4*dim
numbins = 32 #corner plot parameter - how many bins you want numbins = 32 #corner plot parameter - how many bins you want
print('model:',model) print('model:',model)
print('nmax:',nmax) print('nmax:',nmax)
print('nm_mock:',nm_mock) print('nm_mock:',nm_mock)
print('tshift:',tshift) print('tshift:',tshift)
print('error:', error_str) print('error:', error_str)
print('error value:',error_val) print('error value:',error_val)
print('export:',export) print('export:',export)
print('nr code:',nr_code) print('nr code:',nr_code)
print('fit noise:',fitnoise) print('fit noise:',fitnoise)
``` ```
%% Output %% Output
model: w-tau model: w-tau
nmax: 0 nmax: 0
nm_mock: 1 nm_mock: 1
tshift: 0.0 tshift: 0.2
error: True error: True
error value: 0.002 error value: 0.002
export: True export: True
nr code: Mock-data nr code: Mock-data
fit noise: False fit noise: False
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Create output directories # Create output directories
if not os.path.exists(output_folder): if not os.path.exists(output_folder):
os.mkdir(output_folder) os.mkdir(output_folder)
print("Directory " , output_folder , " Created ") print("Directory " , output_folder , " Created ")
if nr_code == 'Mock-data': if nr_code == 'Mock-data':
nm_mock_str = 'rec_with'+parser.get('rd-mock-parameters','nm_mock')+'_' nm_mock_str = 'rec_with'+parser.get('rd-mock-parameters','nm_mock')+'_'
else: else:
nm_mock_str='' nm_mock_str=''
if error_str: if error_str:
output_folder_1=(output_folder+'/'+model+'-nmax'+str(nmax)+'_'+nm_mock_str+str(error_str)+'_'+str(error_type)+'_fitnoise_'+str(fitnoise)) output_folder_1=(output_folder+'/'+model+'-nmax'+str(nmax)+'_'+nm_mock_str+str(error_str)+'_'+str(error_type)+'_fitnoise_'+str(fitnoise))
else: else:
output_folder_1=output_folder+'/'+model+'-nmax'+str(nmax)+'_'+nm_mock_str+str(error_str)+'_'+'fitnoise_'+str(fitnoise) output_folder_1=output_folder+'/'+model+'-nmax'+str(nmax)+'_'+nm_mock_str+str(error_str)+'_'+'fitnoise_'+str(fitnoise)
if not os.path.exists(output_folder_1): if not os.path.exists(output_folder_1):
os.mkdir(output_folder_1) os.mkdir(output_folder_1)
print("Directory " , output_folder_1 , " Created ") print("Directory " , output_folder_1 , " Created ")
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Define output files # Define output files
pars = [simulation_number,model,nmax,tshift,npoints] pars = [simulation_number,model,nmax,tshift,npoints]
corner_plot = rdata.create_output_files(output_folder_1,pars,'corner_plot') corner_plot = rdata.create_output_files(output_folder_1,pars,'corner_plot')
corner_plot_extra = rdata.create_output_files(output_folder_1,pars,'corner_plot_extra') corner_plot_extra = rdata.create_output_files(output_folder_1,pars,'corner_plot_extra')
diagnosis_plot = rdata.create_output_files(output_folder_1,pars,'diagnosis') diagnosis_plot = rdata.create_output_files(output_folder_1,pars,'diagnosis')
fit_plot = rdata.create_output_files(output_folder_1,pars,'fit') fit_plot = rdata.create_output_files(output_folder_1,pars,'fit')
samples_file = rdata.create_output_files(output_folder_1,pars,'post_samples') samples_file = rdata.create_output_files(output_folder_1,pars,'post_samples')
results_file = rdata.create_output_files(output_folder_1,pars,'sampler_results') results_file = rdata.create_output_files(output_folder_1,pars,'sampler_results')
sumary_data = rdata.create_output_files(output_folder_1,pars,'log_z') sumary_data = rdata.create_output_files(output_folder_1,pars,'log_z')
best_data = rdata.create_output_files(output_folder_1,pars,'best_vals') best_data = rdata.create_output_files(output_folder_1,pars,'best_vals')
files = [corner_plot,corner_plot_extra,diagnosis_plot,fit_plot,samples_file,results_file,sumary_data,best_data]
```
%% Cell type:code id: tags:
``` python
# Remove old files if overwrite = True
if overwrite:
rd_ut.rm_files(files)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
#Load NR data, align in time and resize. Plot real part and amplitude. Finally compute the mismatch and the snr estimate #Load NR data, align in time and resize. Plot real part and amplitude. Finally compute the mismatch and the snr estimate
data = rdata.read_data(nr_code,simulation_path_1,RD=True,tshift=tshift,tend = tend,metadata_file=metadata_file,parser=parser) data = rdata.read_data(nr_code,simulation_path_1,RD=True,tshift=tshift,tend = tend,metadata_file=metadata_file,parser=parser)
data_l = rdata.read_data(nr_code,simulation_path_2,RD=True,tshift=tshift,tend = tend,metadata_file=metadata_file,parser=parser) data_l = rdata.read_data(nr_code,simulation_path_2,RD=True,tshift=tshift,tend = tend,metadata_file=metadata_file,parser=parser)
data_r, data_lr = rdata.nr_resize(data,data_l,tshift=tshift,tend=tend) data_r, data_lr = rdata.nr_resize(data,data_l,tshift=tshift,tend=tend)
times_rd = data_r[:,0] times_rd = data_r[:,0]
plt.figure(figsize = (12, 8)) plt.figure(figsize = (12, 8))
plt.plot(times_rd, data_r[:,1].real, "r", alpha=0.3, lw=3, label=r'$Lev6$: real') plt.plot(times_rd, data_r[:,1].real, "r", alpha=0.3, lw=3, label=r'$Lev6$: real')
plt.plot(times_rd, np.sqrt((data_r[:,1].real)**2+(data_r[:,1].imag)**2), "r", alpha=0.3, lw=3, label=r'$Lev5\,amp$') plt.plot(times_rd, np.sqrt((data_r[:,1].real)**2+(data_r[:,1].imag)**2), "r", alpha=0.3, lw=3, label=r'$Lev5\,amp$')
plt.plot(times_rd, data_lr[:,1].real, "b", alpha=0.3, lw=3, label=r'$Lev5: real$') plt.plot(times_rd, data_lr[:,1].real, "b", alpha=0.3, lw=3, label=r'$Lev5: real$')
plt.plot(times_rd, np.sqrt((data_lr[:,1].real)**2+(data_lr[:,1].imag)**2), "b", alpha=0.3, lw=3, label=r'$Lev5\,amp$') plt.plot(times_rd, np.sqrt((data_lr[:,1].real)**2+(data_lr[:,1].imag)**2), "b", alpha=0.3, lw=3, label=r'$Lev5\,amp$')
if error_str and error_val==0: if error_str and error_val==0:
error = np.sqrt(data_r[:,1]*data_r[:,1]-2*data_r[:,1]*data_lr[:,1]+data_lr[:,1]*data_lr[:,1]) error = np.sqrt(data_r[:,1]*data_r[:,1]-2*data_r[:,1]*data_lr[:,1]+data_lr[:,1]*data_lr[:,1])
error_est=np.sqrt(error.imag**2+error.real**2) error_est=np.sqrt(error.imag**2+error.real**2)
plt.plot(times_rd, error_est, "g", alpha=0.3, lw=2, label='error') plt.plot(times_rd, error_est, "g", alpha=0.3, lw=2, label='error')
plt.legend() plt.legend()
mismatch=1-rd_ut.EasyMatchT(times_rd,data_r[:,1],data_lr[:,1],tshift,tend) mismatch=1-rd_ut.EasyMatchT(times_rd,data_r[:,1],data_lr[:,1],tshift,tend)
error=np.sqrt(2*mismatch) error=np.sqrt(2*mismatch)
print('error estimate:',error) print('error estimate:',error)
print('mismatch:', mismatch) print('mismatch:', mismatch)
print('snr:', rd_ut.EasySNRT(times_rd,data_r[:,1],data_lr[:,1],tshift,tend)/error**2) print('snr:', rd_ut.EasySNRT(times_rd,data_r[:,1],data_lr[:,1],tshift,tend)/error**2)
``` ```
%% Output %% Output
/work/francisco.jimenez/venv/lib/python3.7/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part /work/francisco.jimenez/venv/lib/python3.7/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
return array(a, dtype, copy=False, order=order) return array(a, dtype, copy=False, order=order)
/work/francisco.jimenez/venv/lib/python3.7/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part /work/francisco.jimenez/venv/lib/python3.7/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
return array(a, dtype, copy=False, order=order) return array(a, dtype, copy=False, order=order)
/work/francisco.jimenez/venv/lib/python3.7/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part /work/francisco.jimenez/venv/lib/python3.7/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
return array(a, dtype, copy=False, order=order) return array(a, dtype, copy=False, order=order)
/work/francisco.jimenez/venv/lib/python3.7/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part /work/francisco.jimenez/venv/lib/python3.7/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
return array(a, dtype, copy=False, order=order) return array(a, dtype, copy=False, order=order)
/work/francisco.jimenez/venv/lib/python3.7/site-packages/ipykernel_launcher.py:19: RuntimeWarning: invalid value encountered in sqrt /work/francisco.jimenez/venv/lib/python3.7/site-packages/ipykernel_launcher.py:19: RuntimeWarning: invalid value encountered in sqrt
error estimate: nan error estimate: nan
mismatch: -2.220446049250313e-16 mismatch: -2.220446049250313e-16
snr: nan snr: nan
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Phase alignement # Phase alignement
if parser.has_option('rd-model','phase_alignment'): if parser.has_option('rd-model','phase_alignment'):
phase_alignment=eval(parser.get('rd-model','phase_alignment')) phase_alignment=eval(parser.get('rd-model','phase_alignment'))
else: else:
phase_alignment=False phase_alignment=False
if phase_alignment: if phase_alignment:
datar_al = rdata.phase_align(data_r,data_lr) datar_al = rdata.phase_align(data_r,data_lr)
gwdatanew5 = data_lr[:,1] gwdatanew5 = data_lr[:,1]
gwdatanew = datar_al[:,1] gwdatanew = datar_al[:,1]
timesrd_final = datar_al[:,0] timesrd_final = datar_al[:,0]
mismatch=1-rd_ut.EasyMatchT(timesrd_final,gwdatanew,gwdatanew5,tshift,tend) mismatch=1-rd_ut.EasyMatchT(timesrd_final,gwdatanew,gwdatanew5,tshift,tend)
error=np.sqrt(2*mismatch) error=np.sqrt(2*mismatch)
print('error estimate:',error) print('error estimate:',error)
print('mismatch:', mismatch) print('mismatch:', mismatch)
print('snr:', rd_ut.EasySNRT(timesrd_final,gwdatanew,gwdatanew5,tshift,tend)/error) print('snr:', rd_ut.EasySNRT(timesrd_final,gwdatanew,gwdatanew5,tshift,tend)/error)
if error_str: if error_str:
error = np.sqrt(gwdatanew*gwdatanew-2*gwdatanew*gwdatanew5+gwdatanew5*gwdatanew5) error = np.sqrt(gwdatanew*gwdatanew-2*gwdatanew*gwdatanew5+gwdatanew5*gwdatanew5)
error_est=np.sqrt(error.imag**2+error.real**2) error_est=np.sqrt(error.imag**2+error.real**2)
else : else :
error = 1 error = 1
else: else:
datar_al = data_r datar_al = data_r
timesrd_final = datar_al[:,0] timesrd_final = datar_al[:,0]
#Test the new interpolated data #Test the new interpolated data
if error_str and error_val==0: if error_str and error_val==0:
plt.figure(figsize = (12, 8)) plt.figure(figsize = (12, 8))
plt.plot(timesrd_final, datar_al[:,1].real, "r", alpha=0.3, lw=2, label='Original') plt.plot(timesrd_final, datar_al[:,1].real, "r", alpha=0.3, lw=2, label='Original')
plt.plot(timesrd_final, data_lr[:,1].real, "b", alpha=0.3, lw=2, label='Aligned') plt.plot(timesrd_final, data_lr[:,1].real, "b", alpha=0.3, lw=2, label='Aligned')
plt.plot(timesrd_final, error_est, "b", alpha=0.3, lw=2, label='error') plt.plot(timesrd_final, error_est, "b", alpha=0.3, lw=2, label='error')
plt.legend() plt.legend()
``` ```
%% Output %% Output
error estimate: nan error estimate: nan
mismatch: -2.220446049250313e-16 mismatch: -2.220446049250313e-16
snr: nan snr: nan
/work/francisco.jimenez/venv/lib/python3.7/site-packages/ipykernel_launcher.py:13: RuntimeWarning: invalid value encountered in sqrt /work/francisco.jimenez/venv/lib/python3.7/site-packages/ipykernel_launcher.py:13: RuntimeWarning: invalid value encountered in sqrt
del sys.path[0] del sys.path[0]
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Define your noise depending on the noise configuration. Load priors and setup the likelihood with rd_pe.Ringdown_PE. # Define your noise depending on the noise configuration. Load priors and setup the likelihood with rd_pe.Ringdown_PE.
if error_str and error_val==0: if error_str and error_val==0:
error_final = error_est error_final = error_est
norm_factor = 100*len(error_final)/2*np.log(2*np.pi) norm_factor = 100*len(error_final)/2*np.log(2*np.pi)
elif error_str and error_val!=0: elif error_str and error_val!=0:
datar_al[:,1]+=random.uniform(0, error_val) datar_al[:,1]+=random.uniform(0, error_val)
datar_al[:,1]+=1j*random.uniform(0, error_val) datar_al[:,1]+=1j*random.uniform(0, error_val)
error_tsh = error_val error_tsh = error_val
error_final=(error_tsh.real**2+error_tsh.imag**2) error_final=(error_tsh.real**2+error_tsh.imag**2)
norm_factor = 0 norm_factor = 0
else: else:
error_tsh=1 error_tsh=1
error_final=(error_tsh.real**2+error_tsh.imag**2) error_final=(error_tsh.real**2+error_tsh.imag**2)
norm_factor = 0 norm_factor = 0
priors = rd_pe.load_priors(model,parser,nmax,fitnoise=fitnoise) priors = rd_pe.load_priors(model,parser,nmax,fitnoise=fitnoise)
rdown=rd.Ringdown_Spectrum(mf,af,2,2,n=nmax,s=-2,time=timesrd_final) rdown=rd.Ringdown_Spectrum(mf,af,2,2,n=nmax,s=-2,time=timesrd_final)
rdown_pe = rd_pe.Ringdown_PE(rdown,datar_al,dim,priors,errors2=error_final,norm_factor=norm_factor,model=model,l_int=l_int) rdown_pe = rd_pe.Ringdown_PE(rdown,datar_al,dim,priors,errors2=error_final,norm_factor=norm_factor,model=model,l_int=l_int)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Get a first estimate by trying to fit the data. # Get a first estimate by trying to fit the data.
nll = lambda *args: -rdown_pe.log_likelihood(*args) nll = lambda *args: -rdown_pe.log_likelihood(*args)
if model == 'w-tau-fixed-m-af': if model == 'w-tau-fixed-m-af':
if fitnoise: if fitnoise:
initial = np.concatenate((np.ones(2*dim),[0.8,0.9,1])) initial = np.concatenate((np.ones(2*dim),[0.8,0.9,1]))
soln = minimize(nll, initial,bounds=priors) soln = minimize(nll, initial,bounds=priors)
vars_ml=soln.x vars_ml=soln.x
else: else:
initial = np.concatenate((np.ones(2*dim),[0.8,0.9])) initial = np.concatenate((np.ones(2*dim),[0.8,0.9]))
soln = minimize(nll, initial,bounds=priors) soln = minimize(nll, initial,bounds=priors)
vars_ml=soln.x vars_ml=soln.x
elif model == 'w-tau-fixed': elif model == 'w-tau-fixed':
if fitnoise: if fitnoise:
initial = np.concatenate((np.ones(2*dim),[0.2])) initial = np.concatenate((np.ones(2*dim),[0.2]))
soln = minimize(nll, initial,bounds=priors) soln = minimize(nll, initial,bounds=priors)
vars_ml=soln.x vars_ml=soln.x
else: else:
initial = np.ones(2*dim) initial = np.ones(2*dim)
soln = minimize(nll, initial,bounds=priors) soln = minimize(nll, initial,bounds=priors)
vars_ml=soln.x vars_ml=soln.x
else: else:
if fitnoise: if fitnoise:
initial = np.concatenate((np.ones(ndim),[1])) initial = np.concatenate((np.ones(ndim),[1]))
soln = minimize(nll, initial,bounds=priors) soln = minimize(nll, initial,bounds=priors)
vars_ml=soln.x vars_ml=soln.x
else: else:
initial = np.ones(ndim) initial = np.ones(ndim)
soln = minimize(nll, initial,bounds=priors) soln = minimize(nll, initial,bounds=priors)
vars_ml=soln.x vars_ml=soln.x
print("best fit pars from fit: ",vars_ml) print("best fit pars from fit: ",vars_ml)
``` ```
%% Output %% Output
best fit pars from fit: [ 0.2 10.52057151 0.08836181 6.28318531] best fit pars from fit: [ 0.2 10.52057151 0.08836181 6.28318531]
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
mypool = Pool(nbcores) mypool = Pool(nbcores)
mypool.size = nbcores mypool.size = nbcores
start = time.process_time() start = time.process_time()
f2=dynesty.NestedSampler(rdown_pe.log_likelihood,rdown_pe.prior_transform, len(priors), nlive=npoints,sample=sampler,pool=mypool) f2=dynesty.NestedSampler(rdown_pe.log_likelihood,rdown_pe.prior_transform, len(priors), nlive=npoints,sample=sampler,pool=mypool)
if parser.has_option('setup','dlogz'): if parser.has_option('setup','dlogz'):
dlogz=np.float(parser.get('setup','dlogz')) dlogz=np.float(parser.get('setup','dlogz'))
f2.run_nested(dlogz=dlogz,print_progress=False) f2.run_nested(dlogz=dlogz,print_progress=False)
else: else:
f2.run_nested(print_progress=False) f2.run_nested(print_progress=False)
print(time.process_time() - start) print(time.process_time() - start)
``` ```
%% Output %% Output
46474it [21:02, 36.81it/s, +1000 | bound: 286 | nc: 1 | ncall: 1106811 | eff(%): 4.289 | loglstar: -inf < -959286.172 < inf | logz: -959328.014 +/- 0.280 | dlogz: 0.000 > 0.010] 46474it [21:02, 36.81it/s, +1000 | bound: 286 | nc: 1 | ncall: 1106811 | eff(%): 4.289 | loglstar: -inf < -959286.172 < inf | logz: -959328.014 +/- 0.280 | dlogz: 0.000 > 0.010]
645.863601611 645.863601611
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
res = f2.results res = f2.results
res.samples_u.shape res.samples_u.shape
res.summary() res.summary()
samps=f2.results.samples samps=f2.results.samples
postsamps = rd_ut.posterior_samples(f2) postsamps = rd_ut.posterior_samples(f2)
samps_tr=np.transpose(samps) samps_tr=np.transpose(samps)
half_points=int(round((len(samps_tr[0])/1.25))) half_points=int(round((len(samps_tr[0])/1.25)))
evidence = res.logz[-1] evidence = res.logz[-1]
evidence_error = res.logzerr[-1] evidence_error = res.logzerr[-1]
if export: if export:
rd_ut.save_object(res, results_file) rd_ut.save_object(res, results_file)
``` ```
%% Output %% Output
Summary Summary
======= =======
nlive: 1000 nlive: 1000
niter: 46474 niter: 46474
ncall: 1106811 ncall: 1106811
eff(%): 4.289 eff(%): 4.289
logz: -959328.014 +/- 0.280 logz: -959328.014 +/- 0.280
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
pars = nmax,model,samps_tr, half_points pars = nmax,model,samps_tr, half_points
npamps = rd_ut.get_best_amps(pars,parser=parser,nr_code=nr_code) npamps = rd_ut.get_best_amps(pars,parser=parser,nr_code=nr_code)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
if export: if export:
pars = simulation_number, nmax, tshift, evidence, evidence_error pars = simulation_number, nmax, tshift, evidence, evidence_error
rd_ut.export_logz_files(sumary_data,pars) rd_ut.export_logz_files(sumary_data,pars)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
labels = rd_ut.define_labels(dim,model,fitnoise) labels = rd_ut.define_labels(dim,model,fitnoise)
if export: if export:
pars = tshift, len(priors), labels pars = tshift, len(priors), labels
rd_ut.export_bestvals_files(best_data,postsamps,pars) rd_ut.export_bestvals_files(best_data,postsamps,pars)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
w, tau = rdown.QNM_spectrum() w, tau = rdown.QNM_spectrum()
pars = w, tau, mf, af, npamps pars = w, tau, mf, af, npamps
truths = rd_ut.get_truths(model,pars,fitnoise) truths = rd_ut.get_truths(model,pars,fitnoise)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fg=corner.corner(postsamps,quantiles=[0.05,0.5,0.95],show_titles=True,max_n_ticks = 4,bins=50,truths=truths,labels=labels,truth_color='red') fg=corner.corner(postsamps,quantiles=[0.05,0.5,0.95],show_titles=True,max_n_ticks = 4,bins=50,truths=truths,labels=labels,truth_color='red')
plt.show() plt.show()
if export: if export:
fg.savefig(corner_plot, format = 'png', bbox_inches = 'tight') fg.savefig(corner_plot, format = 'png', bbox_inches = 'tight')
``` ```
%% Output %% Output
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
from importlib import reload from importlib import reload
reload(rd_ut) reload(rd_ut)
if model == 'w-tau-fixed-m-af' and export == True: if model == 'w-tau-fixed-m-af' and export == True:
truths=np.concatenate((w,tau)) truths=np.concatenate((w,tau))
labels_mf = np.concatenate((w_lab,tau_lab)) labels_mf = np.concatenate((w_lab,tau_lab))
new_samples = rd_ut.convert_m_af_2_w_tau_post(res,fitnoise=False) new_samples = rd_ut.convert_m_af_2_w_tau_post(res,fitnoise=False)
figure = corner.corner(new_samples,truths=truths,quantiles=[0.05,0.95],labels=labels_mf,smooth=True,color='b',truth_color='r',show_titles=True) figure = corner.corner(new_samples,truths=truths,quantiles=[0.05,0.95],labels=labels_mf,smooth=True,color='b',truth_color='r',show_titles=True)
figure.savefig(corner_plot_extra, format = 'png', bbox_inches = 'tight') figure.savefig(corner_plot_extra, format = 'png', bbox_inches = 'tight')
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
#lnz_truth = ndim * -np.log(2 * 10.) # analytic evidence solution #lnz_truth = ndim * -np.log(2 * 10.) # analytic evidence solution
fig, axes = dyplot.runplot(res) fig, axes = dyplot.runplot(res)
fig.tight_layout() fig.tight_layout()
if export: if export:
fig.savefig(diagnosis_plot, format = 'png', dpi = 384, bbox_inches = 'tight') fig.savefig(diagnosis_plot, format = 'png', dpi = 384, bbox_inches = 'tight')
``` ```
%% Output %% Output
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
if export: if export:
dict = {'w-tau':rdown.rd_model_wtau , 'w-q': rdown.rd_model_wq, 'w-tau-fixed':rdown.rd_model_wtau_fixed,'w-tau-fixed-m-af': rdown.rd_model_wtau_m_af} dict = {'w-tau':rdown.rd_model_wtau , 'w-q': rdown.rd_model_wq, 'w-tau-fixed':rdown.rd_model_wtau_fixed,'w-tau-fixed-m-af': rdown.rd_model_wtau_m_af}
figband = plt.figure(figsize = (12, 9)) figband = plt.figure(figsize = (12, 9))
plt.plot(datar_al[:,0].real,datar_al[:,1].real, "green", alpha=0.9, lw=3, label=r'$res_{240}$') plt.plot(datar_al[:,0].real,datar_al[:,1].real, "green", alpha=0.9, lw=3, label=r'$res_{240}$')
onesig_bounds = np.array([np.percentile(postsamps[:, i], [5, 95]) for i in range(len(postsamps[0]))]).T onesig_bounds = np.array([np.percentile(postsamps[:, i], [5, 95]) for i in range(len(postsamps[0]))]).T
samples_1sigma = filter(lambda sample: np.all(onesig_bounds[0] <= sample) and np.all(sample <= onesig_bounds[1]), postsamps) samples_1sigma = filter(lambda sample: np.all(onesig_bounds[0] <= sample) and np.all(sample <= onesig_bounds[1]), postsamps)
samples_1sigma_down = list(samples_1sigma)[::downfactor] samples_1sigma_down = list(samples_1sigma)[::downfactor]
for sample in samples_1sigma_down: for sample in samples_1sigma_down:
plt.plot(datar_al[:,0].real, dict[model](sample).real, "r-", alpha=0.1, lw=1) plt.plot(datar_al[:,0].real, 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.title(r'Comparison of the MC fit data and the $1-\sigma$ error band')
plt.legend() plt.legend()
plt.xlabel('t') plt.xlabel('t')
plt.ylabel('h') plt.ylabel('h')
plt.show() plt.show()
figband.savefig(fit_plot) figband.savefig(fit_plot)
``` ```
%% Output %% Output
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
if export: if export:
with open(samples_file,'w') as file: with open(samples_file,'w') as file:
writer = csv.writer(file) writer = csv.writer(file)
writer.writerow(labels) writer.writerow(labels)
writer.writerows(samps[::downfactor]) writer.writerows(samps[::downfactor])
``` ```
......
code_new/source_code/RD_Fits.py
+ 18
− 8
View file @ 68665660
...@@ -8,7 +8,7 @@ ...@@ -8,7 +8,7 @@
""" """
# In[39]: # In[1]:
#Import relevant modules, import data and all that #Import relevant modules, import data and all that
...@@ -37,14 +37,14 @@ import rdown_utilities as rd_ut ...@@ -37,14 +37,14 @@ import rdown_utilities as rd_ut
import read_data as rdata import read_data as rdata
# In[40]: # In[2]:
## Loading and running data tested with NR data ## Loading and running data tested with NR data
## Loading and running data tested with Mock data ## Loading and running data tested with Mock data
# In[41]: # In[3]:
# Cell that calls the arguments from your 'config.ini' file. # Cell that calls the arguments from your 'config.ini' file.
...@@ -58,12 +58,12 @@ try: ...@@ -58,12 +58,12 @@ try:
parser.sections() parser.sections()
except SystemExit: except SystemExit:
parser = ConfigParser() parser = ConfigParser()
parser.read('config_n0_to_1_mock.ini') parser.read('./run_0_mock/config_n0_to_1_mock.ini')
parser.sections() parser.sections()
pass pass
# In[42]: # In[4]:
# Load variables from config file # Load variables from config file
...@@ -73,7 +73,7 @@ except SystemExit: ...@@ -73,7 +73,7 @@ except SystemExit:
error_type, error_val, af, mf,tau_var_str,nm_mock)=rdata.read_config_file(parser) error_type, error_val, af, mf,tau_var_str,nm_mock)=rdata.read_config_file(parser)
# In[43]: # In[5]:
# Show configuration options # Show configuration options
...@@ -92,7 +92,7 @@ print('nr code:',nr_code) ...@@ -92,7 +92,7 @@ print('nr code:',nr_code)
print('fit noise:',fitnoise) print('fit noise:',fitnoise)
# In[44]: # In[6]:
# Create output directories # Create output directories
...@@ -115,7 +115,7 @@ if not os.path.exists(output_folder_1): ...@@ -115,7 +115,7 @@ if not os.path.exists(output_folder_1):
print("Directory " , output_folder_1 , " Created ") print("Directory " , output_folder_1 , " Created ")
# In[45]: # In[7]:
# Define output files # Define output files
...@@ -129,6 +129,16 @@ results_file = rdata.create_output_files(output_folder_1,pars,'sampler_results') ...@@ -129,6 +129,16 @@ results_file = rdata.create_output_files(output_folder_1,pars,'sampler_results')
sumary_data = rdata.create_output_files(output_folder_1,pars,'log_z') sumary_data = rdata.create_output_files(output_folder_1,pars,'log_z')
best_data = rdata.create_output_files(output_folder_1,pars,'best_vals') best_data = rdata.create_output_files(output_folder_1,pars,'best_vals')
files = [corner_plot,corner_plot_extra,diagnosis_plot,fit_plot,samples_file,results_file,sumary_data,best_data]
# In[8]:
# Remove old files if overwrite = True
if overwrite:
rd_ut.rm_files(files)
# In[46]: # In[46]:
......
code_new/source_code/Run_Setup.ipynb 0 → 100644
+ 167
− 0
View file @ 68665660
%% Cell type:code id:square-toddler tags:
``` python
"""
Created by Sumit Kumar on 2020-03-08 && modified by Xisco on 2021-04
Last modified:
"""
import os, sys, numpy, glob, argparse
from datetime import date
import subprocess
from subprocess import call
import re
from configparser import ConfigParser
today = date.today()
date = today.strftime("%Y%m%d")
runname='run_0_mock'
nmax = 0
config_file ='config_n0_to_1_mock.ini'
overwrite = True
######
times = '(0 0.2 0.4 0.8 1.2 2 2.5 5 7.5 10 12 15 18 20)'
accounting_group = 'cbc.test.pe_ringdown'
cpus=8
nlive_points = 2000
req_memory="16GB"
not_user='frjifo@aei.mpg.de'
pythonfile='/work/francisco.jimenez/sio/git/rdstackingproject/code_new/RD_Fits.py'
pythonscript='/work/francisco.jimenez/venv/bin/python'
#######################################################
pwd = os.getcwd()
run_dir = '%s/%s'%(pwd,runname)
logs_dir = '%s/logs'%run_dir
os.system('mkdir -p %s'%logs_dir)
os.system('cp %s %s/'%(config_file,run_dir))
###########################################################################
# Creating Condor submit file
###########################################################################
filename1 = '%s/%s'%(run_dir,'condor_submit')
text_file1 = open(filename1 + ".sub", "w")
text_file1.write("universe = vanilla\n")
text_file1.write("getenv = true\n")
text_file1.write("# run script -- make sure that condor has execute permission for this file (chmod a+x script.py)\n")
text_file1.write("executable = "'%s/%s'%(run_dir,runname+'.sh \n'))
text_file1.write("# file to dump stdout (this directory should exist)\n")
text_file1.write("output = %s/%s-$(Process).out\n"%(logs_dir,runname))
text_file1.write("# file to dump stderr\n")
text_file1.write("error = %s/%s-$(Process).err\n"%(logs_dir,runname))
text_file1.write("# condor logs\n")
text_file1.write("log = %s/%s-$(Process).log\n"%(logs_dir,runname))
text_file1.write("initialdir = %s \n"%run_dir)
text_file1.write("notify_user ="+not_user+' \n')
text_file1.write("notification = Complete\n")
text_file1.write('''arguments = "-processid $(Process)" \n''')
text_file1.write("request_memory = "+str(req_memory)+"\n")
text_file1.write("request_cpus = "+str(cpus)+"\n")
text_file1.write("on_exit_remove = (ExitBySignal == False) || ((ExitBySignal == True) && (ExitSignal != 11))\n")
text_file1.write("accounting_group = %s\n"%accounting_group)
text_file1.write("queue 1\n")
text_file1.write("\n")
text_file1.close()
###########################################################
# Creating python executable file
############################################################
filename2 = run_dir+'/'+runname+'.sh'
text_file2 = open(filename2, "w")
text_file2.write("#! /bin/bash \n")
text_file2.write("\n")
text_file2.write("times="+times+"\n")
text_file2.write("config_file="+config_file+"\n")
text_file2.write("pythonfile="+pythonfile+"\n")
text_file2.write("pythonscript="+pythonscript+"\n")
text_file2.write("\n")
text_file2.write("for i in ${times[@]}; do\n")
text_file2.write(" awk -v a=\"$i\" '/^tshift/ && $3 != \"supplied\" { $3=a } { print }' $config_file > tmp && mv tmp $config_file\n")
text_file2.write(" $pythonscript $pythonfile -c $config_file \n")
text_file2.write("done\n")
text_file2.write("awk -v a=\"0\" '/^tshift/ && $3 != \"supplied\" { $3=a } { print }' $config_file > tmp && mv tmp $config_file \n")
text_file2.write("\n")
text_file2.close()
os.system('chmod u+x %s'%filename2)
os.system('cp '+runname+'.sh %s/'%run_dir)
os.system('chmod u+x ./'+runname+'.sh')
os.system('cd '+run_dir)
###########################################################
# Checking the configuration file and adding some important replacements
############################################################
filename3 = '%s/%s'%(run_dir,config_file)
# change the number of cores
bashCommand = "awk -v a="+str(cpus)+" '/^nb_cores/ && $3 != \"supplied\" { $3=a } { print }' "+str(config_file)+" > tmp && mv tmp "+str(config_file);
subprocess.call(bashCommand,shell=True)
filename3 = '%s/%s'%(run_dir,config_file)
# change the number of nmax
bashCommand = "awk -v a="+str(nmax)+" '/^nmax/ && $3 != \"supplied\" { $3=a } { print }' "+str(config_file)+" > tmp && mv tmp "+str(config_file);
subprocess.call(bashCommand,shell=True)
filename3 = '%s/%s'%(run_dir,config_file)
# change the number of nmax
bashCommand = "awk -v a="+str(nlive_points)+" '/^npoints/ && $3 != \"supplied\" { $3=a } { print }' "+str(config_file)+" > tmp && mv tmp "+str(config_file);
subprocess.call(bashCommand,shell=True)
# change the overwrite parameter
bashCommand = "awk -v a="+str(overwrite)+" '/^overwrite/ && $3 != \"supplied\" { $3=a } { print }' "+str(config_file)+" > tmp && mv tmp "+str(config_file);
subprocess.call(bashCommand,shell=True)
###########################################################
# Submit the job
############################################################
filename4 = '%s/%s'%(run_dir,'condor_submit.sub')
bashCommand = ['condor_submit', str(filename4)]
process = subprocess.Popen(bashCommand, stdout=subprocess.PIPE,stderr=subprocess.PIPE);
output,error = process.communicate()
error
```
%% Output
b''
code_new/source_code/rdown.py 0 → 100644
+ 208
− 0
View file @ 68665660
# Copyright (C) 2021 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
#
# =============================================================================
#
# Module to generate RD waveforms.
import numpy as np
import qnm
import os
f_fpars= [[2.95845, -2.58697, 0.0533469], [2.12539, -1.78054, 0.0865503], [1.74755, -1.44776, 0.123666], [1.78287, -1.53203, 0.129475], [2.04028, -1.83224, 0.112497]]
q_fpars=[[0.584077, 1.52053, -0.480658], [0.00561441, 0.630715, -0.432664], [-0.197965, 0.515956, -0.369706], [-0.275097, 0.455691, -0.331543], [-0.287596, 0.398514, -0.309799]]
c=2.99792458*10**8;G=6.67259*10**(-11);MS=1.9885*10**30;
class Ringdown_Spectrum:
"""RDown model generator"""
def __init__(self,mf,af,l,m,n=4,s=-2,time=[],fixed=False,qnm_model='berti'):
self.mf = mf
self.af = af
self.l = l
self.m = m
self.n = n
self.time = time
self.grav_220 = [qnm.modes_cache(s=s,l=self.l,m=self.m,n=i) for i in range (0,self.n+1)]
self.dim = self.n+1
self.fixed = fixed
self.qnm_model = qnm_model
dict_omega = {'berti': self.QNM_Berti , 'qnm': self.QNM_spectrum}
dic = {'w-tau':self.rd_model_wtau , 'w-q': self.rd_model_wq, 'w-tau-fixed':self.rd_model_wtau_fixed,'w-tau-fixed-m-af': self.rd_model_wtau_m_af}
if len(self.time)==0:
self.time = np.arange(0,100,0.1)
if self.fixed:
omegas_new=np.asarray([self.grav_220[i](a=self.af)[0] for i in range (0,self.dim)])
self.w = (np.real(omegas_new))/self.mf
self.tau=-1/(np.imag(omegas_new))*self.mf
def QNM_spectrum(self):
""" It computes the RD frequencies and damping times in NR units.
"""
omegas_new=np.asarray([self.grav_220[i](a=self.af)[0] for i in range (0,self.n+1)])
w_m_a = (np.real(omegas_new))/self.mf
tau_m_a=-1/(np.imag(omegas_new))*self.mf
return (w_m_a, tau_m_a)
def QNM_Berti(self,rdowndata):
""" It computes the RD frequencies and damping times in NR units.
"""
position=np.argmax(rdowndata[0,0] >= (self.af))
#w_m_a=f1+f2*(1-af)**f3
w_m_a=[None]*(self.n+1)
tau_ma_a=[None]*(self.n+1)
for i in range(self.n+1):
qnm=rdowndata[i,1:3,position]
w_m_a[i] = qnm[0]/self.mf
tau_ma_a[i] = -1/(qnm[1])*self.mf
return w_m_a, tau_ma_a
def w_fpars_Berti(self,n):
return f_fpars[n]
def tau_qpars_Berti(self,n):
return q_fpars[n]
def mass_from_wtau(self,n,w,tau):
f1,f2,f3 = w_fpars_Berti(n)
q1,q2,q3 = tau_qpars_Berti(n)
res=(f1 + f2*(2**(-1/q3)*((-2*q1 + w*tau)/q2)**(1/q3))**f3)/w
return res
def spin_from_wtau(self,n,w,tau):
f1,f2,f3 = w_fpars_Berti(n)
q1,q2,q3 = tau_qpars_Berti(n)
res=1 - 2**(-1/q3)*((-2*q1 + w*tau)/q2)**(1/q3)
return res
def mass_from_wtau_loop(self,w,tau,l,m):
res=[None]*dim
for n in range (0,dim):
f1,f2,f3 = w_fpars_Berti(n)
q1,q2,q3 = tau_qpars_Berti(n)
res[n]=(f1 + f2*(2**(-1/q3)*((-2*q1 + w[n]*tau[n])/q2)**(1/q3))**f3)/w[n]
return res
def spin_from_wtau_loop(self,w,tau,l,m):
res=[None]*dim
for n in range (0,dim):
f1,f2,f3 = w_fpars_Berti(n)
q1,q2,q3 = tau_qpars_Berti(n)
res[n]= 1 - 2**(-1/q3)*((-2*q1 + w[n]*tau[n])/q2)**(1/q3)
return res
def rd_model_wtau(self,theta):
"""RD model parametrized with the damping time tau.
"""
assert int(len(theta)/4) == self.dim, 'Please recheck your n and parameters'
wvars = theta[ : (self.dim)]
tvars = theta[(self.dim) : 2*(self.dim)]
xvars = theta[2*(self.dim) : 3*(self.dim)]
yvars = theta[3*(self.dim) : ]
ansatz = 0
for i in range (0,self.dim):
ansatz += (xvars[i]*np.exp(1j*yvars[i]))*np.exp(-self.time/tvars[i]) * (np.cos(wvars[i]*self.time)-1j*np.sin(wvars[i]*self.time))
# -1j to agree with SXS convention
return ansatz
def rd_model_wtau_m_af(theta):
"""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 = theta[index_mass]
spin_vars = theta[index_spin]
w_m_a , tau_m_a = dict_omega[self.qnm_model](mass_vars,spin_vars,2,2)
ansatz = 0
for i in range (0,dim):
ansatz += (xvars[i]*np.exp(1j*yvars[i]))*np.exp(-timesrd_final_tsh/tau_m_a[i]) * (np.cos(w_m_a[i]*timesrd_final_tsh)-1j*np.sin(w_m_a[i]*timesrd_final_tsh))
# -1j to agree with SXS convention
return ansatz
def rd_model_wtau_fixed(theta):
"""RD model parametrized with the damping time tau and with the QNM spectrum fixd to GR.
"""
xvars = theta[ : (dim)]
yvars = theta[(dim) : 2*(dim)]
ansatz = 0
for i in range (0,dim):
ansatz += (xvars[i]*np.exp(1j*yvars[i]))*np.exp(-timesrd_final_tsh/tau[i]) * (np.cos(w[i]*timesrd_final_tsh)-1j*np.sin(w[i]*timesrd_final_tsh))
# -1j to agree with SXS convention
return ansatz
def rd_model_wq(self,theta):
"""RD model parametrized with the quality factor q.
"""
assert int(len(theta)/4) == self.dim, 'Please recheck your n and parameters'
wvars = theta[ : (self.dim)]
qvars = theta[(self.dim) : 2*(self.dim)]
xvars = theta[2*(self.dim) : 3*(self.dim)]
yvars = theta[3*(self.dim) : ]
ansatz = 0
for i in range (0,self.dim):
ansatz += (xvars[i]*np.exp(1j*yvars[i]))*np.exp(-self.time*np.pi*wvars[i]/qvars[i])*(np.cos(wvars[i]*self.time)-1j*np.sin(wvars[i]*self.time))
# -1j to agree with SXS convention
return ansatz
def rd_model_wq_fixed(self,theta):
"""RD model parametrized with the damping time tau and with the QNM spectrum fixd to GR.
"""
xvars = theta[ : (self.dim)]
yvars = theta[(self.dim) : 2*(self.dim)]
ansatz = 0
for i in range (0,self.dim):
ansatz += (xvars[i]*np.exp(1j*yvars[i]))*np.exp(-self.time/self.tau[i]) * (np.cos(self.w[i]*self.time)-1j*np.sin(self.w[i]*self.time))
# -1j to agree with SXS convention
return ansatz
def rd_model_wq_m_a(self,theta):
"""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[ : (self.dim)]
yvars = theta[(self.dim) : 2*(self.dim)]
mass_vars = theta[-2]
spin_vars = theta[-1]
w_m_a , tau_m_a = QNM_spectrum()
ansatz = 0
for i in range (0,dim):
ansatz += (xvars[i]*np.exp(1j*yvars[i]))*np.exp(-timesrd_final_tsh/tau_m_a[i]) * (np.cos(w_m_a[i]*timesrd_final_tsh)-1j*np.sin(w_m_a[i]*timesrd_final_tsh))
# -1j to agree with SXS convention
return ansatz
\ No newline at end of file
code_new/source_code/rdown_utilities.py 0 → 100644
+ 318
− 0
View file @ 68665660
# Copyright (C) 2021 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
#
# =============================================================================
#
# Module to run PE on RD data
import numpy as np
from dynesty.utils import resample_equal
from dynesty import utils as dyfunc
import os
import csv
import pandas as pd
import pickle
def posterior_samples(sampler):
"""
Returns posterior samples from nested samples and weights
given by dynsety sampler
"""
dynesty_samples = sampler.results['samples']
wt = np.exp(sampler.results['logwt'] -
sampler.results['logz'][-1])
# Make sure that sum of weights equal to 1
weights = wt/np.sum(wt)
posterior_dynesty = dyfunc.resample_equal(dynesty_samples, weights)
return posterior_dynesty
def FFT_FreqBins(times):
Len = len(times)
DeltaT = times[-1]- times[0]
dt = DeltaT/(Len-1)
dnu = 1/(Len*dt)
maxfreq = 1/(2*dt)
add = dnu/4
p = np.arange(0.0,maxfreq+add,dnu)
m = np.arange(p[-1]-(2*maxfreq)+dnu,-dnu/2+add,dnu)
res=np.concatenate((p,m))
return res
def hFromPsi4FFI(tpsi4,f0):
timecheck1=tpsi4[-2,0]-tpsi4[-1,0]
timecheck2=tpsi4[1,0]-tpsi4[0,0]
if np.abs(timecheck1-timecheck2)>=0.0001:
print("The data might not be equally sampled!!")
times,data= tpsi4[:,0],tpsi4[:,1]
freqs = FT_FreqBins(xaxis.real).real
position = np.argmax(freqs >= f0)
freqs[:position]=f0*np.ones(len(freqs[:position]))
freqs=2*np.pi*freqs
fdata=fft(data)
len(myTable)*ifft(- fdata/floor**2);
np.stack((times,data)).T
def twopoint_autocovariance(t,n):
""" It computes the two-point autocovariance function.
"""
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)])
def save_object(obj, filename):
with open(filename, 'wb') as output: # Overwrites any existing file.
pickle.dump(obj, output, pickle.HIGHEST_PROTOCOL)
def EasyMatchT(t,h1,h2,tmin,tmax):
""" It computes the time-domain match for (h1|h2) complex waveforms.
"""
pos = np.argmax(t >= (tmin));
h1red=h1[pos:];
h2red=h2[pos:];
norm1=np.sum(np.abs(h1red)**2)
norm2=np.sum(np.abs(h2red)**2)
myTable=h1red*np.conjugate(h2red)
res=((np.sum(myTable)/np.sqrt(norm1*norm2))).real
return res
def EasySNRT(t,h1,h2,tmin,tmax):
""" It computes the time-domain snr for (h1|h2) 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 FindTmaximum(y):
""" 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)
timemax=y[index,0]
return timemax
def export_logz_files(output_file,pars):
sim_num, nmax, tshift, evidence, evidence_error = pars
"""
Generate the logz.csv files you want to export the data to.
file_type must be one of this options: [corner_plot,corner_plot_extra,diagnosis,fit,post_samples,sampler_results,log_z]
"""
summary_titles=['n','id','t_shift','dlogz','dlogz_err']
if os.path.exists(output_file):
outvalues = np.array([[nmax, sim_num, tshift, evidence,evidence_error]])
else:
outvalues = np.array([summary_titles,[nmax, sim_num, tshift, evidence,evidence_error]])
with open(output_file, 'a') as file:
writer = csv.writer(file)
if (outvalues.shape)[0]>1 :
writer.writerows(outvalues)
else:
writer.writerow(outvalues[0])
return
def export_bestvals_files(best_data_file,postsamps,pars):
tshift, lenpriors, labels = pars
sigma_vars_m = np.empty(lenpriors)
sigma_vars_p = np.empty(lenpriors)
sigma_vars = np.empty(lenpriors)
sigma_vars_ml = np.empty(lenpriors)
for i in range(lenpriors):
amps_aux = postsamps[:,i]
sigma_vars_m[i] = np.quantile(amps_aux, 0.05)
sigma_vars[i] = np.quantile(amps_aux, 0.5)
sigma_vars_ml[i] = postsamps[-1,i]
sigma_vars_p[i] = np.quantile(amps_aux, 0.95)
sigma_vars_all = [sigma_vars,sigma_vars_ml,sigma_vars_m,sigma_vars_p]
sigma_vars_all=np.stack([sigma_vars,sigma_vars_ml,sigma_vars_m,sigma_vars_p], axis=0)
key =['max val','max val ml','lower bound','higher bound']
dfslist = [pd.DataFrame(np.concatenate(([tshift],sigma_vars_all[i])).reshape((-1,lenpriors+1)), columns=np.concatenate((['tshift'],labels)), index = [key[i]]) for i in range(4)]
df2 = pd.concat(dfslist)
if os.path.exists(best_data_file):
df2.to_csv(best_data_file, mode='a', header=False,index = True)
else:
df2.to_csv(best_data_file, index = True)
def define_labels(dim,model,fitnoise):
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))
if fitnoise:
noise_lab = ['noise']
labels = np.concatenate((labels,noise_lab))
return labels
def get_truths(model,pars,fitnoise):
w, tau, mf, af , npamps = pars
if model == 'w-q':
tau_val = np.pi*w*tau
truths = np.concatenate((w,tau_val,npamps))
elif model == 'w-tau':
tau_val = tau
truths = np.concatenate((w,tau_val,npamps))
elif model == 'w-tau-fixed':
truths = npamps
elif model == 'w-tau-fixed-m-af':
truths = np.concatenate((npamps,[mf],[af]))
if fitnoise:
truths = np.concatenate((truths,[1]))
return truths
def get_best_amps(pars,parser=None,nr_code=None):
nmax,model,samps_tr,half_points = pars
if model=='w-tau-fixed':
rg = (nmax+1)
elif model=='w-tau-fixed':
rg = (nmax+1)+2
else:
rg = (nmax+1)*2
if model=='w-tau-fixed-a-mf':
npamps = np.empty((nmax+1))
for i in range(0,(nmax+1)):
amps_aux = samps_tr[i+rg][half_points:-1]
npamps[i] = np.quantile(amps_aux, 0.5)
else :
npamps = np.empty((nmax+1)*2)
for i in range(0,(nmax+1)*2):
amps_aux = samps_tr[i][half_points:-1]
npamps[i] = np.quantile(amps_aux, 0.5)
if nr_code == 'Mock-data':
nm_mock = parser.get('rd-mock-parameters','nm_mock')
nm_mock = np.int(nm_mock)
amp_mock=np.empty(nm_mock+1)
ph_mock=np.empty(nm_mock+1)
for i in range(nm_mock+1):
amp_mockp = parser.get('rd-mock-parameters','amp'+str(i))
amp_mock[i] = np.float(amp_mockp)
ph_mockp=parser.get('rd-mock-parameters','phase'+str(i))
ph_mock[i] = np.float(ph_mockp)
npamps = np.concatenate((amp_mock,ph_mock))
return npamps
def convert_m_af_2_w_tau_post(res,fitnoise=False):
samples_2=res.samples
samps=f2.results.samples
if fitnoise:
fmass_spin=(samps.T)[-3:-1].T
else:
fmass_spin=(samps.T)[-2:].T
#fmass_spin=new_samples[-2:]
fmass_spin_dist=[None]*len(fmass_spin)
weight=np.exp(res.logwt - res.logz[-1])
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))
fmass_spin_dist_v2=np.asarray(fmass_spin_dist)
new_samples = dyfunc.resample_equal(fmass_spin_dist_v2, weight)
return new_samples
def save_object(obj, filename):
with open(filename, 'wb') as output: # Overwrites any existing file.
pickle.dump(obj, output, pickle.HIGHEST_PROTOCOL)
def rm_files(files):
""" rm all old files """
for i in files:
if os.path.exists(i):
os.remove(i)
\ No newline at end of file
code_new/source_code/read_data.py
+ 3
− 3
View file @ 68665660
...@@ -226,11 +226,11 @@ def create_output_files(output_folder,pars,file_type): ...@@ -226,11 +226,11 @@ def create_output_files(output_folder,pars,file_type):
""" """
if file_type=='corner_plot': if file_type=='corner_plot':
outfile = output_folder+'/Dynesty_'+str(sim_num)+'_'+model+'_nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'corner_plot.png' outfile = output_folder+'/Dynesty_'+str(sim_num)+'_'+model+'_nmax_'+str(nmax)+'_tshift_'+str(tshift)+'_'+str(npoints)+'corner_plot.png'
elif file_type=='corner_plot_extra': elif file_type=='corner_plot_extra':
outfile = output_folder+'/Dynesty_'+str(sim_num)+'_'+model+'_nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'corner_plot_extra.png' outfile = output_folder+'/Dynesty_'+str(sim_num)+'_'+model+'_nmax_'+str(nmax)+'_tshift_'+str(tshift)+'_'+str(npoints)+'corner_plot_extra.png'
elif file_type=='diagnosis': elif file_type=='diagnosis':
outfile = output_folder+'/Dynesty_diagnosis'+str(sim_num)+'_'+model+'_nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'.png' outfile = output_folder+'/Dynesty_diagnosis'+str(sim_num)+'_'+model+'_nmax_'+str(nmax)+'_tshift_'+str(tshift)+'_'+str(npoints)+'.png'
elif file_type=='fit': elif file_type=='fit':
outfile = output_folder+'/Fit_results_'+str(sim_num)+'tshift_'+str(tshift)+'_'+model+'_nmax_'+str(nmax)+'.png' outfile = output_folder+'/Fit_results_'+str(sim_num)+'tshift_'+str(tshift)+'_'+model+'_nmax_'+str(nmax)+'.png'
elif file_type=='post_samples': elif file_type=='post_samples':
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment

AEI Privacy Policy Imprint