Complexity improvements and Polytope Normalization

GeomScale · Jun 14, 2024 · 6dc3108 · 6dc3108
1 parent c623e6f
commit 6dc3108
Show file tree

Hide file tree

Showing 2 changed files with 16 additions and 10 deletions.
diff --git 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
@@ -89,6 +89,9 @@ struct Walk
         unsigned int n = P.dimension();
         NT T;
 
+        //normalize the Polyope
+        P.normalize();
+
         for (auto j=0u; j<walk_length; ++j)
         {
             T = rng.sample_urdist() * _Len;
@@ -105,8 +108,7 @@ struct Walk
                 _lambda_prev = pbpair.first;
                 T -= _lambda_prev;
                 update_position(_p, _v, _lambda_prev, _omega);
-                if(P.is_in(_p) != -1)
-                    P.nudge_in(_p);
+                P.nudge_in(_p);
                 P.compute_reflection(_v, _p, pbpair.second);
                 it++;
             }
@@ -203,6 +205,9 @@ private :
         _p = p;
         _v = GetDirection<Point>::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++;