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

using get_inj_parameters_and_analysis_dict()

parent e05968e8
......@@ -135,19 +135,16 @@ if __name__=='__main__':
# SET THE INJECTION PARAMETERS AND THE ANALYSIS PARAMETERS
injection_parameters = set_inj.ini_file_to_dict(opts.inj_file)
injection_parameters, ext_analysis_dict = set_inj.get_inj_parameters_and_analysis_dict(opts.inj_file, **config_dict['analysis'])
pu.save_dict_to_json(injection_parameters, outdir + 'injection_parameters.json')
print("\nInjected parameters derived from {} are:".format(opts.inj_file))
pu.print_dict(injection_parameters)
ext_analysis_dict = set_inj.compute_extended_analysis_dict(injection_parameters['mass_1'], injection_parameters['mass_2'], injection_parameters['geocent_time'], injection_parameters['chirp_mass'], **config_dict['analysis'])
pu.save_dict_to_json(ext_analysis_dict, outdir + 'ext_analysis_dict.json')
print("\nAnalysis parameters derived from the component masses are:")
pu.print_dict(ext_analysis_dict)
start_time = ext_analysis_dict['start_time']
# import IPython; IPython.embed(); sys.exit()
# INITIALIZE THE THREE INTERFEROMETERS
interferometers = bilby.gw.detector.InterferometerList(ifo_chosen)
......@@ -493,7 +490,7 @@ if __name__=='__main__':
fileOutChain.write("{} \t {:.6f} \t {:.6f} \t {:.6f} \t {:.6f} \t {:14.12e} \t {:14.12e} \t {:.6f} \t {:.6f} \t {:14.12e} \t {:.6f} \t {:.6f} \t{:.6f} \t {:.6f} \t {:e}\n".format(traj_index, q_pos[0] , q_pos[1], q_pos[2], q_pos[3]/MPC, q_pos[4], q_pos[5], q_pos[6], q_pos[7], q_pos[8], logL, 100.0*Acc, q_pos[9], q_pos[10], epsilon))
if n_traj_fit > 20 and traj_index%int(n_traj_fit/4)==0:
if n_traj_fit >= 20 and traj_index%int(n_traj_fit/4)==0:
# Make also a corner plot of all the points visited during phase1's accepted trajectories which hence will be used during phase2 for fitting
injection_fparam = injection_parameters.copy()
injection_fparam['cos(theta_jn)'] = np.cos(injection_parameters['theta_jn'])
......@@ -521,27 +518,6 @@ if __name__=='__main__':
stop = time.time()
duration1 += stop - start
# Make a walker and corner plots of the chain in phase1
samples = plots.BNS_HMC_Chain_dat_to_walkers(file_path=fileOutChain_path)
# Condition meant for short runs of less 9 trajectories where the corner plot would bug because there are less samples than dimensions...
if samples.shape[0] > samples.shape[1]:
plots.plot_walkers_and_corner_from_samples(samples, injection_parameters, search_parameter_keys, save=True, label='{}_samples'.format(len(samples)), outdir=outdir_phase1+'plots')
# Make also a corner plot of all the points visited during phase1's accepted trajectories which hence will be used during phase2 for fitting
injection_fparam = injection_parameters.copy()
injection_fparam['cos(theta_jn)'] = np.cos(injection_parameters['theta_jn'])
injection_fparam['ln(D)'] = np.log(injection_parameters['luminosity_distance'])
injection_fparam['ln(Mc)'] = np.log(injection_parameters['chirp_mass'])
injection_fparam['ln(mu)'] = np.log(injection_parameters['reduced_mass'])
injection_fparam['sin(dec)'] = np.sin(injection_parameters['dec'])
injection_fparam['ln(tc)'] = np.log(injection_parameters['geocent_time']-start_time)
plots.plot_walkers_and_corner_from_samples(np.asarray(qpos_global_fit)[:, np.asarray(conf.param_indices)], injection_fparam, search_fparam_keys, save=True, label='{}_samples'.format(len(qpos_global_fit)), outdir=outdir_phase1+'plots')
# from scipy import stats
# z = np.abs(stats.zscore(qpos_global_fit))
# threshold = 3
# qpos_global_fit_no_outliers = qpos_global_fit[(z < threshold).all(axis=1)]
# plots.plot_walkers_and_corner_from_samples(qpos_global_fit_no_outliers[:, np.asarray(conf.param_indices)], injection_fparam, search_fparam_keys, save=True, label='{}_samples'.format(len(qpos_global_fit_no_outliers)), outdir=outdir_phase1+'plots')
# Saving data at the end of phase1
phase1_status = dict(count=count, Acc=Acc, duration=duration1, n_phase1_already_run=l_end - 1)
pu.save_dict_to_json(phase1_status, outdir_phase1 + 'status.json')
......@@ -557,7 +533,6 @@ if __name__=='__main__':
if os.path.exists(fileOutChainPhase1_path):
subprocess.run('rm ' + fileOutChainPhase1_path, shell=True)
subprocess.run('cp ' + fileOutChain_path + ' ' + fileOutChainPhase1_path, shell=True)
if conf.debug:
H_p_logL_filename = outdir_phase1 + "H_p_logL.dat"
pmom_trajs_filename = outdir_phase1 + "pmom_trajs.dat"
......@@ -565,6 +540,31 @@ if __name__=='__main__':
np.savetxt(pmom_trajs_filename, pmom_trajs)
# Make a walker and corner plots of the chain in phase1
samples = plots.BNS_HMC_Chain_dat_to_walkers(file_path=fileOutChain_path)
# Condition meant for short runs of less 9 trajectories where the corner plot would bug because there are less samples than dimensions...
if samples.shape[0] > samples.shape[1]:
# 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)
injection_parameters['tc_3p5PN'] = ext_analysis_dict['tc_3p5PN']
plots.plot_walkers_and_corner_from_samples(samples, injection_parameters, search_parameter_keys, save=True, label='{}_samples'.format(len(samples)), outdir=outdir_phase1+'plots')
# Make also a corner plot of all the points visited during phase1's accepted trajectories which hence will be used during phase2 for fitting
injection_fparam = injection_parameters.copy()
injection_fparam['cos(theta_jn)'] = np.cos(injection_parameters['theta_jn'])
injection_fparam['ln(D)'] = np.log(injection_parameters['luminosity_distance'])
injection_fparam['ln(Mc)'] = np.log(injection_parameters['chirp_mass'])
injection_fparam['ln(mu)'] = np.log(injection_parameters['reduced_mass'])
injection_fparam['sin(dec)'] = np.sin(injection_parameters['dec'])
injection_fparam['ln(tc)'] = np.log(injection_parameters['geocent_time']-start_time)
plots.plot_walkers_and_corner_from_samples(np.asarray(qpos_global_fit)[:, np.asarray(conf.param_indices)], injection_fparam, search_fparam_keys, save=True, label='{}_samples'.format(len(qpos_global_fit)), outdir=outdir_phase1+'plots')
# from scipy import stats
# z = np.abs(stats.zscore(qpos_global_fit))
# threshold = 3
# qpos_global_fit_no_outliers = qpos_global_fit[(z < threshold).all(axis=1)]
# plots.plot_walkers_and_corner_from_samples(qpos_global_fit_no_outliers[:, np.asarray(conf.param_indices)], injection_fparam, search_fparam_keys, save=True, label='{}_samples'.format(len(qpos_global_fit_no_outliers)), outdir=outdir_phase1+'plots')
print("Total time for numerical derivative part of the algorithm : {:.3} hours".format(duration1/3600))
print(f'Average time to compute 1 numerical trajectory: {duration1 / n_traj_fit:.2} sec')
......
......@@ -90,18 +90,17 @@ def compute_time_from_flow_to_fhigh_3p5PN(f_low, f_high, mass_1, mass_2):
return tc_flow - tc_fhigh
def ini_file_to_dict(file_path='../examples/GW170817.ini'):
def ini_file_to_dict(inj_file_path='../examples/GW170817.ini'):
config = ConfigParser()
config.read(file_path)
config.read(inj_file_path)
injection_parameters = pu.config_parser_to_dict(config)['parameters']
injection_parameters = conv.generate_mass_parameters(injection_parameters)
injection_parameters['reduced_mass'] = paru.component_masses_to_reduced_mass(injection_parameters['mass_1'], injection_parameters['mass_2'])
# 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)
# injection_parameters['phi_c'] = (injection_parameters['phase'] * 2) % (2*np.pi)
# injection_parameters['tc_3p5PN'] = tc_3p5PN
return injection_parameters
def compute_extended_analysis_dict(mass_1, mass_2, geocent_time, chirp_mass, **analysis_kwargs):
......@@ -232,6 +231,18 @@ def compute_extended_analysis_dict(mass_1, mass_2, geocent_time, chirp_mass, **a
return ext_analysis_dict
def get_inj_parameters_and_analysis_dict(inj_file_path='../examples/GW170817.ini', **analysis_kwargs):
injection_parameters = ini_file_to_dict(inj_file_path)
ext_analysis_dict = compute_extended_analysis_dict(injection_parameters['mass_1'], injection_parameters['mass_2'], injection_parameters['geocent_time'], injection_parameters['chirp_mass'], **analysis_kwargs)
# 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)
injection_parameters['tc_3p5PN'] = ext_analysis_dict['tc_3p5PN']
return injection_parameters, ext_analysis_dict
if __name__=='__main__':
from optparse import OptionParser
......
......@@ -269,9 +269,12 @@ d2_other_param_2 = (Matrix_ClosestPoint[other_param_index_2] - q_pos_fit[other_p
- `[DONE]` compute covariance matrix of the scattered points from phase1 and compare to the scales we use
## TO DO RECAP
- `[TO DO !]` Run on Binary1.in and check if I am getting a bimodality in phic-psi. Am I getting the same distributions as Yann's?
- `[TO DO !]` run phase1 with IMRPhenomD. IMRPhenomP has apparently a hack on to define phic because of the way they handle precession.
- `[TO DO !]` prepare the steps to install gwhmc and bilby on Ed's laptop
- `[TO DO !]` Run on Binary1.in and check if I am getting a bimodality in phic-psi. Am I getting the same distributions as Yann's? Do it with IMRPhenomD
- `[TO DO !]` Compare forward difference wrt central difference in the gradients computations:
- Firstly to know which offset I should use, I should compare on ~100 points (which ones?) which offset gives the best results
- Then generate randomly ~10000 points in the 9 dimensional parameter space. For each of these compute the 9 gradients both with forward and central-difference and compute the relative error and absolute error between the two. Or should it be the absolute error divided by the scale of the gradient for that parameter ?
- `[DONE]` run phase1 with IMRPhenomD. IMRPhenomP has apparently a hack on to define phic because of the way they handle precession.
- `[DONE]` prepare the steps to install gwhmc and bilby on Ed's laptop
- `[TO DO !]` marginalize over `phic`
- `[TO DO !]` start adding the spins !
- `[TO DO !]` In BNS_HMC_9D.py, split the initialization into callable chunks
......
[analysis]
# Interferometers to run the analysis on.
ifos = H1,L1,V1
approx = IMRPhenomD
approximant = IMRPhenomD
minimum_frequency = 30.0
maximum_frequency = 2048.0
reference_frequency = 20.0
......
......@@ -22,9 +22,9 @@ phase_marginalization = False
[hmc]
# Total number of trajectories of HMC
n_traj_hmc_tot = 10
n_traj_hmc_tot = 125
# Number of trajectories using numericalgradients in phase 1.
n_traj_fit = 10
n_traj_fit = 125
# Number of iterations if phase1 loop to go over for this run. This permits to split phase1 into several chunks. If set to 0, then the value will be set to that of n_traj_fit
n_phase1_for_this_run = 800
# Number of points used to sub-select points in Ordered Look-Up Tables for local fit bimodal parameters like (cos(inc), psi, logD)
......
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