add new files

code_new/Sumit/rdown_utilities.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
+from dynesty.utils import resample_equal
+from dynesty import utils as dyfunc
+import os
+import csv
+import pandas as pd
+
+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
\ No newline at end of file

code_new/Sumit/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