source code

code_new/source_code/RD_Fits.py

+#!/usr/bin/env python
+# coding: utf-8
+
+# In[38]:
+
+
+"""Generate ringdown templates in the time and perform parameter estimation on them.
+"""
+
+
+# In[39]:
+
+
+#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
+
+
+# In[40]:
+
+
+## Loading and running data tested with NR data
+## Loading and running data tested with Mock data
+
+
+# In[41]:
+
+
+# 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_n0_to_1_mock.ini')
+    parser.sections()
+    pass
+
+
+# In[42]:
+
+
+# Load variables from config file
+(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)
+
+
+# In[43]:
+
+
+# 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)
+
+
+# In[44]:
+
+
+# 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 ")
+
+
+# In[45]:
+
+
+# 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')
+
+
+# In[46]:
+
+
+#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)
+
+
+# In[47]:
+
+
+# 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()
+
+
+# In[48]:
+
+
+# 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
+
+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_pe = rd_pe.Ringdown_PE(rdown,datar_al,dim,priors,errors2=error_final,norm_factor=norm_factor,model=model,l_int=l_int)
+
+
+# In[49]:
+
+
+# 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)
+
+
+# In[50]:
+
+
+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)
+
+
+# In[52]:
+
+
+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)
+
+
+# In[53]:
+
+
+pars = nmax,model,samps_tr, half_points
+npamps = rd_ut.get_best_amps(pars,parser=parser,nr_code=nr_code)
+
+
+# In[54]:
+
+
+if export:
+    pars = simulation_number, nmax, tshift, evidence, evidence_error
+    rd_ut.export_logz_files(sumary_data,pars)
+
+
+# In[55]:
+
+
+labels = rd_ut.define_labels(dim,model,fitnoise)
+if export:
+    pars = tshift, len(priors), labels 
+    rd_ut.export_bestvals_files(best_data,postsamps,pars)
+
+
+# In[56]:
+
+
+w, tau = rdown.QNM_spectrum()
+pars = w, tau, mf, af, npamps 
+truths = rd_ut.get_truths(model,pars,fitnoise)
+
+
+# In[57]:
+
+
+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')
+
+
+# In[58]:
+
+
+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')
+
+
+# In[151]:
+
+
+#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')
+
+
+# In[166]:
+
+
+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)
+
+
+# In[162]:
+
+
+if export:
+    with open(samples_file,'w') as file:
+        writer = csv.writer(file)
+        writer.writerow(labels)
+        writer.writerows(samps[::downfactor])
+

code_new/source_code/read_data.py

