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

correcting rd_model_wq_m_a

parent cc60ebec
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
code_new/Sumit/RD_Fits.ipynb
+ 75
− 6
View file @ 1714269b
%% Cell type:code id: tags:
``` python
"""Generate ringdown templates in the time and perform parameter estimation on them.
"""
```
%% Output
'Generate ringdown templates in the time and perform parameter estimation on them.\n'
%% Cell type:code id: tags:
``` python
#Import relevant modules, import data and all that
import time
import numpy as np
import corner
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
from matplotlib import rc
from configparser import ConfigParser
plt.rcParams.update({'font.size': 16.5})
from multiprocessing import Pool
import random
import dynesty
from dynesty import plotting as dyplot
from dynesty.utils import resample_equal
from dynesty import utils as dyfunc
import os
import argparse
import scipy.optimize as optimization
from scipy.optimize import minimize
import rdown as rd
import rdown_pe as rd_pe
import rdown_utilities as rd_ut
import read_data as rdata
```
%% Cell type:code id: tags:
``` python
## Loading and running data tested with NR data
## Loading and running data tested with Mock data
```
%% Cell type:code id: tags:
``` python
# 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")
result = parser.parse_args()
config_file=result.config_file
parser = ConfigParser()
parser.read(config_file)
parser.sections()
except SystemExit:
parser = ConfigParser()
parser.read('config_fixed_n6_m_af.ini')
parser.sections()
pass
```
%% Output
usage: ipykernel_launcher.py [-h] [-c CONFIG_FILE]
ipykernel_launcher.py: error: unrecognized arguments: -f /work/francisco.jimenez/.local/share/jupyter/runtime/kernel-433f529b-f75c-4320-8450-e97305857525.json
ipykernel_launcher.py: error: unrecognized arguments: -f /work/francisco.jimenez/.local/share/jupyter/runtime/kernel-b3311339-24b7-4dfe-a2d4-890fff02f1cb.json
%% Cell type:code id: tags:
``` python
# Load variables from config file
(rootpath, simulation_path_1,simulation_path_2, metadata_file , simulation_number, output_folder,
export, overwrite, sampler,nr_code, nbcores,tshift,tend,t_align,
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)
```
%% Cell type:code id: tags:
``` python
# Show configuration options
dim = nmax+1
ndim = 4*dim
numbins = 32 #corner plot parameter - how many bins you want
print('model:',model)
print('nmax:',nmax)
print('nm_mock:',nm_mock)
print('tshift:',tshift)
print('error:', error_str)
print('error value:',error_val)
print('export:',export)
print('nr code:',nr_code)
print('fit noise:',fitnoise)
```
%% Output
model: w-tau-fixed-m-af
nmax: 6
nm_mock: None
tshift: 0.0
error: True
error value: 0.02
export: True
nr code: SXS
fit noise: True
%% Cell type:code id: tags:
``` python
# Create output directories
if not os.path.exists(output_folder):
os.mkdir(output_folder)
print("Directory " , output_folder , " Created ")
if nr_code == 'Mock-data':
nm_mock_str = 'rec_with'+parser.get('rd-mock-parameters','nm_mock')+'_'
else:
nm_mock_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))
else:
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):
os.mkdir(output_folder_1)
print("Directory " , output_folder_1 , " Created ")
```
%% Cell type:code id: tags:
``` python
# Define output files
pars = [simulation_number,model,nmax,tshift,npoints]
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')
diagnosis_plot = rdata.create_output_files(output_folder_1,pars,'diagnosis')
fit_plot = rdata.create_output_files(output_folder_1,pars,'fit')
samples_file = rdata.create_output_files(output_folder_1,pars,'post_samples')
results_file = rdata.create_output_files(output_folder_1,pars,'sampler_results')
sumary_data = rdata.create_output_files(output_folder_1,pars,'log_z')
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:
``` python
#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_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)
times_rd = data_r[:,0]
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, 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, 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:
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)
plt.plot(times_rd, error_est, "g", alpha=0.3, lw=2, label='error')
plt.legend()
mismatch=1-rd_ut.EasyMatchT(times_rd,data_r[:,1],data_lr[:,1],tshift,tend)
error=np.sqrt(2*mismatch)
print('error estimate:',error)
print('mismatch:', mismatch)
print('snr:', rd_ut.EasySNRT(times_rd,data_r[:,1],data_lr[:,1],tshift,tend)/error**2)
```
%% 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
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
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
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
return array(a, dtype, copy=False, order=order)
error estimate: 0.2005842760402416
mismatch: 0.020117025897293916
snr: 228.11029832920406
%% Cell type:code id: tags:
``` python
# Phase alignement
if parser.has_option('rd-model','phase_alignment'):
phase_alignment=eval(parser.get('rd-model','phase_alignment'))
else:
phase_alignment=False
if phase_alignment:
datar_al = rdata.phase_align(data_r,data_lr)
gwdatanew5 = data_lr[:,1]
gwdatanew = datar_al[:,1]
timesrd_final = datar_al[:,0]
mismatch=1-rd_ut.EasyMatchT(timesrd_final,gwdatanew,gwdatanew5,tshift,tend)
error=np.sqrt(2*mismatch)
print('error estimate:',error)
print('mismatch:', mismatch)
print('snr:', rd_ut.EasySNRT(timesrd_final,gwdatanew,gwdatanew5,tshift,tend)/error)
if error_str:
error = np.sqrt(gwdatanew*gwdatanew-2*gwdatanew*gwdatanew5+gwdatanew5*gwdatanew5)
error_est=np.sqrt(error.imag**2+error.real**2)
else :
error = 1
else:
datar_al = data_r
timesrd_final = datar_al[:,0]
#Test the new interpolated data
if error_str and error_val==0:
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, 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.legend()
```
%% Output
error estimate: 0.007775922554703448
mismatch: 3.0232485788372898e-05
snr: 1192.3206504219345
%% Cell type:code id: tags:
``` python
# 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:
error_final = error_est
norm_factor = 100*len(error_final)/2*np.log(2*np.pi)
elif error_str and error_val!=0:
datar_al[:,1]+=random.uniform(0, error_val)
datar_al[:,1]+=1j*random.uniform(0, error_val)
error_tsh = error_val
error_final=(error_tsh.real**2+error_tsh.imag**2)
norm_factor = 0
else:
error_tsh=1
error_final=(error_tsh.real**2+error_tsh.imag**2)
norm_factor = 0
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))])
else:
qnm_model='qnm'
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, rdowndata = rdowndata)
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:
``` python
rdown.QNM_spectrum(mf, af)
```
%% Output
(array([0.55578191, 0.54351639, 0.52079511, 0.49045859, 0.45795517,
0.43991139, 0.43963542]),
array([11.74090006, 3.88312743, 2.29980777, 1.62209641, 1.25654938,
1.03071732, 0.86348675]))
%% Cell type:code id: tags:
``` python
# Get a first estimate by trying to fit the data.
nll = lambda *args: -rdown_pe.log_likelihood(*args)
if model == 'w-tau-fixed-m-af':
if fitnoise:
initial = np.concatenate((np.ones(2*dim),[0.8,0.9,1]))
soln = minimize(nll, initial,bounds=priors)
vars_ml=soln.x
else:
initial = np.concatenate((np.ones(2*dim),[0.8,0.9]))
soln = minimize(nll, initial,bounds=priors)
vars_ml=soln.x
elif model == 'w-tau-fixed':
if fitnoise:
initial = np.concatenate((np.ones(2*dim),[0.2]))
soln = minimize(nll, initial,bounds=priors)
vars_ml=soln.x
else:
initial = np.ones(2*dim)
soln = minimize(nll, initial,bounds=priors)
vars_ml=soln.x
else:
if fitnoise:
initial = np.concatenate((np.ones(ndim),[1]))
soln = minimize(nll, initial,bounds=priors)
vars_ml=soln.x
else:
initial = np.ones(ndim)
soln = minimize(nll, initial,bounds=priors)
vars_ml=soln.x
print("best fit pars from fit: ",vars_ml)
```
%% Output
best fit pars from fit: [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 0.8 0.9 1. ]
/work/francisco.jimenez/sio/git/rdstackingproject/code_new/rdown.py:81: RuntimeWarning: divide by zero encountered in double_scalars
w_m_a[i] = qnm[0]/mass
/work/francisco.jimenez/sio/git/rdstackingproject/code_new/rdown.py:150: RuntimeWarning: divide by zero encountered in true_divide
ansatz += (xvars[i]*np.exp(1j*yvars[i]))*np.exp(-self.time/tau_m_a[i]) * (np.cos(w_m_a[i]*self.time)-1j*np.sin(w_m_a[i]*self.time))
/work/francisco.jimenez/sio/git/rdstackingproject/code_new/rdown.py:150: RuntimeWarning: invalid value encountered in true_divide
ansatz += (xvars[i]*np.exp(1j*yvars[i]))*np.exp(-self.time/tau_m_a[i]) * (np.cos(w_m_a[i]*self.time)-1j*np.sin(w_m_a[i]*self.time))
/work/francisco.jimenez/sio/git/rdstackingproject/code_new/rdown.py:150: RuntimeWarning: invalid value encountered in multiply
ansatz += (xvars[i]*np.exp(1j*yvars[i]))*np.exp(-self.time/tau_m_a[i]) * (np.cos(w_m_a[i]*self.time)-1j*np.sin(w_m_a[i]*self.time))
/work/francisco.jimenez/venv/lib/python3.7/site-packages/scipy/optimize/_numdiff.py:497: RuntimeWarning: invalid value encountered in subtract
df = fun(x) - f0
%% Cell type:code id: tags:
``` python
mypool = Pool(nbcores)
mypool.size = nbcores
start = time.process_time()
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'):
dlogz=np.float(parser.get('setup','dlogz'))
f2.run_nested(dlogz=dlogz,print_progress=False)
else:
f2.run_nested(print_progress=False)
print(time.process_time() - start)
```
%% Output
Process ForkPoolWorker-2:
Process ForkPoolWorker-6:
Process ForkPoolWorker-5:
Process ForkPoolWorker-1:
Process ForkPoolWorker-7:
Process ForkPoolWorker-8:
Process ForkPoolWorker-3:
Process ForkPoolWorker-4:
Traceback (most recent call last):
File "/work/francisco.jimenez/local/python-3.7.7/lib/python3.7/multiprocessing/process.py", line 297, in _bootstrap
self.run()
Traceback (most recent call last):
Traceback (most recent call last):
Traceback (most recent call last):
Traceback (most recent call last):
%% Cell type:code id: tags:
``` python
res = f2.results
res.samples_u.shape
res.summary()
samps=f2.results.samples
postsamps = rd_ut.posterior_samples(f2)
samps_tr=np.transpose(samps)
half_points=int(round((len(samps_tr[0])/1.25)))
evidence = res.logz[-1]
evidence_error = res.logzerr[-1]
if export:
rd_ut.save_object(res, results_file)
```
%% Output
Summary
=======
nlive: 1000
niter: 46474
ncall: 1106811
eff(%): 4.289
logz: -959328.014 +/- 0.280
%% Cell type:code id: tags:
``` python
pars = nmax,model,samps_tr, half_points
npamps = rd_ut.get_best_amps(pars,parser=parser,nr_code=nr_code)
```
%% Cell type:code id: tags:
``` python
if export:
pars = simulation_number, nmax, tshift, evidence, evidence_error
rd_ut.export_logz_files(sumary_data,pars)
```
%% Cell type:code id: tags:
``` python
labels = rd_ut.define_labels(dim,model,fitnoise)
if export:
pars = tshift, len(priors), labels
rd_ut.export_bestvals_files(best_data,postsamps,pars)
```
%% Cell type:code id: tags:
``` python
w, tau = rdown.QNM_spectrum()
w, tau = rdown.QNM_spectrum(mf, af)
pars = w, tau, mf, af, npamps
truths = rd_ut.get_truths(model,pars,fitnoise)
```
%% Cell type:code id: tags:
``` 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')
plt.show()
if export:
fg.savefig(corner_plot, format = 'png', bbox_inches = 'tight')
```
%% Output
%% Cell type:code id: tags:
``` python
from importlib import reload
reload(rd_ut)
if model == 'w-tau-fixed-m-af' and export == True:
truths=np.concatenate((w,tau))
labels_mf = np.concatenate((w_lab,tau_lab))
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.savefig(corner_plot_extra, format = 'png', bbox_inches = 'tight')
```
%% Cell type:code id: tags:
``` python
#lnz_truth = ndim * -np.log(2 * 10.) # analytic evidence solution
fig, axes = dyplot.runplot(res)
fig.tight_layout()
if export:
fig.savefig(diagnosis_plot, format = 'png', dpi = 384, bbox_inches = 'tight')
```
%% Output
%% Cell type:code id: tags:
``` python
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}
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}$')
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_down = list(samples_1sigma)[::downfactor]
for sample in samples_1sigma_down:
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.legend()
plt.xlabel('t')
plt.ylabel('h')
plt.show()
figband.savefig(fit_plot)
```
%% Output
%% Cell type:code id: tags:
``` python
if export:
with open(samples_file,'w') as file:
writer = csv.writer(file)
writer.writerow(labels)
writer.writerows(samps[::downfactor])
```
......
code_new/Sumit/rdown.py
+ 8
− 8
View file @ 1714269b
......@@ -59,12 +59,12 @@ class Ringdown_Spectrum:
self.tau=-1/(np.imag(omegas_new))*self.mf
def QNM_spectrum(self):
def QNM_spectrum(self,mass,spin):
""" 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
omegas_new=np.asarray([self.grav_220[i](a=spin)[0] for i in range (0,self.n+1)])
w_m_a = (np.real(omegas_new))/mass
tau_m_a=-1/(np.imag(omegas_new))*mass
return (w_m_a, tau_m_a)
......@@ -159,7 +159,7 @@ class Ringdown_Spectrum:
ansatz = 0
for i in range (0,self.dim):
ansatz += (xvars[i]*np.exp(1j*yvars[i]))*np.exp(-self.time/tau[i]) * (np.cos(w[i]*self.time)-1j*np.sin(w[i]*self.time))
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
......@@ -200,10 +200,10 @@ class Ringdown_Spectrum:
mass_vars = theta[-2]
spin_vars = theta[-1]
w_m_a , tau_m_a = QNM_spectrum()
w_m_a , tau_m_a = self.dict_omega[self.qnm_model](mass_vars,spin_vars)
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))
for i in range (0,self.dim):
ansatz += (xvars[i]*np.exp(1j*yvars[i]))*np.exp(-self.time/tau_m_a[i]) * (np.cos(w_m_a[i]*self.time)-1j*np.sin(w_m_a[i]*self.time))
# -1j to agree with SXS convention
return ansatz
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment

AEI Privacy Policy Imprint