1function [O, logW] = extend_orbit_forward(logmu, gradlogmu, theta_b, rho_b, logw_b, h, delta, L)
2 O = {};
3 logW = [];
4 theta = theta_b;
5 rho = rho_b;
6 logw = logw_b;
8 for iter = 1:L
9 ell_f = micro(logmu, gradlogmu, theta, rho, h, delta);
10 ell = p_micro(ell_f);
11 [theta_1, rho_1] = leapfrog(logmu, gradlogmu, theta, rho, h * (2^(-ell)), 2^ell);
12 ell_b = micro(logmu, gradlogmu, theta_1, -rho_1, h, delta);
14 if pmf_p_micro(ell, ell_b) == 0
15 logw = -Inf;
16 else
17 logw = logw + logmu(theta_1) - 0.5 * norm(rho_1)^2 - logmu(theta) + 0.5 * norm(rho)^2 + log(pmf_p_micro(ell, ell_b))- log(pmf_p_micro(ell, ell_f));
18 end
20 theta = theta_1;
21 rho = rho_1;
22 O = [O; {theta, rho}];
23 logW = [logW; logw];
24 end
25end