+# 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
+import rdown_utilities as rd_ut
+import romspline
+import rdown as rd
+import h5py
+import json
+from scipy import interpolate
+from scipy.interpolate import interp1d
+
+
+def read_data(nr_code,sim_path,mf=1,af=0,parser=None,RD=True,tshift=0,tend = 100,metadata_file=None):
+    if nr_code=='SXS': 
+        gw = {}
+        gw = h5py.File(sim_path, 'r') 
+        gw_data = gw["Extrapolated_N3.dir"]["Y_l2_m2.dat"] 
+        times = gw_data[:,0]
+        
+        metadata = {}
+        with open(metadata_file) as file:
+            metadata = json.load(file)
+
+        af = metadata['remnant_dimensionless_spin'][-1]
+        mf = metadata['remnant_mass']
+
+        tmax=rd_ut.FindTmaximum(gw_data[round(len(gw_data)/2):])
+        times = times - tmax
+ 
+        if RD:
+            position = np.argmax(times >= 0)
+            gw_data = gw_data[:,1][position:]+1j*gw_data[:,2][position:]
+            times = times[times >= 0]            
+            
+    elif nr_code=='Maya': 
+        dt=0.1
+        gw = {}
+        gw = h5py.File(sim_path, 'r') 
+        gw_sxs_bbh_0305_amp = np.asarray(gw['amp_l2_m2/Y'])[6:] 
+        times_1 = np.asarray(gw['amp_l2_m2/X'])[6:] 
+        gw_sxs_bbh_0305_amp_int = romspline.ReducedOrderSpline(times_1, gw_sxs_bbh_0305_amp)
+        gw_sxs_bbh_0305_pha = np.asarray(gw['phase_l2_m2/Y'])[6:]
+        times = np.asarray(gw['phase_l2_m2/X'])[6:]
+        gw_sxs_bbh_0305_pha_int = romspline.ReducedOrderSpline(times, gw_sxs_bbh_0305_pha)
+
+        tmin=max(times_1[0],times[0])
+        tmax=min(times_1[-1],times[-1])
+        times=np.arange(tmin,tmax,dt)
+        amps=gw_sxs_bbh_0305_amp_int(times)
+        phs=gw_sxs_bbh_0305_pha_int(times)
+        gw_sxs_bbh_0305 = np.asarray([times,amps*np.cos(phs),amps*np.sin(phs)]).T
+        gw5_sxs_bbh_0305 = gw_sxs_bbh_0305
+
+        times = gw_sxs_bbh_0305[:,0]
+        tmax=rd_ut.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=rd_ut.FindTmaximum(gw5_sxs_bbh_0305[round(len(gw_sxs_bbh_0305)/2):])
+        times5 = times5 - tmax5
+            #Select the data from 0 onwards
+        position = np.argmax( times >= (t_align))
+        position5 = np.argmax(times5 >= (t_align))
+        gw_sxs_bbh_0305rd=gw_sxs_bbh_0305[position+1:]
+        gw_sxs_bbh_0305rd5=gw5_sxs_bbh_0305[position5+1:]
+        timesrd=gw_sxs_bbh_0305[position:-1][:,0][:]
+        timesrd5=gw5_sxs_bbh_0305[position5:-1][:,0][:]
+
+    elif nr_code=='LaZeV': 
+        dt=0.1
+        gw = {}
+        gw = h5py.File(simulation_path_1, 'r') 
+        gw_sxs_bbh_0305_amp = np.asarray(gw['amp_l2_m2/Y'])[6:] 
+        times_1 = np.asarray(gw['amp_l2_m2/X'])[6:] 
+        gw_sxs_bbh_0305_amp_int = romspline.ReducedOrderSpline(times_1, gw_sxs_bbh_0305_amp)
+        gw_sxs_bbh_0305_pha = np.asarray(gw['phase_l2_m2/Y'])[6:]
+        times = np.asarray(gw['phase_l2_m2/X'])[6:]
+        gw_sxs_bbh_0305_pha_int = romspline.ReducedOrderSpline(times, gw_sxs_bbh_0305_pha)
+
+        tmin=max(times_1[0],times[0])
+        tmax=min(times_1[-1],times[-1])
+        times=np.arange(tmin,tmax,dt)
+        amps=gw_sxs_bbh_0305_amp_int(times)
+        phs=gw_sxs_bbh_0305_pha_int(times)
+        gw_sxs_bbh_0305 = np.asarray([times,amps*np.cos(phs),amps*np.sin(phs)]).T
+        gw5_sxs_bbh_0305 = gw_sxs_bbh_0305
+
+        times = gw_sxs_bbh_0305[:,0]
+        tmax=rd_ut.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=rd_ut.FindTmaximum(gw5_sxs_bbh_0305[round(len(gw_sxs_bbh_0305)/2):])
+        times5 = times5 - tmax5
+
+            #Select the data from 0 onwards
+        position = np.argmax( times >= (t_align))
+        position5 = np.argmax(times5 >= (t_align))
+        gw_sxs_bbh_0305rd=gw_sxs_bbh_0305[position+1:]
+        gw_sxs_bbh_0305rd5=gw5_sxs_bbh_0305[position5+1:]
+        timesrd=gw_sxs_bbh_0305[position:-1][:,0][:]
+        timesrd5=gw5_sxs_bbh_0305[position5:-1][:,0][:]
+
+    elif nr_code=='Mock-data': 
+        times = np.arange(tshift,tend+10,0.1)
+
+        nm_mock = parser.get('rd-mock-parameters','nm_mock')
+        nm_mock = np.int(nm_mock)
+        mf = parser.get('rd-mock-parameters','mf')
+        mf = np.float(mf)
+        af = np.float(parser.get('rd-mock-parameters','af'))
+        af = np.float(af)
+
+        rdown=rd.Ringdown_Spectrum(mf,af,2,2,n=nm_mock,s=-2,time=times)
+
+        w_mock=np.empty(nm_mock+1)
+        tau_mock=np.empty(nm_mock+1)
+        amp_mock=np.empty(nm_mock+1)
+        ph_mock=np.empty(nm_mock+1)
+
+        for i in range(nm_mock+1): 
+            wp_mock = parser.get('rd-mock-parameters','w'+str(i))
+            w_mock[i] = np.float(wp_mock)
+            tp_mock=parser.get('rd-mock-parameters','tau'+str(i))
+            tau_mock[i] = np.float(tp_mock)
+            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)
+
+
+        pars = np.concatenate((w_mock,tau_mock,amp_mock,ph_mock))
+        gw_data=rdown.rd_model_wtau(pars)
+        
+    return  np.stack((times,gw_data)).T
+
+def nr_resize(data_1,data_2,tshift=0,tend=100):
+    times_1 = data_1[:,0].real
+    times_2 = data_2[:,0].real
+    data_1_re = data_1[:,1].real
+    data_1_im = data_1[:,1].imag
+    data_2_re = data_2[:,1].real
+    data_2_im = data_2[:,1].imag
+
+    gwnew_re = interpolate.interp1d(times_1, data_1_re, kind = 'cubic')
+    gwnew_im = interpolate.interp1d(times_1, data_1_im, kind = 'cubic')
+    gwnew_re5 = interpolate.interp1d(times_2, data_2_re, kind = 'cubic')
+    gwnew_im5 = interpolate.interp1d(times_2, data_2_im, kind = 'cubic')
+    
+    if times_2[-1]>= times_1[-1]: 
+        times_rd = times_1
+    else:
+        times_rd = times_2
+
+    gwdatanew_re = gwnew_re(times_rd)
+    gwdatanew_im = gwnew_im(times_rd)
+    gwdatanew_re5 = gwnew_re5(times_rd)
+    gwdatanew_im5 = gwnew_im5(times_rd)
+    gwdatanew = gwdatanew_re + 1j*gwdatanew_im
+    gwdatanew5 = gwdatanew_re5 + 1j*gwdatanew_im5
+    
+    position_in = np.argmax(times_rd >= tshift)
+    position_end = np.argmax(times_rd >= tend)
+    times_rd = times_rd[position_in:position_end]
+    gwdatanew = gwdatanew[position_in:position_end]
+    gwdatanew5 = gwdatanew5[position_in:position_end]
+
+    return(np.stack((times_rd,gwdatanew)).T,np.stack((times_rd,gwdatanew5)).T)
+
+
+def phase_align(data_1,data_2,t_align=0):
+    
+    #timesrd_final = data_1[:,0]
+    #phas = np.angle(data_1[:,1])
+    #phas = np.unwrap(phas)
+    #phas5 = np.angle(data_2[:,1])
+    #phas5 = np.unwrap(phas5)
+    #position = np.argmax(timesrd_final >= (t_align))
+    #dphase = phas5[position]-phas[position]
+    
+    #gwdatanew = data_1[:,1]*np.exp(1j*dphase)
+    
+    
+    phas = np.angle(data_1[:,1])
+    phas = np.unwrap(phas)
+    phas5 = np.angle(data_2[:,1])
+    phas5 = np.unwrap(phas5)
+    position = np.argmax(data_1[:,0] >= (0))
+    dphase = phas5[position]-phas[position]
+    gwdatanew = data_1[:,1]*np.exp(1j*dphase)
+    timesrd_final = data_1[:,0]
+    
+    return np.stack((timesrd_final,gwdatanew)).T
+
+def create_output_files(output_folder,pars,file_type):
+    sim_num, model, nmax, tshift, npoints = pars
+    
+    """
+    Generate the output 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]
+    """
+    
+    if file_type=='corner_plot':
+        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':
+        outfile = output_folder+'/Dynesty_'+str(sim_num)+'_'+model+'_nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'corner_plot_extra.png'
+    elif file_type=='diagnosis':
+        outfile = output_folder+'/Dynesty_diagnosis'+str(sim_num)+'_'+model+'_nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'.png'
+    elif file_type=='fit':
+        outfile = output_folder+'/Fit_results_'+str(sim_num)+'tshift_'+str(tshift)+'_'+model+'_nmax_'+str(nmax)+'.png'
+    elif file_type=='post_samples':
+        outfile = output_folder+'/posterior_samples-'+str(sim_num)+'tshift_'+str(tshift)+'_'+model+'_nmax_'+str(nmax)+'.csv'
+    elif file_type=='sampler_results':
+        outfile = output_folder+'/results_'+str(sim_num)+'tshift_'+str(tshift)+'_'+model+'_nmax_'+str(nmax)+'.pkl'
+    elif file_type=='log_z':
+        outfile = output_folder+'/summary'+str(sim_num)+'_'+model+'_nmax_'+str(nmax)+'.csv'
+    elif file_type=='best_vals':
+        outfile = output_folder+'/best_values_'+str(sim_num)+'_'+model+'_nmax_'+str(nmax)+'.csv'
+    else:
+        print ('Something went wrong')
+        return 
+    
+    return outfile
+
+def read_config_file(parser):
+    # Setup path and output folders
+    rootpath=parser.get('nr-paths','rootpath')
+    simulation_path_1 = parser.get('nr-paths','simulation_path_1')
+
+    if parser.get('nr-paths','simulation_path_2'):
+        simulation_path_2 = parser.get('nr-paths','simulation_path_2')
+    else:
+        simulation_path_2 = simulation_path_1
+
+    metadata_file = parser.get('nr-paths','metadata_json')
+    simulation_number = parser.get('nr-paths','simulation_number')
+    simulation_number = np.int(simulation_number)
+
+    output_folder = parser.get('output-folder','output-folder')
+    if parser.has_option('setup','export'):
+        export=eval(parser.get('setup','export'))
+    else:
+        export=True
+        
+    # Setup sampler and output options
+    overwrite = eval(parser.get('setup','overwrite'))
+    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','nb_cores'):
+        nbcores = np.int(parser.get('setup','nb_cores'))
+    else:
+        nbcores = 1
+        
+    # time shift , end and align options
+    tshift=parser.get('time-setup','tshift')
+    tshift = np.float(tshift)
+
+    tend=parser.get('time-setup','tend')
+    tend = np.float(tend)
+
+    t_align=parser.get('time-setup','t_align')
+    t_align = np.float(t_align)
+    
+    # n-tones & nlive points
+    nmax=parser.get('n-tones','nmax')
+    nmax = np.int(nmax)
+
+    npoints=parser.get('n-live-points','npoints')
+    npoints = np.int(npoints)
+    
+    # setup the RD model
+    model=parser.get('rd-model','model')
+    error_str = eval(parser.get('rd-model','error_str'))
+    fitnoise=eval(parser.get('rd-model','fit_noise'))
+    if fitnoise:
+        l_int=1
+        index_mass=-3
+        index_spin=-2
+           # prior_dim = len(priors_min)
+    else: 
+        index_mass=-2
+        index_spin=-1
+        l_int=0
+
+    if error_str:
+        error_val=np.float(parser.get('rd-model','error_val'))
+        if error_val==0:
+            error_type=''
+        else:
+            error_type=error_val
+    else:
+        error_type='False'
+        error_val =0
+
+    if nr_code == 'SXS':
+        metadata = {}
+        with open(metadata_file) as file:
+            metadata = json.load(file)
+        af = metadata['remnant_dimensionless_spin'][-1]
+        mf = metadata['remnant_mass']
+    else: 
+        mf = parser.get('rd-mock-parameters','mf')
+        mf = np.float(mf)
+        af = np.float(parser.get('rd-mock-parameters','af'))
+        af = np.float(af)
+
+    if model == 'w-q':
+        tau_var_str='q'
+    else:
+        tau_var_str='tau'
+        
+    if nr_code == 'Mock-data':   
+        nm_mock = int(parser.get('rd-mock-parameters','nm_mock'))
+    else:
+        nm_mock = None
+        
+    res = 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
+    return res
\ No newline at end of file