diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index 838209239..e7427074f 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -290,26 +290,26 @@ class HPolytope { return -1; } - //Nudge the Point p inside the Polytope void nudge_in(Point& p, NT tol=NT(0)) const { int m = A.rows(); const NT* b_data = b.data(); for (int i = 0; i < m; i++) { - //Check if corresponding hyperplane is violated - if (*b_data - A.row(i) * p.getCoefficients() < NT(-tol)){ - + + NT dist = *b_data - A.row(i) * p.getCoefficients(); + + if (dist < NT(-tol)){ //Nudging correction NT eps = -1e-7; - NT eps_1 = -(*b_data - A.row(i) * p.getCoefficients()); - MT A_nor = A; - A_nor.normalize(); + NT eps_1 = -dist; + //A.row is already normalized, no need to do it again + VT A_i = A.row(i); NT eps_2 = eps_1 + eps; //Nudge the point inside with respect to the normal its vector - Point shift(A_nor.row(i)); + Point shift(A_i); shift.operator*=(eps_2); p.operator+=(shift); } diff --git a/include/random_walks/gaussian_hamiltonian_monte_carlo_exact_walk.hpp b/include/random_walks/gaussian_hamiltonian_monte_carlo_exact_walk.hpp index 49a46f532..b874258db 100644 --- a/include/random_walks/gaussian_hamiltonian_monte_carlo_exact_walk.hpp +++ b/include/random_walks/gaussian_hamiltonian_monte_carlo_exact_walk.hpp @@ -89,6 +89,9 @@ struct Walk unsigned int n = P.dimension(); NT T; + //normalize the Polyope + P.normalize(); + for (auto j=0u; j::apply(n, rng, false); + //normalize the Polyope + P.normalize(); + NT T = rng.sample_urdist() * _Len; int it = 0; @@ -220,6 +225,7 @@ private : } _lambda_prev = pbpair.first; update_position(_p, _v, _lambda_prev, _omega); + P.nudge_in(_p); T -= _lambda_prev; P.compute_reflection(_v, _p, pbpair.second); it++;