Commit 6a27b6f7 authored by Marc Arene's avatar Marc Arene
Browse files

changed inputs from:

`minimum_frequency, interferometers` to `interferometers, 
waveform_arguments`
parent 7c5dbb56
......@@ -220,8 +220,12 @@ if __name__=='__main__':
start_time=start_time)
# WAVEFORM GENERATION
conf.waveform_arguments['maximum_frequency'] = conf.maximum_frequency_injected_waveform
template_ifos_injected = bilby_wv.WaveForm_ThreeDetectors(injection_parameters, minimum_frequency, interferometers)
waveform_arguments = {}
waveform_arguments['minimum_frequency'] = injection_parameters['meta']['minimum_frequency']
waveform_arguments['waveform_approximant'] = injection_parameters['meta']['approximant']
waveform_arguments['reference_frequency'] = injection_parameters['meta']['reference_frequency']
waveform_arguments['maximum_frequency'] = injection_parameters['meta']['maximum_frequency_injected_waveform']
template_ifos_injected = bilby_wv.WaveForm_ThreeDetectors(injection_parameters, interferometers, waveform_arguments)
# INJECT TEMPLATE INTO NOISE STRAIN DATA
# Add the template signal to each detector's strain data
......@@ -230,13 +234,13 @@ if __name__=='__main__':
interferometers[i].fd_strain = interferometers[i].strain_data.frequency_domain_strain
# COMPUTE AND PRINT SNR
conf.waveform_arguments['maximum_frequency'] = conf.maximum_frequency_search_waveform
template_ifos = bilby_wv.WaveForm_ThreeDetectors(injection_parameters, minimum_frequency, interferometers)
waveform_arguments['maximum_frequency'] = injection_parameters['meta']['maximum_frequency_search_waveform']
template_ifos = bilby_wv.WaveForm_ThreeDetectors(injection_parameters, interferometers, waveform_arguments)
logL, Snr, opt_snr_Best = logL_snr.main(template_ifos, interferometers, True)
# COMPUTE dlogL and print it to check with the plots
dlogL = CdlogL.main(injection_parameters, start_time, minimum_frequency, interferometers, True)
dlogL = CdlogL.main(injection_parameters, start_time, interferometers, waveform_arguments, True)
# DEFINE BOUNDARY CONDITIONS FOR HMC
# The boundary conditions are given both for masses m1, m2 (solar masses) and for the distance (meters)
......@@ -253,7 +257,7 @@ if __name__=='__main__':
boundary=np.array([[boundaries['m1min'],boundaries['m1max']],[boundaries['m2min'],boundaries['m2max']],[boundaries['Dmin'] * MPC, boundaries['Dmax'] * MPC]])
FisherMatrix, scale, scale_SVD = FIM.main(injection_parameters, start_time, minimum_frequency, interferometers)
FisherMatrix, scale, scale_SVD = FIM.main(injection_parameters, start_time, interferometers, waveform_arguments)
# # re-adjust scales if FIM is incorrect
......@@ -458,7 +462,7 @@ if __name__=='__main__':
# Recompute what was logL at the end of phase1. I could have done the same with dlogL but well I preferred storing its value at the end of phase1.
parameters = paru.q_pos_to_dictionary(q_pos, start_time)
template_ifos = bilby_wv.WaveForm_ThreeDetectors(parameters, minimum_frequency, interferometers)
template_ifos = bilby_wv.WaveForm_ThreeDetectors(parameters, interferometers, waveform_arguments)
logL = bilby_utils.loglikelihood(template_ifos, interferometers)
if likelihood_gradient is not None:
......@@ -477,7 +481,7 @@ if __name__=='__main__':
# physical_parameters_end_phase1 = paru.q_pos_to_dictionary(q_pos, start_time)
# parameters_end_phase1 = injection_parameters.copy()
# parameters_end_phase1.update(physical_parameters_end_phase1)
# dlogL = CdlogL.main(parameters_end_phase1, start_time, minimum_frequency, interferometers, True)
# dlogL = CdlogL.main(parameters_end_phase1, start_time, interferometers, waveform_arguments, True)
# Remove former chain file if it existed and copy that from phase 1 to which the new chain will be appended
# if os.path.exists(fileOutChain_path):
......@@ -520,7 +524,7 @@ if __name__=='__main__':
if epsilon <= 0.001: epsilon = 0.001
if epsilon >= 0.01: epsilon = 0.01
q_pos, logL, dlogL, qpos_global_fit, dlogL_global_fit, qpos_olut_fit, dlogL_olut_fit, Accept, H_p_logL, pmom_trajs, a, H, H_0 = bht.NumericalGradientLeapfrog_ThreeDetectors(q_pos, logL, dlogL, start_time, minimum_frequency, interferometers, randGen, epsilon, length_num_traj, scale, qpos_global_fit, dlogL_global_fit, qpos_olut_fit, dlogL_olut_fit, boundary, 1, 1, traj_index, H_p_logL, pmom_trajs, likelihood_gradient)
q_pos, logL, dlogL, qpos_global_fit, dlogL_global_fit, qpos_olut_fit, dlogL_olut_fit, Accept, H_p_logL, pmom_trajs, a, H, H_0 = bht.NumericalGradientLeapfrog_ThreeDetectors(q_pos, logL, dlogL, start_time, interferometers, waveform_arguments, randGen, epsilon, length_num_traj, scale, qpos_global_fit, dlogL_global_fit, qpos_olut_fit, dlogL_olut_fit, boundary, 1, 1, traj_index, H_p_logL, pmom_trajs, likelihood_gradient)
# Evaluation of acceptante rate
traj_status = " rejected "
......@@ -771,7 +775,7 @@ if __name__=='__main__':
# the hybrid trajectories or the full numericaltrajectories.
if False:
length_traj = 50 + int(randGen.uniform() * 100)
q_pos, logL, dlogL, qpos_olut_fit, dlogL_olut_fit, Accept, H_p_logL, pmom_trajs, a, H, H_0 = bht.FullHybridGradientLeapfrog_ThreeDetectors(q_pos, logL, dlogL, start_time, minimum_frequency, interferometers, randGen, epsilon, length_traj, scale, qpos_olut_fit, dlogL_olut_fit, len(qpos_olut_fit), boundary, n_param, fit_coeff_global, flag_update_qpos_olut_fit, oluts_dict, len(qpos_olut_fit), n_fit_1, n_fit_2, SquaredScale, traj_index, H_p_logL, pmom_trajs, likelihood_gradient)
q_pos, logL, dlogL, qpos_olut_fit, dlogL_olut_fit, Accept, H_p_logL, pmom_trajs, a, H, H_0 = bht.FullHybridGradientLeapfrog_ThreeDetectors(q_pos, logL, dlogL, start_time, interferometers, waveform_arguments, randGen, epsilon, length_traj, scale, qpos_olut_fit, dlogL_olut_fit, len(qpos_olut_fit), boundary, n_param, fit_coeff_global, flag_update_qpos_olut_fit, oluts_dict, len(qpos_olut_fit), n_fit_1, n_fit_2, SquaredScale, traj_index, H_p_logL, pmom_trajs, likelihood_gradient)
traj_type = '{:16}'.format('fullhybrid')
# elif False:
# if traj_index == n_traj_fit + 1:
......@@ -780,20 +784,20 @@ if __name__=='__main__':
length_traj = 50 + int(randGen.uniform() * 100)
# length_traj = 200
# QUESTION: flagNumTable seems to be not used anymore
q_pos, logL, dlogL, Accept, a, H, H_0 = bht.AnalyticalGradientLeapfrog_ThreeDetectors(q_pos, logL, dlogL, start_time, minimum_frequency, interferometers, randGen, epsilon, length_traj, scale, fit_coeff_global, oluts_dict, boundary, n_param, flagNumTable, len(qpos_olut_fit), n_fit_1, n_fit_2, SquaredScale, traj_index, likelihood_gradient, qpos_scaler, dlogL_scalers, dlogL_fit_methods)
q_pos, logL, dlogL, Accept, a, H, H_0 = bht.AnalyticalGradientLeapfrog_ThreeDetectors(q_pos, logL, dlogL, start_time, interferometers, waveform_arguments, randGen, epsilon, length_traj, scale, fit_coeff_global, oluts_dict, boundary, n_param, flagNumTable, len(qpos_olut_fit), n_fit_1, n_fit_2, SquaredScale, traj_index, likelihood_gradient, qpos_scaler, dlogL_scalers, dlogL_fit_methods)
atr += 1
traj_type = f'analytical-{opts.fit_method}'
# elif False:
elif 0.5 <= Acc or Acc < 0.65:
length_traj = 50 + int(randGen.uniform() * 50)
q_pos, logL, dlogL, qpos_olut_fit, dlogL_olut_fit, Accept, H_p_logL, pmom_trajs, a, H, H_0 = bht.HybridGradientLeapfrog_ThreeDetectors(q_pos, logL, dlogL, start_time, minimum_frequency, interferometers, randGen, epsilon, length_traj, scale, qpos_olut_fit, dlogL_olut_fit, len(qpos_olut_fit), boundary, n_param, fit_coeff_global, flag_update_qpos_olut_fit, traj_index, H_p_logL, pmom_trajs, likelihood_gradient)
q_pos, logL, dlogL, qpos_olut_fit, dlogL_olut_fit, Accept, H_p_logL, pmom_trajs, a, H, H_0 = bht.HybridGradientLeapfrog_ThreeDetectors(q_pos, logL, dlogL, start_time, interferometers, waveform_arguments, randGen, epsilon, length_traj, scale, qpos_olut_fit, dlogL_olut_fit, len(qpos_olut_fit), boundary, n_param, fit_coeff_global, flag_update_qpos_olut_fit, traj_index, H_p_logL, pmom_trajs, likelihood_gradient)
htr += 1
traj_type = 'hybrid'
if Accept == 1: num_traj_olut += 1
# elif False:
elif Acc <= 0.5:
length_traj = 50 + int(randGen.uniform() * 50)
q_pos, logL, dlogL, qpos_global_fit, dlogL_global_fit, qpos_olut_fit, dlogL_olut_fit, Accept, H_p_logL, pmom_trajs, a, H, H_0 = bht.NumericalGradientLeapfrog_ThreeDetectors(q_pos, logL, dlogL, start_time, minimum_frequency, interferometers, randGen, epsilon, length_traj, scale, qpos_global_fit, dlogL_global_fit, qpos_olut_fit, dlogL_olut_fit, boundary, flag_update_qpos_global_fit, flag_update_qpos_olut_fit, traj_index, H_p_logL, pmom_trajs, likelihood_gradient)
q_pos, logL, dlogL, qpos_global_fit, dlogL_global_fit, qpos_olut_fit, dlogL_olut_fit, Accept, H_p_logL, pmom_trajs, a, H, H_0 = bht.NumericalGradientLeapfrog_ThreeDetectors(q_pos, logL, dlogL, start_time, interferometers, waveform_arguments, randGen, epsilon, length_traj, scale, qpos_global_fit, dlogL_global_fit, qpos_olut_fit, dlogL_olut_fit, boundary, flag_update_qpos_global_fit, flag_update_qpos_olut_fit, traj_index, H_p_logL, pmom_trajs, likelihood_gradient)
ntr += 1
traj_type = 'numerical'
if Accept == 1:
......@@ -804,14 +808,14 @@ if __name__=='__main__':
# length_traj = 50 + int(randGen.uniform() * 100)
length_traj = 1000
# length_traj = 200
q_pos, logL, dlogL, PosFit_traj, DlogLFit_traj, DlogLFit_traj_ana_along_num, Accept = bht.AnaAlongNumGradientLeapfrog_ThreeDetectors(q_pos, logL, dlogL, start_time, minimum_frequency, interferometers, randGen, epsilon, length_traj, scale, fit_coeff_global, oluts_dict, boundary, n_param, len(qpos_olut_fit), n_fit_1, n_fit_2, SquaredScale, likelihood_gradient, None)
q_pos, logL, dlogL, PosFit_traj, DlogLFit_traj, DlogLFit_traj_ana_along_num, Accept = bht.AnaAlongNumGradientLeapfrog_ThreeDetectors(q_pos, logL, dlogL, start_time, interferometers, waveform_arguments, randGen, epsilon, length_traj, scale, fit_coeff_global, oluts_dict, boundary, n_param, len(qpos_olut_fit), n_fit_1, n_fit_2, SquaredScale, likelihood_gradient, None)
plots.plot_num_and_ana_gradients(DlogLFit_traj.T, DlogLFit_traj_ana_along_num.T, opts_dict=config_dict['hmc'], save=True, show_plot=False, outdir=outdir + '/plots', alongNum=True)
# q_pos, logL, dlogL, PosFit_traj, DlogLFit_traj, PosFit_traj_ana, DlogLFit_traj_ana_vs_num, Accept = bht.AnaVsNumGradientLeapfrog_ThreeDetectors(q_pos, logL, dlogL, start_time, minimum_frequency, interferometers, randGen, epsilon, length_traj, scale, fit_coeff_global, oluts_dict, boundary, n_param, len(qpos_olut_fit), n_fit_1, n_fit_2, SquaredScale)
# q_pos, logL, dlogL, PosFit_traj, DlogLFit_traj, PosFit_traj_ana, DlogLFit_traj_ana_vs_num, Accept = bht.AnaVsNumGradientLeapfrog_ThreeDetectors(q_pos, logL, dlogL, start_time, interferometers, waveform_arguments, randGen, epsilon, length_traj, scale, fit_coeff_global, oluts_dict, boundary, n_param, len(qpos_olut_fit), n_fit_1, n_fit_2, SquaredScale)
# plots.plot_num_and_ana_gradients(DlogLFit_traj.T, DlogLFit_traj_ana_vs_num.T, opts_dict=config_dict['hmc'], save=True, show_plot=False, outdir=outdir + '/plots', alongNum=False)
# plots.plot_num_vs_ana_positions(PosFit_traj.T, PosFit_traj_ana.T, opts_dict=config_dict['hmc'], save=True, show_plot=False, outdir=outdir + '/plots')
# p_mom = randGen.normal(loc=0, scale=1.0, size=9)
# q_pos_meth1, DlogLFit_traj_meth1, logL_traj_meth1, q_pos_meth2, DlogLFit_traj_meth2 = bht.NumGradientBothMethods_ThreeDetectors_Trajectory(q_pos, p_mom, dlogL, start_time, minimum_frequency, interferometers, epsilon, length_traj, scale, boundary)
# q_pos_meth1, DlogLFit_traj_meth1, logL_traj_meth1, q_pos_meth2, DlogLFit_traj_meth2 = bht.NumGradientBothMethods_ThreeDetectors_Trajectory(q_pos, p_mom, dlogL, start_time, interferometers, waveform_arguments, epsilon, length_traj, scale, boundary)
# plots.plot_num_gradients_meth12_and_logL(DlogLFit_traj_meth1.T, DlogLFit_traj_meth2.T, logL_traj_meth1, opts_dict=config_dict['hmc'], save=True, show_plot=True, outdir=outdir + '/plots')
sys.exit()
......@@ -848,7 +852,7 @@ if __name__=='__main__':
epsilon_stuck = 2.5e-3
length_traj = 20 + int(randGen.uniform() * 80)
# print(' Before: min = {:10e}'.format(np.abs(np.asarray(dlogL_olut_fit)).min()))
q_pos, logL, dlogL, qpos_olut_fit, dlogL_olut_fit, Accept, H_p_logL, pmom_trajs, a, H, H_0 = bht.HybridGradientLeapfrog_ThreeDetectors(q_pos, logL, dlogL, start_time, minimum_frequency, interferometers, randGen, epsilon_stuck, length_traj, scale, qpos_olut_fit, dlogL_olut_fit, len(qpos_olut_fit), boundary, n_param, fit_coeff_global, flag_update_qpos_olut_fit, traj_index, H_p_logL, pmom_trajs, likelihood_gradient)
q_pos, logL, dlogL, qpos_olut_fit, dlogL_olut_fit, Accept, H_p_logL, pmom_trajs, a, H, H_0 = bht.HybridGradientLeapfrog_ThreeDetectors(q_pos, logL, dlogL, start_time, interferometers, waveform_arguments, randGen, epsilon_stuck, length_traj, scale, qpos_olut_fit, dlogL_olut_fit, len(qpos_olut_fit), boundary, n_param, fit_coeff_global, flag_update_qpos_olut_fit, traj_index, H_p_logL, pmom_trajs, likelihood_gradient)
# print(' After: min = {:10e}'.format(np.abs(np.asarray(dlogL_olut_fit)).min()))
traj_type = 'hybrid'
htr += 1
......@@ -856,7 +860,7 @@ if __name__=='__main__':
else: # if hybrid trajectory did not work, switch to full numericaltrajectory
epsilon_stuck = 2.5e-3
length_traj = 20 + int(randGen.uniform() * 80)
q_pos, logL, dlogL, qpos_global_fit, dlogL_global_fit, qpos_olut_fit, dlogL_olut_fit, Accept, H_p_logL, pmom_trajs, a, H, H_0 = bht.NumericalGradientLeapfrog_ThreeDetectors(q_pos, logL, dlogL, start_time, minimum_frequency, interferometers, randGen, epsilon_stuck, length_traj, scale, qpos_global_fit, dlogL_global_fit, qpos_olut_fit, dlogL_olut_fit, boundary, flag_update_qpos_global_fit, flag_update_qpos_olut_fit, traj_index, H_p_logL, pmom_trajs, likelihood_gradient)
q_pos, logL, dlogL, qpos_global_fit, dlogL_global_fit, qpos_olut_fit, dlogL_olut_fit, Accept, H_p_logL, pmom_trajs, a, H, H_0 = bht.NumericalGradientLeapfrog_ThreeDetectors(q_pos, logL, dlogL, start_time, interferometers, waveform_arguments, randGen, epsilon_stuck, length_traj, scale, qpos_global_fit, dlogL_global_fit, qpos_olut_fit, dlogL_olut_fit, boundary, flag_update_qpos_global_fit, flag_update_qpos_olut_fit, traj_index, H_p_logL, pmom_trajs, likelihood_gradient)
ntr += 1
traj_type = 'numerical'
if Accept == 1:
......
......@@ -15,7 +15,7 @@ import Codes.set_psds as set_psds
def main(injection_parameters, start_time, minimum_frequency, interferometers):
def main(injection_parameters, start_time, interferometers, waveform_arguments):
######################################
# FISHER MATRICES #
......@@ -27,7 +27,7 @@ def main(injection_parameters, start_time, minimum_frequency, interferometers):
# other option : inverse of square root of diagonal FIM terms
# Compute the Fisher Matrix
FisherMatrix = bilby_wv.FisherMatrix_ThreeDetectors(injection_parameters, start_time, minimum_frequency, interferometers)
FisherMatrix = bilby_wv.FisherMatrix_ThreeDetectors(injection_parameters, start_time, interferometers, waveform_arguments)
# If I want to run the hmc on less than 9 dimensions, my 9x9 FIM will be really singular, hence I should remove the rows and columns of the parameters which I keep fixed to get a non-singular sub_FIM which I can invert, otherwise just run the last commented two lines.
ixgrid = np.ix_(conf.param_indices, conf.param_indices)
......@@ -141,7 +141,7 @@ if __name__=='__main__':
# COMPUTE THE FISHER MATRICES
FisherMatrix, scale, scale_SVD = main(injection_parameters, start_time, minimum_frequency, interferometers)
FisherMatrix, scale, scale_SVD = main(injection_parameters, start_time, interferometers, waveform_arguments)
print("\nValues of FIM:")
......
......@@ -16,9 +16,9 @@ import Codes.set_psds as set_psds
import Codes.logL_snr as logL_snr
def main(injection_parameters, start_time, minimum_frequency, interferometers, do_prints=False):
def main(injection_parameters, start_time, interferometers, waveform_arguments, do_prints=False):
dlogL, logL = bilby_wv.dlogL_ThreeDetectors(injection_parameters, start_time, minimum_frequency, interferometers)
dlogL, logL = bilby_wv.dlogL_ThreeDetectors(injection_parameters, start_time, interferometers, waveform_arguments)
prefix = '\ndlogL with bilby '
if do_prints:
......@@ -107,7 +107,7 @@ if __name__=='__main__':
# WAVEFORM GENERATION
injection_parameters['meta']['maximum_frequency_generated_waveform'] = injection_parameters['meta']['maximum_frequency_injected_waveform']
template_ifos_injected = bilby_wv.WaveForm_ThreeDetectors(injection_parameters, minimum_frequency, interferometers)
template_ifos_injected = bilby_wv.WaveForm_ThreeDetectors(injection_parameters, interferometers, waveform_arguments)
# INJECT TEMPLATE INTO NOISE STRAIN DATA
# Add the template signal to each detector's strain data
......@@ -118,5 +118,5 @@ if __name__=='__main__':
# COMPUTE AND PRINT SNR
# Generate a template up to frequency of analysis
injection_parameters['meta']['maximum_frequency_generated_waveform'] = injection_parameters['meta']['maximum_frequency_search_waveform']
template_ifos = bilby_wv.WaveForm_ThreeDetectors(injection_parameters, minimum_frequency, interferometers)
template_ifos = bilby_wv.WaveForm_ThreeDetectors(injection_parameters, interferometers, waveform_arguments)
logLBest, mf_snr_best, opt_snr_Best = main(template_ifos, interferometers, True)
......@@ -228,11 +228,11 @@ def ini_file_to_dict(file_path='../examples/GW170817.ini'):
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.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.maximum_frequency_injected_waveform = maximum_frequency_injected_waveform
......
This diff is collapsed.
This diff is collapsed.
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