new volume structure SoB plus random walks

GeomScale · Mar 28, 2020 · ddcbfeb · ddcbfeb
1 parent 8f24239
commit ddcbfeb
Show file tree

Hide file tree

Showing 5 changed files with 696 additions and 46 deletions.
diff --git a/include/convex_bodies/ballintersectconvex.h b/include/convex_bodies/ballintersectconvex.h
@@ -27,6 +27,16 @@ class BallIntersectPolytope {
     Polytope first() const { return P; }
     CBall second() const { return B; }
 
+    std::pair<Point,NT> InnerBall()
+    {
+        return P.InnerBall();
+    }
+
+    NT ComputeDiameter()
+    {
+        return NT(4) * B.radius();
+    }
+
     int is_in(Point &p){
         if(B.is_in(p)==-1)
             return P.is_in(p);

diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h
@@ -33,40 +33,38 @@ class HPolytope{
     MT A; //matrix A
     VT b; // vector b, s.t.: Ax<=b
     unsigned int            _d; //dimension
+    std::pair<Point,NT> inner_ball;
     //NT maxNT = 1.79769e+308;
     //NT minNT = -1.79769e+308;
     NT maxNT = std::numeric_limits<NT>::max();
     NT minNT = std::numeric_limits<NT>::lowest();
 
 public:
-    HPolytope() {}
-
-    // constructor: cube(d)
-    HPolytope(unsigned int d): _d(d) {
-        A.resize(2 * d, d);
-        b.resize(2 * d);
-        for (unsigned int i = 0; i < d; ++i) {
-            b(i) = 1;
-            for (unsigned int j = 0; j < d; ++j) {
-                if (i == j) {
-                    A(i, j) = 1;
-                } else {
-                    A(i, j) = 0;
-                }
-            }
-        }
-        for (unsigned int i = 0; i < d; ++i) {
-            b(i + d) = 1;
-            for (unsigned int j = 0; j < d; ++j) {
-                if (i == j) {
-                    A(i + d, j) = -1;
-                } else {
-                    A(i + d, j) = 0;
-                }
-            }
-        }
+
+    HPolytope()
+    {
+        typedef Point PolytopePoint;
+        typedef typename Point::FT NT;
+        inner_ball = ComputeChebychevBall<NT, Point>(A, b);
+    }
+
+    std::pair<Point,NT> InnerBall()
+    {
+        return inner_ball;
+    }
+
+    //Compute Chebyshev ball of H-polytope P:= Ax<=b
+    //Use LpSolve library
+    std::pair<Point,NT> ComputeInnerBall()
+    {
+        inner_ball = ComputeChebychevBall<NT, Point>(A, b);
+        return inner_ball;
     }
 
+    NT ComputeDiameter()
+    {
+        return NT(4) * std::sqrt(NT(_d)) * inner_ball.second;
+    }
 
     // return dimension
     unsigned int dimension() const {
@@ -263,15 +261,6 @@ class HPolytope{
     }
 
 
-    //Compute Chebyshev ball of H-polytope P:= Ax<=b
-    //Use LpSolve library
-    std::pair<Point,NT> ComputeInnerBall() {
-
-        //lpSolve lib for the linear program
-        return ComputeChebychevBall<NT, Point>(A, b);
-    }
-
-
     // compute intersection point of ray starting from r and pointing to v
     // with polytope discribed by A and b
     std::pair<NT,NT> line_intersect(Point &r, Point &v) {
@@ -308,8 +297,11 @@ class HPolytope{
 
     // compute intersection points of a ray starting from r and pointing to v
     // with polytope discribed by A and b
-    std::pair<NT,NT> line_intersect(Point &r, Point &v, std::vector<NT> &Ar,
-            std::vector<NT> &Av, bool pos = false) {
+    std::pair<NT,NT> line_intersect(Point &r,
+                                    Point &v,
+                                    std::vector<NT> &Ar,
+                                    std::vector<NT> &Av,
+                                    bool pos = false) {
 
         NT lamda = 0, min_plus = NT(maxNT), max_minus = NT(minNT);
         NT sum_nom, sum_denom, mult;
@@ -346,8 +338,12 @@ class HPolytope{
         return std::pair<NT, NT>(min_plus, max_minus);
     }
 
-    std::pair<NT,NT> line_intersect(Point &r, Point &v, std::vector<NT> &Ar,
-            std::vector<NT> &Av, const NT &lambda_prev, bool pos = false) {
+    std::pair<NT,NT> line_intersect(Point &r,
+                                    Point &v,
+                                    std::vector<NT> &Ar,
+                                    std::vector<NT> &Av,
+                                    const NT &lambda_prev,
+                                    bool pos = false) {
 
         NT lamda = 0, min_plus = NT(maxNT), max_minus = NT(minNT);
         NT sum_nom, sum_denom, mult;

diff --git a/include/samplers/samplers.h b/include/samplers/samplers.h
@@ -236,13 +236,13 @@ void rand_point_generator(Polytope &P,
 
 template <typename BallPoly, typename PointList, typename Parameters, typename Point>
 void rand_point_generator(BallPoly &PBLarge,
-                         Point &p,   // a point to start
-                         const unsigned int rnum,
-                         const unsigned int walk_len,
-                         PointList &randPoints,
-                         const BallPoly &PBSmall,
-                         unsigned int &nump_PBSmall,
-                         const Parameters &var) {  // constants for volume
+                          Point &p,   // a point to start
+                          const unsigned int rnum,
+                          const unsigned int walk_len,
+                          PointList &randPoints,
+                          const BallPoly &PBSmall,
+                          unsigned int &nump_PBSmall,
+                          const Parameters &var) {  // constants for volume
 
     typedef typename Point::FT NT;
     typedef typename Parameters::RNGType RNGType;