Commit 4ac4e73f authored by Marc Arene's avatar Marc Arene
Browse files

Phase1 Acc with ROQ = 96%

parent 2a741153
......@@ -13,12 +13,14 @@ import bilby
from Headers.PN_Coefficients import * # Headers.Constants already imported in PN_Coeff
import Library.python_utils as pu
import Library.param_utils as paru
import Library.roq_utils as roqu
import Library.BNS_HMC_Tools as bht
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 Library.likelihood_gradient as lg
import Codes.set_injection_parameters as set_inj
import Codes.set_psds as set_psds
......@@ -252,23 +254,123 @@ if __name__=='__main__':
FisherMatrix, scale, scale_SVD = FIM.main(injection_parameters, start_time, minimum_frequency, interferometers)
# # re-adjust scales if FIM is incorrect
# if scale[0] >1.0: scale[0] = 1.0 # assume 50% error in cosinc
# if scale[1] > PI: scale[1] = PI # assume 50% error in phi_c
# if scale[2]>PIb2: scale[2] = PIb2 # assume 50% error in psi
# if scale[3]> 0.5: scale[3] = 0.5 # assume 50% error in lnDL
# re-adjust scales if FIM is incorrect
if scale[0] >1.0: scale[0] = 1.0 # assume 50% error in cosinc
if scale[1] > PI: scale[1] = PI # assume 50% error in phi_c
if scale[2]>PIb2: scale[2] = PIb2 # assume 50% error in psi
if scale[3]> 0.5: scale[3] = 0.5 # assume 50% error in lnDL
if scale[0] >1.0: scale[0] = 1.0/2 # assume 50% error in cosinc
if scale[1] > PI: scale[1] = PI/2 # assume 50% error in phi_c
if scale[2]>PIb2: scale[2] = PIb2/2 # assume 50% error in psi
if scale[3]> 0.5: scale[3] = 0.5/2 # assume 50% error in lnDL
scale_filename = outdir + "scale.dat"
np.savetxt(scale_filename, scale)
# ROQ Basis used?
if 'roq' in injection_parameters.keys():
if injection_parameters['meta']['approximant'] != 'IMRPhenomPv2':
raise ValueError("Need to set approximant to IMRPhenomPv2 to work with the ROQ basis, approximant here is {}".format(injection_parameters['meta']['approximant']))
scale_factor = injection_parameters['roq']['scale_factor']
roq_directory = injection_parameters['roq']['directory']
# Load in the pieces for the linear part of the ROQ. Note you will need to
freq_nodes_linear = np.load(roq_directory + "fnodes_linear.npy") * scale_factor
# Load in the pieces for the quadratic part of the ROQ
freq_nodes_quadratic = np.load(roq_directory + "fnodes_quadratic.npy") * scale_factor
# Load the parameters describing the valid parameters for the basis.
params = injection_parameters['roq']['rescaled_params'].copy()
roq_outdir = '../.ROQ_weights/' + inj_name + '/'
weights_file_path = roq_outdir + '{:.0f}Hz_{:.0f}Hz_{:.0f}Hz_{:.0f}s_weights.json'.format(params['flow'], injection_parameters['meta']['maximum_frequency_ifo'], injection_parameters['meta']['maximum_frequency_injected_waveform'], params['seglen'])
# make ROQ waveform generator
search_waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
duration=duration, sampling_frequency=sampling_frequency,
frequency_domain_source_model=bilby.gw.source.binary_neutron_star_roq,
waveform_arguments=dict(
frequency_nodes_linear=freq_nodes_linear,
frequency_nodes_quadratic=freq_nodes_quadratic,
reference_frequency=injection_parameters['meta']['reference_frequency'], waveform_approximant=injection_parameters['meta']['approximant']),
parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_neutron_star_parameters)
# Here we add constraints on chirp mass and mass ratio to the prior, these are
# determined by the domain of validity of the ROQ basis.
priors = bilby.gw.prior.BNSPriorDict()
# for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'theta_jn', 'phase', 'psi', 'ra', 'dec', 'phi_12', 'phi_jl', 'luminosity_distance']:
for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'theta_jn', 'phase', 'psi', 'ra', 'dec', 'luminosity_distance']:
priors[key] = injection_parameters[key]
for key in ['mass_1', 'mass_2']:
priors[key].minimum = max(priors[key].minimum, params['compmin'])
priors['chirp_mass'] = bilby.core.prior.Constraint(name='chirp_mass', minimum=float(params['chirpmassmin']), maximum=float(params['chirpmassmax']))
priors['mass_ratio'] = bilby.core.prior.Constraint(0.125, 1, name='mass_ratio')
priors['geocent_time'] = bilby.core.prior.Uniform(injection_parameters['geocent_time'] - 0.1, injection_parameters['geocent_time'] + 0.1, latex_label='$t_c$', unit='s')
pu.mkdirs(outdir)
if os.access(weights_file_path, os.W_OK):
# load the weights from the file
likelihood = bilby.gw.likelihood.ROQGravitationalWaveTransient(
interferometers=interferometers, waveform_generator=search_waveform_generator,
weights=weights_file_path, priors=priors)
else:
# Load in the pieces for the linear part of the ROQ. Note you will need to
basis_matrix_linear = np.load(roq_directory + "B_linear.npy").T
# Load in the pieces for the quadratic part of the ROQ
basis_matrix_quadratic = np.load(roq_directory + "B_quadratic.npy").T
# Redoing here a bit what's already being done in set_inj because `bilby.gw.likelihood.ROQGravitationalWaveTransient()` does not accept dictionnary for roq_params...
roq_params_ndarray = np.genfromtxt(injection_parameters['roq']['params_path'], names=True)
roq_params_ndarray = roqu.rescale_params(roq_params_ndarray, scale_factor)
# roq_params_ndarray = None
likelihood = bilby.gw.likelihood.ROQGravitationalWaveTransient(
interferometers=interferometers, waveform_generator=search_waveform_generator,
linear_matrix=basis_matrix_linear, quadratic_matrix=basis_matrix_quadratic,
priors=priors, roq_params=roq_params_ndarray)
# remove the basis matrices as these are big for longer bases
del basis_matrix_linear, basis_matrix_quadratic
# write the weights to file so they can be loaded multiple times
likelihood.save_weights(weights_file_path)
search_parameters_keys = [conf.fparam_names[i] for i in conf.param_indices]
likelihood.parameters = injection_parameters
logL = likelihood.log_likelihood_ratio()
print("")
print("likelihood.log_likelihood_ratio() = {}".format(logL))
likelihood_gradient = lg.GWTransientLikelihoodGradient(likelihood, search_parameters_keys)
dlogL_dict = likelihood_gradient.calculate_log_likelihood_gradient()
pu.print_dict({"dlogL_ROQ": dlogL_dict})
dlogL = likelihood_gradient.convert_dlogL_dict_to_np_array(dlogL_dict)
else:
likelihood_gradient = None
# Print starting quantities of the HMC
print("\n\n********** MAIN code *********\n")
q_pos = paru.injection_parameters_to_q_pos(injection_parameters, start_time)
print("{:11} | {:20} | {:12} | {:20} | {:20} | {:20} | {:14}".format("param name", "param value (q_pos)", "fparam_name", "fparam value", "scale", "dlogL", "offset python"))
print("{:11} | {:20} | {:14} | {:20} | {:20} | {:20} | {:14}".format("param name", "param value (q_pos)", "fparam_name", "fparam value", "scale", "dlogL", "offset python"))
for i in range(9):
print("{:11} | {:20.6e} | {:12} | {:20.6e} | {:20.6e} | {:20.6e} | {:14.0e}\t".format(conf.param_names[i], q_pos[i], conf.fparam_names[i], conf.param_functions[i](q_pos[i]), scale[i], dlogL[i], conf.param_offsets[i]))
print("{:11} | {:20.6e} | {:14} | {:20.6e} | {:20.6e} | {:20.6e} | {:14.0e}\t".format(conf.param_names[i], q_pos[i], conf.fparam_names[i], conf.param_functions[i](q_pos[i]), scale[i], dlogL[i], conf.param_offsets[i]))
print("\n\n")
# import IPython; IPython.embed();
......@@ -358,6 +460,10 @@ if __name__=='__main__':
template_ifos = bilby_wv.WaveForm_ThreeDetectors(parameters, minimum_frequency, interferometers)
logL = bilby_utils.loglikelihood(template_ifos, interferometers)
if likelihood_gradient is not None:
likelihood_gradient.likelihood.parameters = parameters.copy()
if os.path.exists(outdir_phase1 + 'randGen_state.json'):
print("Setting the random number generator as it was at the end of phase1...")
randGen_state_end_phase1_dict = pu.json_to_dict(outdir_phase1 + 'randGen_state.json')
......@@ -414,7 +520,7 @@ if __name__=='__main__':
if epsilon <= 0.001: epsilon = 0.001
if epsilon >= 0.01: epsilon = 0.01
q_pos, logL, dlogL, pt_fit_good, dlogL_fit_good, pt_fit_bad, dlogL_fit_bad, Accept, pt_debug, dlogL_debug, H_p_logL, pmom_trajs, a, H, H_0 = bht.NumericalGradientLeapfrog_ThreeDetectors(q_pos, logL, dlogL, start_time, minimum_frequency, interferometers, randGen, epsilon, lengthTnum, scale, pt_fit_good, dlogL_fit_good, pt_fit_bad, dlogL_fit_bad, boundary, 1, 1, traj_index, pt_debug, dlogL_debug, H_p_logL, pmom_trajs)
q_pos, logL, dlogL, pt_fit_good, dlogL_fit_good, pt_fit_bad, dlogL_fit_bad, Accept, pt_debug, dlogL_debug, H_p_logL, pmom_trajs, a, H, H_0 = bht.NumericalGradientLeapfrog_ThreeDetectors(q_pos, logL, dlogL, start_time, minimum_frequency, interferometers, randGen, epsilon, lengthTnum, scale, pt_fit_good, dlogL_fit_good, pt_fit_bad, dlogL_fit_bad, boundary, 1, 1, traj_index, pt_debug, dlogL_debug, H_p_logL, pmom_trajs, likelihood_gradient)
# Evaluation of acceptante rate
traj_status = " rejected "
......@@ -647,7 +753,7 @@ if __name__=='__main__':
else:
# length_traj = 50 + int(randGen.uniform() * 100)
length_traj = 1000
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_good, PointFit_SortCosInc, PointFit_SortPsi, PointFit_SortLogD, boundary, n_param, len(pt_fit_bad), n_fit_1, n_fit_2, SquaredScale)
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_good, PointFit_SortCosInc, PointFit_SortPsi, PointFit_SortLogD, boundary, n_param, len(pt_fit_bad), n_fit_1, n_fit_2, SquaredScale, likelihood_gradient)
plots.plot_num_and_ana_gradients(DlogLFit_traj.T, DlogLFit_traj_ana_along_num.T, opts_dict=opts.__dict__, 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_good, PointFit_SortCosInc, PointFit_SortPsi, PointFit_SortLogD, boundary, n_param, len(pt_fit_bad), n_fit_1, n_fit_2, SquaredScale)
......
......@@ -115,16 +115,16 @@ def ini_file_to_dict(file_path='../examples/GW170817.ini'):
R_isco = 6. # Orbital separation at ISCO, in geometric units. 6M for PN ISCO; 2.8M for EOB
m1_min = 1
m2_min = 1
approximant = 'TaylorF2'
# approximant = 'TaylorF2'
# approximant = 'SpinTaylorF2'
# approximant = 'IMRPhenomD'
# approximant = 'IMRPhenomPv2'
# approximant = 'IMRPhenomPv2_ROQ'
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)
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':
f_high = f_ISCO
f_high_run = f_ISCO_run
......@@ -150,10 +150,11 @@ def ini_file_to_dict(file_path='../examples/GW170817.ini'):
rescaled_params = roqu.rescale_params(params_best_basis, scale_factor)
minimum_frequency = rescaled_params['flow']
# maximum_frequency_injected_waveform = rescaled_params['fhigh']
maximum_frequency_injected_waveform = min(2048, rescaled_params['fhigh'])
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_search_waveform = maximum_frequency_ifo
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'])
......
......@@ -224,10 +224,11 @@ def Table_GradientFit_LogD(q_pos_fit, PointFit_SortLogD, n_pt_fit, n_param, n_1_
"""
return Table_GradientFit(q_pos_fit, PointFit_SortLogD, n_pt_fit, n_param, n_1_fit, n_2_fit, SquaredScale, 3)
def NumericalGradientLeapfrog_ThreeDetectors(q_pos_0, logL_0, dlogL_0, start_time, minimum_frequency, interferometers, randGen, epsilon, lengthT, scale, pt_fit_good, dlogL_fit_good, pt_fit_bad, dlogL_fit_bad, boundary, FlagUpdatePtFit_good, FlagUpdatePtFit_bad, traj_number_debug, pt_debug, dlogL_debug, H_p_logL, pmom_trajs):
def NumericalGradientLeapfrog_ThreeDetectors(q_pos_0, logL_0, dlogL_0, start_time, minimum_frequency, interferometers, randGen, epsilon, lengthT, scale, pt_fit_good, dlogL_fit_good, pt_fit_bad, dlogL_fit_bad, boundary, FlagUpdatePtFit_good, FlagUpdatePtFit_bad, traj_number_debug, pt_debug, dlogL_debug, H_p_logL, pmom_trajs, likelihood_gradient):
"""
"""
# breakpoint()
# Initialization of the momentum: draw from a gaussian distribution with mean 0 and sigma=1
p_mom_0 = randGen.normal(loc=0, scale=1.0, size=9)
......@@ -235,17 +236,21 @@ def NumericalGradientLeapfrog_ThreeDetectors(q_pos_0, logL_0, dlogL_0, start_tim
H_0 = computeHamiltonian(p_mom_0, logL_0)
# Main HMC routine with numerical gradients
q_pos, p_mom, dlogL_pos, pt_fit_traj, dlogL_fit_traj, BreakIndic, H_p_logL_traj, pmom_traj = NumericalGradientLeapfrog_ThreeDetectors_Trajectory(q_pos_0, p_mom_0, dlogL_0, logL_0, H_0, start_time, minimum_frequency, interferometers, epsilon, lengthT, scale, boundary, traj_number_debug)
q_pos, p_mom, dlogL_pos, pt_fit_traj, dlogL_fit_traj, BreakIndic, H_p_logL_traj, pmom_traj = NumericalGradientLeapfrog_ThreeDetectors_Trajectory(q_pos_0, p_mom_0, dlogL_0, logL_0, H_0, start_time, minimum_frequency, interferometers, epsilon, lengthT, scale, boundary, traj_number_debug, likelihood_gradient)
# Compute the Log-likelihood at the end if there was no NaN
if BreakIndic==0:
# breakpoint()
# Compute waveforms
# QUESTION: should I pass as an input the injection_parameters dictionary and update it here with the following parameters dictionary or is it sufficient to just give as input to bilby/lal this parameters dictionary?
parameters = paru.q_pos_to_dictionary(q_pos, start_time)
template_ifos = bilby_wv.WaveForm_ThreeDetectors(parameters, minimum_frequency, interferometers)
if likelihood_gradient is None:
parameters = paru.q_pos_to_dictionary(q_pos, start_time)
template_ifos = bilby_wv.WaveForm_ThreeDetectors(parameters, minimum_frequency, interferometers)
# Calculate Hamiltonian at final position, and evaluate Hastings ratio with appropriate priors
logL_final = bilby_utils.loglikelihood(template_ifos, interferometers)
else:
logL_final = likelihood_gradient.likelihood.log_likelihood_ratio()
# Calculate Hamiltonian at final position, and evaluate Hastings ratio with appropriate priors
logL_final = bilby_utils.loglikelihood(template_ifos, interferometers)
H = computeHamiltonian(p_mom, logL_final)
a = randGen.uniform(low=0, high=1) # draw random number
......@@ -281,7 +286,7 @@ def NumericalGradientLeapfrog_ThreeDetectors(q_pos_0, logL_0, dlogL_0, start_tim
return q_pos_0, logL_0, dlogL_0, pt_fit_good, dlogL_fit_good, pt_fit_bad, dlogL_fit_bad, Accept, pt_debug, dlogL_debug, H_p_logL, pmom_trajs, a, H, H_0
def NumericalGradientLeapfrog_ThreeDetectors_Trajectory(q_pos_0, p_mom_0, dlogL_0, logL_0, H_0, start_time, minimum_frequency, interferometers, epsilon, lengthT, scale, boundary, traj_number_debug):
def NumericalGradientLeapfrog_ThreeDetectors_Trajectory(q_pos_0, p_mom_0, dlogL_0, logL_0, H_0, start_time, minimum_frequency, interferometers, epsilon, lengthT, scale, boundary, traj_number_debug, likelihood_gradient):
"""
HMC routine for a single trajectory where the gradients are all computed numerically.
......@@ -320,6 +325,9 @@ def NumericalGradientLeapfrog_ThreeDetectors_Trajectory(q_pos_0, p_mom_0, dlogL_
q_pos_new[7] = q_pos[7] + epsilon * p_mom_new[7] * scale[7] # long
q_pos_new[8] = np.exp( np.log(q_pos[8]) + epsilon * p_mom_new[8] * scale[8]) # ln(tc)
if q_pos_new[4] < 0.8356464705882353 or 1.5306876470588235 < q_pos_new[4]:
print("Chirp mass going outside of ROQ boundaries...!")
q_pos, p_mom, BreakIndic = BoundaryCheck_On_Position(q_pos_new, q_pos, p_mom_new, p_mom, dlogL_pos, boundary, epsilon, scale, cos_inc, sin_theta)
if BreakIndic==1:
......@@ -327,8 +335,13 @@ def NumericalGradientLeapfrog_ThreeDetectors_Trajectory(q_pos_0, p_mom_0, dlogL_
print("TRAJECTORIES WENT TO NAN\n")
break
parameters = paru.q_pos_to_dictionary(q_pos, start_time)
dlogL_pos, logL = bilby_wv.dlogL_ThreeDetectors(parameters, start_time, minimum_frequency, interferometers)
if likelihood_gradient is None:
parameters = paru.q_pos_to_dictionary(q_pos, start_time)
dlogL_pos, logL = bilby_wv.dlogL_ThreeDetectors(parameters, start_time, minimum_frequency, interferometers)
else:
likelihood_gradient.likelihood.parameters = paru.q_pos_to_dictionary(q_pos, start_time)
dlogL_pos = likelihood_gradient.calculate_dlogL_np()
logL = likelihood_gradient.likelihood.log_likelihood_ratio()
# 1/2 step in momentum (leapfrog algorithm)
p_mom += 0.5 * epsilon * dlogL_pos * scale
......@@ -896,7 +909,7 @@ def AnaVsNumGradientLeapfrog_ThreeDetectors_Trajectory(q_pos_0, p_mom_0, dlogL_0
return q_pos, p_mom, dlogL_pos, pt_fit_traj, dlogL_fit_traj, q_pos_ana, p_mom_ana, pt_fit_traj_ana, dlogL_fit_traj_ana, BreakIndic
def AnaAlongNumGradientLeapfrog_ThreeDetectors(q_pos_0, logL_0, dlogL_0, start_time, minimum_frequency, interferometers, randGen, epsilon, lengthT, scale, fit_coeff_good, PointFit_SortCosInc, PointFit_SortPsi, PointFit_SortLogD, boundary, n_param, n_pt_fit, n_1_fit, n_2_fit, SquaredScale):
def AnaAlongNumGradientLeapfrog_ThreeDetectors(q_pos_0, logL_0, dlogL_0, start_time, minimum_frequency, interferometers, randGen, epsilon, lengthT, scale, fit_coeff_good, PointFit_SortCosInc, PointFit_SortPsi, PointFit_SortLogD, boundary, n_param, n_pt_fit, n_1_fit, n_2_fit, SquaredScale, likelihood_gradient):
"""
"""
......@@ -904,7 +917,7 @@ def AnaAlongNumGradientLeapfrog_ThreeDetectors(q_pos_0, logL_0, dlogL_0, start_t
p_mom_0 = randGen.normal(loc=0, scale=1.0, size=9)
# Main HMC routine with numerical gradients
q_pos, p_mom, dlogL_pos, pt_fit_traj, dlogL_fit_traj, dlogL_fit_traj_ana, BreakIndic = AnaAlongNumGradientLeapfrog_ThreeDetectors_Trajectory(q_pos_0, p_mom_0, dlogL_0, start_time, minimum_frequency, interferometers, epsilon, lengthT, scale, boundary, fit_coeff_good, PointFit_SortCosInc, PointFit_SortPsi, PointFit_SortLogD, n_param, n_pt_fit, n_1_fit, n_2_fit, SquaredScale)
q_pos, p_mom, dlogL_pos, pt_fit_traj, dlogL_fit_traj, dlogL_fit_traj_ana, BreakIndic = AnaAlongNumGradientLeapfrog_ThreeDetectors_Trajectory(q_pos_0, p_mom_0, dlogL_0, start_time, minimum_frequency, interferometers, epsilon, lengthT, scale, boundary, fit_coeff_good, PointFit_SortCosInc, PointFit_SortPsi, PointFit_SortLogD, n_param, n_pt_fit, n_1_fit, n_2_fit, SquaredScale, likelihood_gradient)
# Compute the Log-likelihood at the end if there was no NaN
if BreakIndic==0:
......@@ -938,7 +951,7 @@ def AnaAlongNumGradientLeapfrog_ThreeDetectors(q_pos_0, logL_0, dlogL_0, start_t
return q_pos_0, logL_0, dlogL_0, pt_fit_traj, dlogL_fit_traj, dlogL_fit_traj_ana, Accept
def AnaAlongNumGradientLeapfrog_ThreeDetectors_Trajectory(q_pos_0, p_mom_0, dlogL_0, start_time, minimum_frequency, interferometers, epsilon, lengthT, scale, boundary, fit_coeff_good, PointFit_SortCosInc, PointFit_SortPsi, PointFit_SortLogD, n_param, n_pt_fit, n_1_fit, n_2_fit, SquaredScale):
def AnaAlongNumGradientLeapfrog_ThreeDetectors_Trajectory(q_pos_0, p_mom_0, dlogL_0, start_time, minimum_frequency, interferometers, epsilon, lengthT, scale, boundary, fit_coeff_good, PointFit_SortCosInc, PointFit_SortPsi, PointFit_SortLogD, n_param, n_pt_fit, n_1_fit, n_2_fit, SquaredScale, likelihood_gradient):
"""
HMC routine for a single trajectory where the gradients are all computed numerically.
......@@ -1005,9 +1018,13 @@ def AnaAlongNumGradientLeapfrog_ThreeDetectors_Trajectory(q_pos_0, p_mom_0, dlog
break
# NUMERICAL GRADIENT
parameters = paru.q_pos_to_dictionary(q_pos, start_time)
dlogL_pos, logL = bilby_wv.dlogL_ThreeDetectors(parameters, start_time, minimum_frequency, interferometers)
if likelihood_gradient is None:
parameters = paru.q_pos_to_dictionary(q_pos, start_time)
dlogL_pos, logL = bilby_wv.dlogL_ThreeDetectors(parameters, start_time, minimum_frequency, interferometers)
else:
likelihood_gradient.likelihood.parameters = paru.q_pos_to_dictionary(q_pos, start_time)
dlogL_pos = likelihood_gradient.calculate_dlogL_np()
logL = likelihood_gradient.likelihood.log_likelihood_ratio()
# ANALYTICAL gradient along the way
q_pos_fit = np.empty(9)
......@@ -1021,11 +1038,11 @@ def AnaAlongNumGradientLeapfrog_ThreeDetectors_Trajectory(q_pos_0, p_mom_0, dlog
q_pos_fit[7] = q_pos[7]
q_pos_fit[8] = np.log(q_pos[8])
dlogL_pos_ana = GradientLogL_AnalyticalFit_Cubic(q_pos_fit, fit_coeff_good, n_param)
if 0 is conf.param_indices:
if 0 in conf.param_indices:
dlogL_pos_ana[0] = Table_GradientFit(q_pos_fit, PointFit_SortCosInc, n_pt_fit, n_param, n_1_fit, n_2_fit, SquaredScale, 0)
if 2 is conf.param_indices:
if 2 in conf.param_indices:
dlogL_pos_ana[2] = Table_GradientFit(q_pos_fit, PointFit_SortPsi, n_pt_fit, n_param, n_1_fit, n_2_fit, SquaredScale, 2)
if 3 is conf.param_indices:
if 3 in conf.param_indices:
dlogL_pos_ana[3] = Table_GradientFit(q_pos_fit, PointFit_SortLogD, n_pt_fit, n_param, n_1_fit, n_2_fit, SquaredScale, 3)
......
......@@ -121,9 +121,9 @@ class GWTransientLikelihoodGradient(LikelihoodGradient):
self.likelihood.parameters = paru.parameters_new_dlnmu(self.likelihood.parameters_copy, -self.offset)
log_likelihood_minus = self.likelihood.log_likelihood_ratio()
parameters_plus = paru.parameters_new_dlnMc(self.likelihood.parameters_copy, self.offset)
parameters_plus = paru.parameters_new_dlnmu(self.likelihood.parameters_copy, self.offset)
# If parameters['symmetric_mass_ratio'] = 0.25: a forward step will give parameters_plus['symmetric_mass_ratio'] > 0.25 which is unphysical, hence we only use backward difference
if parameters_plus['symmetric_mass_ratio'] < 0.25:
if parameters_plus['symmetric_mass_ratio'] > 0.25:
self.log_likelihood_gradient[key] = (self.log_likelihood_at_point - log_likelihood_minus) / self.offset
else:
self.likelihood.parameters = parameters_plus
......@@ -166,4 +166,24 @@ class GWTransientLikelihoodGradient(LikelihoodGradient):
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
def convert_dlogL_dict_to_np_array(self, dlogL_dict):
fparam_names = ["cos(theta_jn)", "phi_c", "psi", "ln(D)", "ln(Mc)", "ln(mu)", "sin(dec)", "ra", "ln(tc)"]
dlogL_np = np.zeros(9)
# dlogL_np = np.zeros(len(fparam_names))
dict_keys = dlogL_dict.keys()
for i, key in enumerate(fparam_names):
if key in dict_keys:
dlogL_np[i] = dlogL_dict[key]
return dlogL_np
def calculate_dlogL_np(self, dlogL_dict=None):
if dlogL_dict is None:
dlogL_dict = self.calculate_log_likelihood_gradient()
return self.convert_dlogL_dict_to_np_array(dlogL_dict)
......@@ -1032,8 +1032,10 @@ logL on gpu = 4ms
- In ROQ paper, apparently the minimum Mchirp is 1.4 when it is equal to 1.18 for GW170817
- should I implement a method to compute `dlogL` which only takes as input `(logL, parameters)` where `logL` is a likelihood method?
- with the GW150914 ROQ example I got:
```
In [5]: %timeit likelihood.log_likelihood_ratio() 1.29 ms ± 8.67 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [5]: %timeit likelihood.log_likelihood_ratio()
1.29 ms ± 8.67 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
```
- be careful that in ... they use `interferometer.power_spectral_density_array` and `interferometer.frequency_domain_strain` which as I saw might reinterpolate each time the psd and multiplies by the mask for freq_array
- code properly what the scale factor should be
......@@ -1061,4 +1063,31 @@ logL on gpu = 4ms
instead of my `Compute_fISCO()` actually changes a bit the SNR and logL value
- => Try running phase1 with this configuration and see how the fitting works
AVANT LE DEPART EN DROME: AVOIR UN RUN PHASE 1AVEC IMRPHENOMPV2 ROQ sur GW170817 !
AVANT LE DEPART EN DROME: AVOIR UN RUN PHASE 1 AVEC IMRPHENOMPV2 ROQ sur GW170817 !
- ROQ run completely diverging...
- try to run with IMRPhenomPv2 but without ROQ
- try to run on a simpler case than GW170817, like GW150914 where the ROQ basis won't be rescaled etc...
- Why aren't my TaylorF2 runs working as well as before? That might be a clue why ROQ is not working as well...
- [DONE] change the self.offset for likelihood_gradient
- Am I going outside of the ROQ boundaries? Like chirp-mass etc ?
- In plot_traj.py, plot individual 0.5*p**2 to spot which parameter(s) is making the kinetic energy blow off
- it's definitely `ln(mu)`
- I noticed that in my dlogL_roq computation for `ln(mu)` I was using:
`paru.parameters_new_dlnMc(self.likelihood.parameters_copy, self.offset)`
instead of:
`paru.parameters_new_dlnmu(self.likelihood.parameters_copy, self.offset)`
- I had inverted the superior condition `parameters_plus['symmetric_mass_ratio'] > 0.25:` to inferior...
- After correcting these silly mistakes, my first ROQ trajectories get accepted. When plotting the Hamiltonian along the traj, it has a weird oscillatory behavior which seems highly correlated with the kinetic energies of phi_c and psi
- Can I run phase 1 with ROQ on a subset of parameters?
- in `calculate_log_likelihood_gradient()` I had forgotten to copy back the original values of the parameters:
`self.likelihood.parameters = self.likelihood.parameters_copy.copy()`
- Analyse `python plot_trajectory.py ../.output_data/GW170817/bilbyoriented/SUB-D_012345678/PSD_1/0_100_2000_200_200_0.005_dlogL1/ROQ/ --traj_nb=25` at step ~175 there is a clear bump in the Hamiltonian clearly due to a bump in dlogL/dsindec
- Analyse `python plot_trajectory.py ../.output_data/GW170817/bilbyoriented/SUB-D_012345678/PSD_1/0_15_2000_200_200_0.005_dlogL1/ROQ_scales_div_by_two/ --traj_nb=9`: there is a peak in dlogL/dlnmu
- Dividing the scales of phi_c and psi by two makes the Hamiltonian less oscillatory and improves the acceptance rate
IMPROVING CODE STRUCTURE:
- Let the approximant be an option `--approximant=`, and also part of the config.ini
- Create a sub-directory based on the approximant name
- Include the ROQ boundaries in my boundaries
- Use the `likelihood_gradient` class even for TaylorF2 so that then I can easily switch to an object oriented structure?
......@@ -39,6 +39,7 @@ if __name__=='__main__':
parser.add_option("--dlogL", default='dlogL1', action="store", type="string", help=""" Method to use to compute the gradient of the log-likelihood. Use `dlogL1` to use dlogL = <s|dh> - <h|dh>, `dlogL2` for dlogL = logL(p+offset) - logL(p-offset).""")
parser.add_option("--ifos", default='H1,L1,V1', action="store", type="string", help=""" Interferometers to run the analysis on. Default is 'H1,L1,V1' """)
parser.add_option("--param_offsets", default='707077777', action="store", type="string", help=""" -log10() of the offset used for central differencing on each parameter. Default is '707077777' meaning conf.param_offsets = [1e-07, 0, 1e-07, 0, 1e-07, 1e-07, 1e-07, 1e-07, 1e-07]. 0 Means analytical formula is used if available.""")
parser.add_option("--param_indices", default='012345678', action="store", type="string", help=""" Indices of the parameters to run the hmc on. Other are kept fixed. Default are values already stored in config.py file.""")
parser.add_option("--inj_file", default='../examples/massive_GW170817.ini', action="store", type="string", help="""Injection file path leading to the `.ini` file specifying the parameters of the injection.""")
......@@ -56,6 +57,8 @@ if __name__=='__main__':
conf.param_offsets[i] = 10**(-exponents[i])
conf.dlogL = opts.dlogL
conf.param_indices = [int(i) for i in list(opts.param_indices)]
inj_file_name = opts.inj_file.split('/')[-1]
inj_name = inj_file_name.split('.')[0]
......@@ -104,12 +107,12 @@ if __name__=='__main__':
# WAVEFORM GENERATION
injection_parameters['meta']['maximum_frequency_generated_waveform'] = injection_parameters['meta']['maximum_frequency_injected_waveform']
template_ifos = bilby_wv.WaveForm_ThreeDetectors(injection_parameters, minimum_frequency, interferometers)
template_ifos_injected = bilby_wv.WaveForm_ThreeDetectors(injection_parameters, minimum_frequency, interferometers)
# INJECT TEMPLATE INTO NOISE STRAIN DATA
# Add the template signal to each detector's strain data
for i in range(len(interferometers)):
interferometers[i].strain_data.frequency_domain_strain += template_ifos[i]
interferometers[i].strain_data.frequency_domain_strain += template_ifos_injected[i]
interferometers[i].fd_strain = interferometers[i].strain_data.frequency_domain_strain
......@@ -122,10 +125,6 @@ if __name__=='__main__':
if injection_parameters['meta']['approximant'] != 'IMRPhenomPv2':
raise ValueError("Need to set approximant to IMRPhenomPv2 to work with the ROQ basis, approximant here is {}".format(injection_parameters['meta']['approximant']))
duration_from_flow = set_inj.ComputeChirpTime3p5PN(minimum_frequency, injection_parameters['mass_1'], injection_parameters['mass_2'])
print('Duration of the signal from {}Hz is {} seconds.'.format(minimum_frequency, duration_from_flow))
scale_factor = injection_parameters['roq']['scale_factor']
roq_directory = injection_parameters['roq']['directory']
# Load in the pieces for the linear part of the ROQ. Note you will need to
......@@ -138,28 +137,6 @@ if __name__=='__main__':
weights_file_path = outdir + '{:.0f}Hz_{:.0f}Hz_{:.0f}Hz_{:.0f}s_weights.json'.format(params['flow'], injection_parameters['meta']['maximum_frequency_ifo'], injection_parameters['meta']['maximum_frequency_injected_waveform'], params['seglen'])
# injection_parameters = dict(
# mass_1=36.0, mass_2=29.0, a_1=0.4, a_2=0.3, tilt_1=0.0, tilt_2=0.0,
# phi_12=1.7, phi_jl=0.3, luminosity_distance=1000., theta_jn=0.4, psi=0.659,
# phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
# waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
# reference_frequency=20. * scale_factor,
# minimum_frequency=20. * scale_factor)
#
# waveform_generator = bilby.gw.WaveformGenerator(
# duration=duration, sampling_frequency=sampling_frequency,
# frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
# waveform_arguments=waveform_arguments,
# parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters)
# ifos = bilby.gw.detector.InterferometerList(['H1', 'L1', 'V1'])
# ifos.set_strain_data_from_power_spectral_densities(
# sampling_frequency=sampling_frequency, duration=duration,
# start_time=injection_parameters['geocent_time'] - 3)
# ifos.inject_signal(waveform_generator=waveform_generator,
# parameters=injection_parameters)
# make ROQ waveform generator
search_waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
duration=duration, sampling_frequency=sampling_frequency,
......@@ -213,13 +190,16 @@ if __name__=='__main__':
likelihood.save_weights(weights_file_path)
search_parameters_keys = conf.fparam_names
search_parameters_keys = [conf.fparam_names[i] for i in conf.param_indices]
likelihood.parameters = injection_parameters
print("")
print("likelihood.log_likelihood_ratio() = {}".format(likelihood.log_likelihood_ratio()))
likelihood_gradient = lg.GWTransientLikelihoodGradient(likelihood, search_parameters_keys)
dlogL = likelihood_gradient.calculate_log_likelihood_gradient()
pu.print_dict({"dlogL": dlogL})
dlogL_dict = likelihood_gradient.calculate_log_likelihood_gradient()
pu.print_dict({"dlogL_dict": dlogL_dict})
dlogL_np = likelihood_gradient.calculate_dlogL_np()
print("dlogL_np = {}".format(dlogL_np))
import IPython; IPython.embed()
......@@ -190,8 +190,9 @@ for param_index in param_indices:
ax111 = axs[1][1].twinx()
ax111.plot(steps, H_p_logL[1], linewidth=0.5, color='blue', label="0.5*p^2")
ax111.plot(steps, 0.5 * pmom_trajs[param_index]**2, linewidth=0.8, color='orange', label="0.5*p_param^2")
ax111.tick_params('y', colors='blue')
ax111.legend(loc=(0.8,0.8))
ax111.legend(loc=(0.8,0.7))
fig.tight_layout()
# import IPython; IPython.embed();
......
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