corner plots and walkers

Codes/BNS_HMC_9D.py

@@ -17,6 +17,7 @@ import Library.Fit_NR as fnr
 import Library.config as conf
 import Library.bilby_waveform as bilby_wv
 import Library.bilby_utils as bilby_utils
+import Library.plots as plots
 
 import Codes.set_injection_parameters as set_inj
 import Codes.set_psds as set_psds
@@ -350,12 +351,16 @@ if __name__=='__main__':
                 print("{:3} - {} - Acc = {:.1%}".format(l, traj_status, Acc), end='')
                 if opts.stdout_file:
                     print("{:3} - {} - Acc = {:.1%}".format(l, traj_status, Acc), end='', file=stdout_terminal) # print also in terminal in order to know in real time how many trajectories have been computed.
-
+                print('')
         fileOutChain.closed
 
         stop = time.time()
         duration1 = stop - start
 
+        samples = plots.BNS_HMC_Chain_dat_to_walkers(file_path=fileOutChain_path)
+        search_parameter_keys = [conf.parameter_keys[i] for i in conf.param_indices]
+        plots.plot_walkers_and_corner_from_samples(samples, injection_parameters, search_parameter_keys, save=True, outdir=outdir+'plots')
+
         pt_Fit_good_filename = outdir + "pt_Fit_good.dat"
         dlogL_Fit_good_filename = outdir + "dlogL_Fit_good.dat"
         pt_Fit_bad_filename = outdir + "pt_Fit_bad.dat"

Codes/set_injection_parameters.py

@@ -96,6 +96,12 @@ def ini_file_to_dict(file_path='../examples/GW170817.ini'):
 
     injection_parameters['meta'] = meta_params
 
+
+    # 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*np.pi)
+    injection_parameters['tc_3p5PN'] = tc_3p5PN
+
+
     return injection_parameters
 
 

Library/config.py

-import bilby.gw.detector as det
 import numpy as np
 
 random_numbers = []
 
-H1 = det.InterferometerList(['H1'])[0]
-H1_vertex = H1.vertex
-
 
 # indices of interests to run the hmc on a subset of the 9 parameters:
 # 0 = cos(iota)
@@ -20,6 +16,7 @@ H1_vertex = H1.vertex
 param_names = ["iota", "phi_c", "psi", "D", "Mc", "mu", "dec", "ra", "tc"]
 fparam_names = ["cos(iota)", "phi_c", "psi", "ln(D)", "ln(Mc)", "ln(mu)", "sin(dec)", "ra", "ln(tc)"]
 param_dict_keys = ["iota", "phase", "psi", "luminosity_distance", "chirp_mass", "reduced_mass", "dec", "ra", "geocent_time"]
+parameter_keys=['theta_jn', 'phi_c', 'psi', 'luminosity_distance', 'chirp_mass', 'reduced_mass', 'dec', 'ra', 'tc_3p5PN']
 
 def identity(param):
     return param

Library/plots.py

 import matplotlib.pyplot as plt
 import numpy as np
 import ipdb
+import corner
 
-def plot_chain(chain_array, parameter_name, figsize=(15,5)):
+import sys
+sys.path.append('../')
+import Library.config as conf
 
-    chain_iterations = np.arange(len(chain_array))
+def perso_plot_walker(walk, parameter_name, figsize=(15,5)):
+
+    idxs = np.arange(len(walk))
     plt.figure(figsize=figsize)
-    plt.plot(chain_iterations, chain_array)
+    plt.plot(idxs, walk)
     plt.xlabel("Iterations")
     plt.ylabel(parameter_name)
     plt.show()
 
-def plot_chains(chains_arrays, parameters_names, figsize=(15,15), title=""):
 
-    nb_of_parameters = len(chains_arrays)
+def perso_plot_walkers(walks, parameters_names, figsize=(15,15), title=""):
+
+    nb_of_parameters = len(walks)
     nb_plots_first_column = int((nb_of_parameters + 1) / 2)
 
-    chain_iterations = np.arange(len(chains_arrays[0]))
+    idxs = np.arange(len(walks[0]))
     fig, axs = plt.subplots(nb_plots_first_column, 2, sharex=True, figsize=figsize)
 
     fig.suptitle(title)
 
     for i in range(nb_plots_first_column):
-        axs[i][0].plot(chain_iterations, chains_arrays[i], linewidth=1)
+        axs[i][0].plot(idxs, walks[i], linewidth=1)
         # axs[i].set_xlabel("Iterations")
         axs[i][0].set_ylabel(parameters_names[i])
 
     for i in range(nb_plots_first_column, nb_of_parameters):
-        axs[i - nb_plots_first_column][1].plot(chain_iterations, chains_arrays[i], linewidth=1)
+        axs[i - nb_plots_first_column][1].plot(idxs, walks[i], linewidth=1)
         # axs[i].set_xlabel("Iterations")
         axs[i - nb_plots_first_column][1].set_ylabel(parameters_names[i])
 
@@ -41,11 +47,67 @@ def plot_chains(chains_arrays, parameters_names, figsize=(15,15), title=""):
 
     plt.show()
 
