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

plot_trajectory from the real start with correct Acceptance proba

parent 859e93d9
......@@ -679,7 +679,7 @@ if __name__ == '__main__':
else:
cut.logger.info(" There are no trajectories left to run in phase1.")
# sys.exit()
sys.exit()
#************************************************************************************************/
#***************************** Phase 2 : Fitting routines *************************************/
......
......@@ -32,6 +32,7 @@ def NumericalGradientLeapfrog_ThreeDetectors(q_pos_0, logL_0, dlogL_0, randGen,
# Main HMC routine with numerical gradients
q_pos, p_mom, dlogL_pos, pt_fit_traj, dlogL_fit_traj, BreakIndic, H_p_logL_logP_traj, pmom_traj, dlogTrajPi_traj = NumericalGradientLeapfrog_ThreeDetectors_Trajectory(q_pos_0, p_mom_0, dlogL_0, logL_0, epsilon, lengthT, traj_index, sampler)
# Compute the Log-likelihood at the end if there was no NaN
if BreakIndic==0:
# Hamiltonian at the start of the trajectory
......@@ -39,6 +40,19 @@ def NumericalGradientLeapfrog_ThreeDetectors(q_pos_0, logL_0, dlogL_0, randGen,
logTrajPi_0 = sampler.traj_prior_gradients.traj_priors.ln_prob(qpos_dict_0)
H_0 = computeHamiltonian(p_mom_0, logL_0, logTrajPi_0)
# Code snippet only used when the `--plot_traj` option was set for the run.
# It records data for each trajectory in a `/plot_traj` specific folder which
# can then be plotted using the /Tests/plot_trajectory.py script.
# Here we simply add the starting point of the trajectory.
if RUN_CONST.plot_traj:
H_p_logL_logP_0 = [H_0, 0.5 * (p_mom_0**2).sum(), logL_0, logTrajPi_0]
sampler.dict_of_trajectories['H_p_logL_logP'].append(H_p_logL_logP_0)
sampler.dict_of_trajectories['pmom_trajs'].append(p_mom_0.tolist())
sampler.dict_of_trajectories['qpos'].append(q_pos_0.tolist())
sampler.dict_of_trajectories['dlogL'].append(dlogL_0)
dlogTrajPi_0 = sampler.traj_prior_gradients.calculate_dlogPi_np(qpos_dict_0)
sampler.dict_of_trajectories['dlogTrajPi_traj'].append(dlogTrajPi_0.tolist())
# Hamiltonian at the end of the trajectory
# sampler.update_likelihood_parameters_from_q_pos(q_pos)
qpos_dict_final = sampler.q_pos_to_parameters_dict(q_pos)
......@@ -60,13 +74,6 @@ def NumericalGradientLeapfrog_ThreeDetectors(q_pos_0, logL_0, dlogL_0, randGen,
logL_0 = logL_final # update initial log-likelihood to final log-likelihood
Accept = 1 # set acceptation flag to 1
# Record both the positions in matrix pt[][] and the value of the likelihood gradient in dlogL[][] for future fitting
# if flag_update_qpos_global_fit:
# for i in range(len(pt_fit_traj)):
# qpos_global_fit.append(pt_fit_traj[i])
# dlogL_global_fit.append(dlogL_fit_traj[i])
# oluts_fit.update_oluts(pt_fit_traj, dlogL_fit_traj)
else: # Reject position
Accept = 0
else: # Reject position because of divergence and/or NaN
......@@ -74,7 +81,18 @@ def NumericalGradientLeapfrog_ThreeDetectors(q_pos_0, logL_0, dlogL_0, randGen,
proba_acc = 0
a = 1
# Code snippet only used when the `--plot_traj` option was set for the run.
# It records data for each trajectory in a `/plot_traj` specific folder which
# can then be plotted using the /Tests/plot_trajectory.py script
if RUN_CONST.plot_traj:
if BreakIndic!=0:
H_p_logL_logP_0 = [0, 0, 0, 0]
sampler.dict_of_trajectories['H_p_logL_logP'].append(0)
sampler.dict_of_trajectories['pmom_trajs'].append([0]*sampler.n_dim)
sampler.dict_of_trajectories['qpos'].append([0]*sampler.n_dim)
sampler.dict_of_trajectories['dlogL'].append(0)
sampler.dict_of_trajectories['dlogTrajPi_traj'].append([0]*sampler.n_dim)
for i in range(len(pt_fit_traj)):
sampler.dict_of_trajectories['H_p_logL_logP'].append(H_p_logL_logP_traj[i])
sampler.dict_of_trajectories['pmom_trajs'].append(pmom_traj[i])
......
......@@ -106,8 +106,8 @@ else:
# q_pos = pt_fit_traj[search_parameter_indices]
# dlogL = dlogL_all[search_parameter_indices]
# import IPython; IPython.embed();
idx_start = trajectory_index * lengthT
idx_end = (trajectory_index +1) * lengthT
idx_start = trajectory_index * (lengthT + 1)
idx_end = (trajectory_index +1) * (lengthT + 1)
pt_fit_traj = pt_fit_traj[:, idx_start: idx_end]
dlogL = dlogL[:, idx_start: idx_end]
H_p_logL_logP = H_p_logL_logP[:, idx_start: idx_end]
......
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