Gaussian HMC Walk Polytope Normalization

Author: vgnecula

include/convex_bodies/hpolytope.h

@@ -56,6 +56,7 @@ class HPolytope {
     MT                   A; //matrix A
     VT                   b; // vector b, s.t.: Ax<=b
     std::pair<Point, NT> _inner_ball;
+    bool normalized = 0;
 
 public:
     //TODO: the default implementation of the Big3 should be ok. Recheck.
@@ -883,15 +884,24 @@ class HPolytope {
 
     void normalize()
     {
-        NT row_norm;
-        for (int i = 0; i < num_of_hyperplanes(); ++i)
+        if(!is_normalized())
         {
-            row_norm = A.row(i).norm();
-            A.row(i) = A.row(i) / row_norm;
-            b(i) = b(i) / row_norm;
+            NT row_norm;
+            for (int i = 0; i < num_of_hyperplanes(); ++i)
+            {
+                row_norm = A.row(i).norm();
+                A.row(i) = A.row(i) / row_norm;
+                b(i) = b(i) / row_norm;
+            }
+            normalized = true;
         }
     }
 
+    bool is_normalized()
+    {
+        return normalized;
+    }
+
     void compute_reflection(Point& v, Point const&, int const& facet) const
     {
         v += -2 * v.dot(A.row(facet)) * A.row(facet);

include/random_walks/gaussian_hamiltonian_monte_carlo_exact_walk.hpp

@@ -55,8 +55,9 @@ struct Walk
     typedef typename Polytope::VT VT;
 
     template <typename GenericPolytope>
-    Walk(GenericPolytope &P, Point const& p, NT const& a_i, RandomNumberGenerator &rng)
+    Walk(GenericPolytope& P, Point const& p, NT const& a_i, RandomNumberGenerator &rng)
     {
+        P.normalize();
         _Len = compute_diameter<GenericPolytope>
                 ::template compute<NT>(P);
         _omega = std::sqrt(NT(2) * a_i);
@@ -65,9 +66,10 @@ struct Walk
     }
 
     template <typename GenericPolytope>
-    Walk(GenericPolytope &P, Point const& p, NT const& a_i, RandomNumberGenerator &rng,
+    Walk(GenericPolytope& P, Point const& p, NT const& a_i, RandomNumberGenerator &rng,
          parameters const& params)
     {
+        P.normalize();
         _Len = params.set_L ? params.m_L
                           : compute_diameter<GenericPolytope>
                             ::template compute<NT>(P);
@@ -80,18 +82,16 @@ struct Walk
     <
         typename GenericPolytope
     >
-    inline void apply(GenericPolytope& P,
+    inline void apply(GenericPolytope const& P,
                       Point& p,
                       NT const& a_i,
                       unsigned int const& walk_length,
                       RandomNumberGenerator &rng)
     {
+
         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;
@@ -124,7 +124,7 @@ struct Walk
     <
         typename GenericPolytope
     >
-    inline void get_starting_point(GenericPolytope& P,
+    inline void get_starting_point(GenericPolytope const& P,
                                    Point const& center,
                                    Point &q,
                                    unsigned int const& walk_length,
@@ -145,7 +145,7 @@ struct Walk
     <
         typename GenericPolytope
     >
-    inline void parameters_burnin(GenericPolytope& P,
+    inline void parameters_burnin(GenericPolytope const& P,
                                   Point const& center,
                                   unsigned int const& num_points,
                                   unsigned int const& walk_length,
@@ -195,7 +195,7 @@ private :
     <
         typename GenericPolytope
     >
-    inline void initialize(GenericPolytope& P,
+    inline void initialize(GenericPolytope const& P,
                            Point const& p,
                            NT const& a_i,
                            RandomNumberGenerator &rng)
@@ -205,9 +205,6 @@ private :
         _p = p;
         _v = GetDirection<Point>::apply(n, rng, false);
 
-        //normalize the Polyope
-        P.normalize();
-
         NT T = rng.sample_urdist() * _Len;
         int it = 0;