+
+def perso_plot_corner(samples, labels=None):
+
+    fig = corner.corner(samples, labels=labels)
+    plt.show()
+
+
+def BNS_HMC_Chain_dat_to_walkers(file_path):
+    HMC_Chain = np.loadtxt(file_path)
+    search_param_indices = np.asarray(conf.param_indices)
+    return HMC_Chain[:, search_param_indices+1]
+
+def init_result_object(samples, injection_parameters, search_parameter_keys=['theta_jn', 'phase', 'psi', 'luminosity_distance', 'chirp_mass', 'reduced_mass', 'dec', 'ra', 'geocent_time'], label='no_label', outdir=''):
+    import pandas as pd
+    import bilby
+    result = bilby.core.result.Result()
+
+    result.label = label
+    result.samples = samples
+    result.injection_parameters = injection_parameters
+    result.search_parameter_keys = search_parameter_keys
+    default_latex_labels = bilby.gw.prior.Prior._default_latex_labels
+    default_latex_labels['reduced_mass'] = '$\\mu$'
+    default_latex_labels['phi_c'] = '$\\phi_c$'
+    default_latex_labels['tc_3p5PN'] = '$\\t_c$'
+    result.parameter_labels = [default_latex_labels[k] for k in result.search_parameter_keys]
+    result.parameter_labels_with_unit = result.parameter_labels
+    data_frame = pd.DataFrame(result.samples, columns=result.search_parameter_keys)
+    result.posterior = data_frame
+
+    # Needed in order to use result.plot_walkers()
+    result.walkers = np.array([samples])
+    result.nburn = 0
+
+    return result
+
+def plot_walkers_and_corner_from_samples(samples, injection_parameters, search_parameter_keys, save=True, outdir='./plots'):
+
+    result = init_result_object(samples, injection_parameters, search_parameter_keys=search_parameter_keys)
+    if save:
+        result.plot_walkers(outdir=outdir)
+    result.plot_corner(save=save, outdir=outdir)
+
+    # result.plot_single_density('ra', truth=injection_parameters['ra'], save=False)
+    # result.plot_marginals()
+
 if __name__=='__main__':
-    HMC_Chain = np.loadtxt('/Users/marcarene/Documents/APC/Code/myPython/YannsCtoPython/newPythonCode/.output_data/GW170817/bilbyoriented/SUB-D_012345678/PSD_1/0e+00_500_2000_200_200_0.005_dlogL1/BNS_HMC_Chain.dat').T
+    import config as conf
+    import sys
+    sys.path.append('../')
+    import Codes.set_injection_parameters as set_inj
 
-    fparam_names = ["cos(iota)", "phi_c", "psi", "ln(D)", "ln(Mc)", "ln(mu)", "sin(dec)", "ra", "ln(tc)"]
+    outdir_500_traj = '/Users/marcarene/Documents/APC/Code/myPython/YannsCtoPython/newPythonCode/.output_data/GW170817/bilbyoriented/SUB-D_012345678/PSD_1/0e+00_500_2000_200_200_0.005_dlogL1'
+    file_path_500_traj = outdir_500_traj + '/BNS_HMC_Chain.dat'
 
-    chains_arrays = HMC_Chain[1:10]
 
-    plot_chains(chains_arrays, fparam_names)
+    injection_parameters = set_inj.ini_file_to_dict()
+    search_parameter_keys=['theta_jn', 'phi_c', 'psi', 'luminosity_distance', 'chirp_mass', 'reduced_mass', 'dec', 'ra', 'tc_3p5PN']
+    # search_parameter_keys=['luminosity_distance']
+    samples = BNS_HMC_Chain_dat_to_walkers(file_path=file_path_500_traj)
+
+    # import IPython; IPython.embed()
+    plot_walkers_and_corner_from_samples(samples, injection_parameters, search_parameter_keys, save=False, outdir=outdir_500_traj+'/plots')
+    plt.show()

Library/utils.py

-def get_analytical_timing_for_dlogL(t_h_source=9.1, t_get_detector_response_geocent=1.1, t_get_fd_time_translation=3.5, t_noise_weighted_inner_product=3.7, t_dh=0.7, nb_of_ifos=3, n_traj_phase1=1000, lengthTnum=200):
+def get_analytical_timing_for_dlogL(t_h_source=1, t_get_detector_response_geocent=0.1, t_get_fd_time_translation=0.1, t_noise_weighted_inner_product=0.1, t_dh=0.07, nb_of_ifos=3, n_traj_phase1=1000, lengthTnum=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():
 
@@ -43,7 +43,7 @@ def get_analytical_timing_for_dlogL(t_h_source=9.1, t_get_detector_response_geoc
     print("dlogL = {:.0f}ms".format(dlogL))
     print("1 trajectory of {} steps = {:.1f}s".format(lengthTnum, one_traj))
     print("Phase 1 with {} trajectories = {:.1f}hrs".format(n_traj_phase1, phase_1/3600))
-    import IPython; IPython.embed()
+    # import IPython; IPython.embed()
 
 
 if __name__=='__main__':