Commit f0a283db authored by Marc Arene's avatar Marc Arene
Browse files

Using a config.ini file

parent 4ac4e73f
This diff is collapsed.
python BNS_HMC_9D.py --n_traj_HMC_tot=801 --n_traj_fit=800 --inj_file=../examples/GW170817.ini --psd=1 --param_indices=45 --stdout_file
python BNS_HMC_9D.py --n_traj_HMC_tot=201 --n_traj_fit=200 --inj_file=../examples/GW170817.ini --psd=1 --param_indices=03 --stdout_file
python BNS_HMC_9D.py --n_traj_hmc_tot=801 --n_traj_fit=800 --inj_file=../examples/GW170817.ini --psd=1 --param_indices=45 --stdout_file
python BNS_HMC_9D.py --n_traj_hmc_tot=201 --n_traj_fit=200 --inj_file=../examples/GW170817.ini --psd=1 --param_indices=03 --stdout_file
[analysis]
ifos = H1,L1,V1
approx = IMRPhenomPv2
minimum_frequency = 30.0
maximum_frequency = 2048.0
reference_frequency = 20.0
roq = True
roq_b_matrix_directory = /Users/marcarene/Documents/APC/Code/ROQ_data/IMRPhenomPv2/
psd = 1
param_indices = 0,1,2,3,4,5,6,7,8
param_offsets = 7,0,7,0,7,7,7,7,7
dlogL = dlogL1
[hmc]
n_traj_hmc_tot = 1501
n_traj_fit = 1500
n_phase1_for_this_run = 1500
n_fit_1 = 2000
n_fit_2 = 200
length_num_traj = 200
epsilon0 = 5
[priors]
# Priors
# For all parameters known to lalinference, the min and max default prior can be overwritten with
# parname-min=MIN
# parname-max=MAX
# The starting point for the MCMC chain(s) can be specified with
# parname=START
# Parameters can be fixed to some value with
# fix-parname=FIXVALUE
# currently known parameters, together with default [min-max] are (radians for angle, Mpc for distance, Msun for masses)
# time Waveform time [trigtime-0.1-trigtime+0.1]
# chirpmass Chirpmass [1,15.3]
# eta Symmetric massratio (needs --use-eta) [0,0.0312]
# q Asymmetric massratio (a.k.a. q=m2/m1 with m1>m2) [0.033,1]
# phase Coalescence phase [0,2Pi]
# costheta_jn Cosine of angle between J and line of sight [-1,1]
# logdistance Log Distance [log(1),log(2000)]
# rightascension Rightascension [0,2Pi]
# declination Declination [-Pi/2,Pi/2]
# polarisation Polarisation angle [0,Pi]
# Spin Parameters:
# a_spin1 Spin1 magnitude [-1,1] for aligned, [0,1] for precessing
# a_spin2 Spin2 magnitude [-1,1] for aligned, [0,1] for precessing
# tilt_spin1 Angle between spin1 and orbital angular momentum [0,Pi]
# tilt_spin2 Angle between spin2 and orbital angular momentum [0, Pi]
# phi_12 Difference between spins' azimuthal angles [0,2Pi]
# phi_jl Difference between total and orbital angular momentum azimuthal angles [0,2Pi]
# Tidal parameters (requires tidal= or tidalT=):
# lambda1 lambda1 [0,3000]
# lambda2 lambda2 [0,3000]
# lambdaT lambdaT [0,3000]
# dLambdaT dLambdaT [-500,500]
# Equation of State parameters (requires 4PolyEOS=):
# logp1 logp1 [33.6,35.4]
# gamma1 gamma1 [2.0,4.5]
# gamma2 gamma2 [1.1,4.5]
# gamma3 gamma3 [1.1,4.5]
# Settings for allowed component masses in Solar Masses, with default values
# comp-max=30.0
# comp-min=1.0
# Allowed total masses in Msun (note, used together with component masses, mc,q,eta priors may lead to inconsistencies. Be careful!)
# mtotal-max=35
# mtotal-min=2
# Setting time prior [seconds]
# dt=0.1
......@@ -106,10 +106,9 @@ def ini_file_to_dict(file_path='../examples/GW170817.ini'):
# Set the duration and sampling frequency of the data segment that we're going
# to inject the signal into. For the
# TaylorF2 waveform, we cut the signal close to the isco frequency
minimum_frequency=30.0
# minimum_frequency=30.0
m1 = injection_parameters['mass_1']
m2 = injection_parameters['mass_2']
minimum_frequency = conf.minimum_frequency
tc_3p5PN = ComputeChirpTime3p5PN(minimum_frequency, m1, m2)
# We integrate the signal up to the frequency of the "Innermost stable circular orbit (ISCO)"
R_isco = 6. # Orbital separation at ISCO, in geometric units. 6M for PN ISCO; 2.8M for EOB
......@@ -119,13 +118,12 @@ def ini_file_to_dict(file_path='../examples/GW170817.ini'):
# approximant = 'SpinTaylorF2'
# approximant = 'IMRPhenomD'
# approximant = 'IMRPhenomPv2'
approximant = 'IMRPhenomPv2_ROQ'
# approximant = 'IMRPhenomPv2_NRTidal'
f_ISCO = Compute_fISCO(R_isco, m1, m2)
f_ISCO_run = Compute_fISCO(R_isco, m1_min, m2_min)
# f_ISCO = Compute_LAL_fISCO(m1 * LAL_MSUN_SI, m2 * LAL_MSUN_SI)
# f_ISCO_run = Compute_LAL_fISCO(m1_min * LAL_MSUN_SI, m2_min * LAL_MSUN_SI)
if approximant == 'TaylorF2':
if conf.approximant == 'TaylorF2':
f_high = f_ISCO
f_high_run = f_ISCO_run
maximum_frequency_injected_waveform = 0
......@@ -135,12 +133,11 @@ def ini_file_to_dict(file_path='../examples/GW170817.ini'):
sampling_frequency = 2 * f_high_run
duration = tc_3p5PN + 2
reference_frequency = 0
elif approximant == 'IMRPhenomPv2_ROQ':
approximant = 'IMRPhenomPv2'
elif conf.approximant == 'IMRPhenomPv2' and conf.roq:
# cf:
# - https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_i_m_r_phenom__c.html#gad3f98acfe9527259a7f73a8ef69a2f7b
# - line 103 of https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/_l_a_l_sim_i_m_r_phenom_d_8c_source.html#l00103
roq_directory = '/Users/marcarene/Documents/APC/Code/ROQ_data/IMRPhenomPv2/'
roq_directory = conf.roq_b_matrix_directory
params_best_basis, scale_factor = roqu.get_best_segment_params_from_chirp_mass(injection_parameters['chirp_mass'], roq_directory=roq_directory)
injection_parameters['roq'] = {}
......@@ -152,14 +149,14 @@ def ini_file_to_dict(file_path='../examples/GW170817.ini'):
minimum_frequency = rescaled_params['flow']
maximum_frequency_injected_waveform = rescaled_params['fhigh']
# maximum_frequency_injected_waveform = min(2048, rescaled_params['fhigh'])
maximum_frequency_ifo = min(2048, rescaled_params['fhigh'])
maximum_frequency_ifo = min(conf.maximum_frequency, rescaled_params['fhigh'])
maximum_frequency_search_waveform = maximum_frequency_injected_waveform
# maximum_frequency_search_waveform = maximum_frequency_ifo
# maximum_frequency_ifo = rescaled_params['fhigh']
sampling_frequency = 2 * max(2048, rescaled_params['fhigh'])
sampling_frequency = 2 * max(conf.maximum_frequency, rescaled_params['fhigh'])
duration = rescaled_params['seglen']
reference_frequency = 20 * scale_factor
reference_frequency = conf.reference_frequency * scale_factor
injection_parameters['roq']['params'] = params_best_basis
......@@ -185,13 +182,13 @@ def ini_file_to_dict(file_path='../examples/GW170817.ini'):
f_high = f_CUT / M_sec
f_high_run = f_CUT / M_min_sec
maximum_frequency_ifo = min(2048, f_high_run)
maximum_frequency_ifo = min(conf.maximum_frequency, f_high_run)
maximum_frequency_injected_waveform = maximum_frequency_ifo
maximum_frequency_search_waveform = maximum_frequency_ifo
sampling_frequency = 2 * maximum_frequency_ifo
duration = 64
reference_frequency = 20
reference_frequency = conf.reference_frequency
# import IPython; IPython.embed();
# maximum_frequency_ifo = min(4096 * 4, f_high_run)
# maximum_frequency_ifo = max(4096, f_high)
......@@ -211,16 +208,31 @@ def ini_file_to_dict(file_path='../examples/GW170817.ini'):
start_time = start_time,
f_high = f_high,
f_high_run = f_high_run,
approximant = approximant,
approximant = conf.approximant,
reference_frequency = reference_frequency
)
injection_parameters['meta'] = meta_params
conf.waveform_arguments = {}
conf.waveform_arguments['minimum_frequency'] = minimum_frequency
conf.waveform_arguments['waveform_approximant'] = conf.approximant
conf.waveform_arguments['reference_frequency'] = reference_frequency
conf.waveform_arguments['maximum_frequency'] = None
conf.minimum_frequency = minimum_frequency
conf.approximant = approximant
conf.maximum_frequency_injected_waveform = maximum_frequency_injected_waveform
conf.maximum_frequency_search_waveform = maximum_frequency_search_waveform
conf.maximum_frequency_ifo = maximum_frequency_ifo
conf.reference_frequency = reference_frequency
conf.sampling_frequency = sampling_frequency
conf.duration = duration
conf.tc_3p5PN = tc_3p5PN
conf.start_time = start_time
conf.f_high = f_high
conf.f_high_run = f_high_run
# These two lines have been added to be able to plot the injected value on phi_c on the corner plots. Otherwise the `phi_c` and `tc_3p5PN` keys are not used in the code.
injection_parameters['phi_c'] = (injection_parameters['phase'] * 2) % (2*np.pi)
......
......@@ -13,11 +13,14 @@ import Library.config as conf
def get_source_frame_polarizations(parameters, minimum_frequency, interferometers):
# breakpoint()
# Fixed arguments passed into the source model. The analysis starts at 40 Hz.
if 'meta' in parameters.keys():
# waveform_arguments = dict(waveform_approximant=parameters['meta']['approximant'], reference_frequency=parameters['meta']['reference_frequency'], minimum_frequency=minimum_frequency, maximum_frequency=0)
waveform_arguments = dict(waveform_approximant=parameters['meta']['approximant'], reference_frequency=parameters['meta']['reference_frequency'], minimum_frequency=minimum_frequency, maximum_frequency=parameters['meta']['maximum_frequency_generated_waveform'])
else:
waveform_arguments = dict(waveform_approximant=conf.approximant, reference_frequency=conf.reference_frequency, minimum_frequency=conf.minimum_frequency, maximum_frequency=conf.maximum_frequency_ifo)
# if 'meta' in parameters.keys():
# # waveform_arguments = dict(waveform_approximant=parameters['meta']['approximant'], reference_frequency=parameters['meta']['reference_frequency'], minimum_frequency=minimum_frequency, maximum_frequency=0)
# waveform_arguments = dict(waveform_approximant=parameters['meta']['approximant'], reference_frequency=parameters['meta']['reference_frequency'], minimum_frequency=minimum_frequency, maximum_frequency=parameters['meta']['maximum_frequency_generated_waveform'])
# else:
# waveform_arguments = dict(waveform_approximant=conf.approximant, reference_frequency=conf.reference_frequency, minimum_frequency=conf.minimum_frequency, maximum_frequency=conf.maximum_frequency_ifo)
waveform_arguments = conf.waveform_arguments
# Create the waveform_generator using a LAL Binary Neutron Star source function
waveform_generator = bilby.gw.WaveformGenerator(
......
......@@ -25,16 +25,16 @@ class LikelihoodGradient(object):
self.likelihood.parameters_copy = self.likelihood.parameters.copy()
for key in self.search_parameters_keys:
self.likelihood.parameters[key] -= self.offset
self.likelihood.parameters[key] = self.likelihood.parameters_copy[key] - self.offset
log_likelihood_minus = self.likelihood.log_likelihood_ratio()
# Reset the parameters to their original value
self.likelihood.parameters = self.likelihood.parameters_copy.copy()
self.likelihood.parameters[key] += self.offset
self.likelihood.parameters[key] = self.likelihood.parameters_copy[key] + self.offset
log_likelihood_plus = self.likelihood.log_likelihood_ratio()
self.log_likelihood_gradient[key] = (log_likelihood_plus - log_likelihood_minus) / (2 * self.offset)
self.likelihood.parameters = self.likelihood.parameters_copy.copy()
return self.log_likelihood_gradient
......
......@@ -108,10 +108,10 @@ def plot_num_and_ana_gradients(dlogL_num, dlogL_ana, opts_dict=None, save=True,
if opts_dict is None:
n_traj_fit = '??'
lengthTnum = '??'
length_num_traj = '??'
else:
n_traj_fit = opts_dict['n_traj_fit']
lengthTnum = opts_dict['lengthTnum']
length_num_traj = opts_dict['length_num_traj']
steps = np.arange(len(dlogL_num[0]))
......@@ -139,7 +139,7 @@ def plot_num_and_ana_gradients(dlogL_num, dlogL_ana, opts_dict=None, save=True,
axs[i].legend()
fig.suptitle(title + " - Fits done after {} trajectories of length {} in phase 1".format(n_traj_fit, lengthTnum))
fig.suptitle(title + " - Fits done after {} trajectories of length {} in phase 1".format(n_traj_fit, length_num_traj))
fig.tight_layout()
fig.subplots_adjust(top=0.95)
......@@ -179,12 +179,12 @@ def plot_num_vs_ana_positions(q_pos_num, q_pos_ana, opts_dict=None, save=True, s
if opts_dict is None:
n_traj_fit = '??'
lengthTnum = '??'
length_num_traj = '??'
else:
n_traj_fit = opts_dict['n_traj_fit']
lengthTnum = opts_dict['lengthTnum']
length_num_traj = opts_dict['length_num_traj']
fig.suptitle("Position comparison for fits done after {} trajectories of length {} in phase 1".format(n_traj_fit, lengthTnum))
fig.suptitle("Position comparison for fits done after {} trajectories of length {} in phase 1".format(n_traj_fit, length_num_traj))
fig.tight_layout()
fig.subplots_adjust(top=0.95)
......
......@@ -133,9 +133,15 @@ def config_parser_to_dict(config):
dictionary[section] = {}
for option in config.options(section):
try:
dictionary[section][option] = config.getfloat(section, option)
dictionary[section][option] = config.getint(section, option)
except ValueError:
dictionary[section][option] = config.get(section, option)
try:
dictionary[section][option] = config.getboolean(section, option)
except ValueError:
try:
dictionary[section][option] = config.getfloat(section, option)
except ValueError:
dictionary[section][option] = config.get(section, option)
return dictionary
......
def get_analytical_timing_for_dlogL(t_h_source=9, t_get_detector_response_geocent=1.3, t_get_fd_time_translation=4, t_noise_weighted_inner_product=0.7, t_dh=0.7, nb_of_ifos=3, n_traj_phase1=1000, lengthTnum=200):
def get_analytical_timing_for_dlogL(t_h_source=9, t_get_detector_response_geocent=1.3, t_get_fd_time_translation=4, t_noise_weighted_inner_product=0.7, t_dh=0.7, nb_of_ifos=3, n_traj_phase1=1000, length_num_traj=200):
"""
Analytical estimation in millisecond of the time it should take to compute one numerical gradient `dlogL` for the 9 parameters. This is based on the following description of what is implemented in bilby_waveform.dlogL_ThreeDetectors():
......@@ -37,11 +37,11 @@ def get_analytical_timing_for_dlogL(t_h_source=9, t_get_detector_response_geocen
dlogL = h_dh + s_dh_minus_h_dh
one_traj = lengthTnum * dlogL / n_traj_phase1 # in seconds
one_traj = length_num_traj * dlogL / n_traj_phase1 # in seconds
phase_1 = n_traj_phase1 * one_traj # in seconds
print("dlogL = {:.0f}ms".format(dlogL))
print("1 trajectory of {} steps = {:.1f}s".format(lengthTnum, one_traj))
print("1 trajectory of {} steps = {:.1f}s".format(length_num_traj, one_traj))
print("Phase 1 with {} trajectories = {:.1f}hrs".format(n_traj_phase1, phase_1/3600))
# import IPython; IPython.embed()
......
......@@ -449,7 +449,7 @@
"Timing results:\n",
" 18.7 sec = total time\n",
" 2.83 sec = initial gradient of the logLikelihood\n",
" 14.2 sec = looping over 1 trajectory, ie 1 NumericalGradientLeapfrog_ThreeDetectors() run with lengthTnum=5, ie 5 dlogL computations, hence 17*5=85 waveform generations\n",
" 14.2 sec = looping over 1 trajectory, ie 1 NumericalGradientLeapfrog_ThreeDetectors() run with length_num_traj=5, ie 5 dlogL computations, hence 17*5=85 waveform generations\n",
" 2.84 sec = 1 dlogL computation + other stuff inside the loop\n",
" 0.167 sec = 1 effective waveform compututation\n",
"```\n",
......@@ -460,7 +460,7 @@
"Timing results:\n",
" 3.89 sec = total time\n",
" 0.378 sec = initial gradient of the logLikelihood\n",
" 1.69 sec = looping over 1 trajectory, ie 1 NumericalGradientLeapfrog_ThreeDetectors() run with lengthTnum=5, ie 5 dlogL computations, hence 17*5=85 waveform generations\n",
" 1.69 sec = looping over 1 trajectory, ie 1 NumericalGradientLeapfrog_ThreeDetectors() run with length_num_traj=5, ie 5 dlogL computations, hence 17*5=85 waveform generations\n",
" 0.338 sec = 1 dlogL computation + other stuff inside the loop\n",
" 0.020 sec = 1 effective waveform compututation\n",
" \n",
......@@ -468,7 +468,7 @@
"Timing results:\n",
" 34.7 sec = total time\n",
" 0.347 sec = initial gradient of the logLikelihood\n",
" 32.5 sec = looping over 1 trajectory, ie 1 NumericalGradientLeapfrog_ThreeDetectors() run with lengthTnum=100, ie 100 dlogL computations, hence 17*100=1700 waveform generations\n",
" 32.5 sec = looping over 1 trajectory, ie 1 NumericalGradientLeapfrog_ThreeDetectors() run with length_num_traj=100, ie 100 dlogL computations, hence 17*100=1700 waveform generations\n",
" 0.325 sec = 1 dlogL computation + other stuff inside the loop\n",
" 0.020 sec = 1 effective waveform compututation\n",
"\n",
......@@ -487,7 +487,7 @@
"Timing results:\n",
"2.8158 sec = total time\n",
"0.0513 sec = initial gradient of the logLikelihood\n",
"2.4537 sec = looping over 1 trajectory, ie 1 NumericalGradientLeapfrog_ThreeDetectors() run with lengthTnum=100, ie 100 dlogL computations, hence 17*100=1700 waveform generations\n",
"2.4537 sec = looping over 1 trajectory, ie 1 NumericalGradientLeapfrog_ThreeDetectors() run with length_num_traj=100, ie 100 dlogL computations, hence 17*100=1700 waveform generations\n",
"0.0245 sec = 1 dlogL computation + other stuff inside the loop\n",
"0.0025 sec = 1 effective waveform compututation: doesn't make sens compared to the 9ms measured!\n",
"\n",
......@@ -508,7 +508,7 @@
"Timing results:\n",
" 5.54 sec = total time\n",
" 0.075 sec = initial gradient of the logLikelihood\n",
" 3.81 sec = looping over 1 trajectory, ie 1 NumericalGradientLeapfrog_ThreeDetectors() run with lengthTnum=100, ie 100 dlogL computations\n",
" 3.81 sec = looping over 1 trajectory, ie 1 NumericalGradientLeapfrog_ThreeDetectors() run with length_num_traj=100, ie 100 dlogL computations\n",
" 0.038 sec = 1 dlogL computation + other stuff inside the loop\n",
" 0.004 sec = 1 effective waveform compututation: doesn't make sens compared to the 9ms measured!\n",
" \n",
......
......@@ -71,7 +71,7 @@ except FileNotFoundError:
# q_pos = pt_fit_traj[param_indices]
# dlogL = dlogL_all[param_indices]
lengthT = conf_dict['lengthTnum']
lengthT = conf_dict['length_num_traj']
pt_fit_traj = pt_fit_traj[:, trajectory_index * lengthT: (trajectory_index +1) * lengthT]
dlogL = dlogL[:, trajectory_index * lengthT: (trajectory_index +1) * lengthT]
H_p_logL = H_p_logL[:, trajectory_index * lengthT: (trajectory_index +1) * lengthT]
......
[analysis]
ifos=['H1','L1','V1']
approx=IMRPhenomPv2
minimum_frequency=30
roq=True
roq_b_matrix_directory = /Users/marcarene/Documents/APC/Code/ROQ_data/IMRPhenomPv2/
[hmc]
n_traj_hmc_tot=1e5
n_traj_fit=1500
n_phase1_for_this_run=1500
n_fit_1=2000
n_fit_2=200
length_num_traj=200
epsilon0=5
[priors]
# Priors
# For all parameters known to lalinference, the min and max default prior can be overwritten with
# parname-min=MIN
# parname-max=MAX
# The starting point for the MCMC chain(s) can be specified with
# parname=START
# Parameters can be fixed to some value with
# fix-parname=FIXVALUE
# currently known parameters, together with default [min-max] are (radians for angle, Mpc for distance, Msun for masses)
# time Waveform time [trigtime-0.1-trigtime+0.1]
# chirpmass Chirpmass [1,15.3]
# eta Symmetric massratio (needs --use-eta) [0,0.0312]
# q Asymmetric massratio (a.k.a. q=m2/m1 with m1>m2) [0.033,1]
# phase Coalescence phase [0,2Pi]
# costheta_jn Cosine of angle between J and line of sight [-1,1]
# logdistance Log Distance [log(1),log(2000)]
# rightascension Rightascension [0,2Pi]
# declination Declination [-Pi/2,Pi/2]
# polarisation Polarisation angle [0,Pi]
# Spin Parameters:
# a_spin1 Spin1 magnitude [-1,1] for aligned, [0,1] for precessing
# a_spin2 Spin2 magnitude [-1,1] for aligned, [0,1] for precessing
# tilt_spin1 Angle between spin1 and orbital angular momentum [0,Pi]
# tilt_spin2 Angle between spin2 and orbital angular momentum [0, Pi]
# phi_12 Difference between spins' azimuthal angles [0,2Pi]
# phi_jl Difference between total and orbital angular momentum azimuthal angles [0,2Pi]
# Tidal parameters (requires tidal= or tidalT=):
# lambda1 lambda1 [0,3000]
# lambda2 lambda2 [0,3000]
# lambdaT lambdaT [0,3000]
# dLambdaT dLambdaT [-500,500]
# Equation of State parameters (requires 4PolyEOS=):
# logp1 logp1 [33.6,35.4]
# gamma1 gamma1 [2.0,4.5]
# gamma2 gamma2 [1.1,4.5]
# gamma3 gamma3 [1.1,4.5]
# Settings for allowed component masses in Solar Masses, with default values
# comp-max=30.0
# comp-min=1.0
# Allowed total masses in Msun (note, used together with component masses, mc,q,eta priors may lead to inconsistencies. Be careful!)
# mtotal-max=35
# mtotal-min=2
# Setting time prior [seconds]
# dt=0.1
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment