Implementation of cb algorithm with H-polytopes in MMC for zonotopes (#1

)

* correct input list algo in Rcpp volume.cpp

* Improvements in rounding and rotating R functions (GeomScale#76)

* add seed in rotate_polytope and update Rd files

* update R rounding function to return the linear map

* fix rounding in V-poly

* overloading rotating() function and remove unused comments

* implement zonotope volume computation with cb algorithm and hpoly in MMC

* fix all bugs in new cooling_hpoly algorithm and improve test function

* all tests for zonotopes with c algortihm passed successfully

R-proj/R/RcppExports.R

@@ -144,7 +144,7 @@ inner_ball <- function(P) {
 #' @param Zono_gen A boolean parameter to declare if the requested polytope has to be a zonotope.
 #' @param dim_gen An integer to declare the dimension of the requested polytope.
 #' @param m_gen An integer to declare the number of generators for the requested random zonotope or the number of vertices for a V-polytope.
-#' @param seed This variable can be used to set a seed for the random polytope generator.
+#' @param seed Optional. A fixed seed for the random polytope generator.
 #'
 #' @section warning:
 #' Do not use this function.
@@ -158,13 +158,14 @@ poly_gen <- function(kind_gen, Vpoly_gen, Zono_gen, dim_gen, m_gen, seed = NULL)
 #'
 #' @param P A convex polytope (H-, V-polytope or a zonotope).
 #' @param T Optional. A rotation matrix.
+#' @param seed Optional. A fixed seed for the random linear map generator.
 #'
 #' @section warning:
 #' Do not use this function.
 #'
 #' @return A matrix that describes the rotated polytope
-rotating <- function(P, T = NULL) {
-    .Call(`_volesti_rotating`, P, T)
+rotating <- function(P, T = NULL, seed = NULL) {
+    .Call(`_volesti_rotating`, P, T, seed)
 }
 
 #' Internal rcpp function for the rounding of a convex polytope

R-proj/R/rotate_polytope.R

@@ -4,6 +4,7 @@
 #' 
 #' @param P A convex polytope. It is an object from class (a) Hpolytope, (b) Vpolytope, (c) Zonotope, (d) intersection of two V-polytopes.
 #' @param T Optional. A \eqn{d\times d} rotation matrix.
+#' @param seed Optional. A fixed seed for the random linear map generator.
 #' 
 #' @return A list that contains the rotated polytope and the matrix \eqn{T} of the linear transformation.
 #'
@@ -27,10 +28,10 @@
 #' Z = gen_rand_zonotope(3,6)
 #' poly_matrix_list = rotate_polytope(Z)
 #' @export
-rotate_polytope <- function(P, T = NULL){
+rotate_polytope <- function(P, T = NULL, seed = NULL){
 
   #call rcpp rotating function
-  Mat = rotating(P, T)
+  Mat = rotating(P, T, seed)
 
   n = P$dimension
   m=dim(Mat)[2]-n
@@ -43,7 +44,7 @@ rotate_polytope <- function(P, T = NULL){
 
   # remove first column
   A = Mat[,-c(1)]
-  #A = Mat[,-c(1)]
+
   type = P$type
   if (type == 2) {
     PP = Vpolytope$new(A)

R-proj/R/round_polytope.R

@@ -36,13 +36,14 @@ round_polytope <- function(P){
 
   # remove first column
   A = Mat[,-c(1)]
+
   type = P$type
   if (type == 2) {
-    PP = list("P" = Vpolytope$new(A), "round_value" = ret_list$round_value)
+    PP = list("P" = Vpolytope$new(A), "T" = ret_list$T, "shift" = ret_list$shift, "round_value" = ret_list$round_value)
   }else if (type == 3) {
-    PP = list("P" = Zonotope$new(A), "round_value" = ret_list$round_value)
+    PP = list("P" = Zonotope$new(A), "T" = ret_list$T, "shift" = ret_list$shift, "round_value" = ret_list$round_value)
   } else {
-    PP = list("P" = Hpolytope$new(A,b), "round_value" = ret_list$round_value)
+    PP = list("P" = Hpolytope$new(A,b), "T" = ret_list$T, "shift" = ret_list$shift, "round_value" = ret_list$round_value)
   }
   return(PP)
 }

R-proj/src/poly_gen.cpp

@@ -31,7 +31,7 @@
 //' @param Zono_gen A boolean parameter to declare if the requested polytope has to be a zonotope.
 //' @param dim_gen An integer to declare the dimension of the requested polytope.
 //' @param m_gen An integer to declare the number of generators for the requested random zonotope or the number of vertices for a V-polytope.
-//' @param seed This variable can be used to set a seed for the random polytope generator.
+//' @param seed Optional. A fixed seed for the random polytope generator.
 //'
 //' @section warning:
 //' Do not use this function.

R-proj/src/rotating.cpp

@@ -27,13 +27,15 @@
 //'
 //' @param P A convex polytope (H-, V-polytope or a zonotope).
 //' @param T Optional. A rotation matrix.
+//' @param seed Optional. A fixed seed for the random linear map generator.
 //'
 //' @section warning:
 //' Do not use this function.
 //'
 //' @return A matrix that describes the rotated polytope
 // [[Rcpp::export]]
-Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable<Rcpp::NumericMatrix> T = R_NilValue){
+Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable<Rcpp::NumericMatrix> T = R_NilValue,
+                              Rcpp::Nullable<int> seed = R_NilValue){
 
     typedef double NT;
     typedef Cartesian<NT>    Kernel;
@@ -51,6 +53,8 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable<Rcpp::NumericMat
     unsigned int n = P.field("dimension");
     int type = P.field("type");
 
+    int seed2 = (!seed.isNotNull()) ? std::chrono::system_clock::now().time_since_epoch().count() : Rcpp::as<int>(seed);
+
     switch (type) {
         case 1: {
             // Hpolytope
@@ -60,7 +64,7 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable<Rcpp::NumericMat
                 TransorfMat = Rcpp::as<MT>(T);
                 HP.linear_transformIt(TransorfMat.inverse());
             } else {
-                TransorfMat = rotating < MT > (HP);
+                TransorfMat = rotating < MT > (HP, seed2);
             }
             Mat = extractMatPoly(HP);
             break;
@@ -73,7 +77,7 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable<Rcpp::NumericMat
                 TransorfMat = Rcpp::as<MT>(T);
                 VP.linear_transformIt(TransorfMat.inverse());
             } else {
-                TransorfMat = rotating < MT > (VP);
+                TransorfMat = rotating < MT > (VP, seed2);
             }
             Mat = extractMatPoly(VP);
             break;
@@ -86,24 +90,13 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable<Rcpp::NumericMat
                 TransorfMat = Rcpp::as<MT>(T);
                 ZP.linear_transformIt(TransorfMat.inverse());
             } else {
-                TransorfMat = rotating < MT > (ZP);
+                TransorfMat = rotating < MT > (ZP, seed2);
             }
             Mat = extractMatPoly(ZP);
             break;
         }
         case 4: {
-            Vpolytope VP1;
-            Vpolytope VP2;
-            InterVP VPcVP;
-            VP1.init(n, Rcpp::as<MT>(P.field("V1")), VT::Ones(Rcpp::as<MT>(P.field("V1")).rows()));
-            VP2.init(n, Rcpp::as<MT>(P.field("V2")), VT::Ones(Rcpp::as<MT>(P.field("V2")).rows()));
-            VPcVP.init(VP1, VP2);
-            if (T.isNotNull()) {
-                TransorfMat = Rcpp::as<MT>(T);
-                VPcVP.linear_transformIt(TransorfMat.inverse());
-            } else {
-                TransorfMat = rotating < MT > (VPcVP);
-            }
+            throw Rcpp::exception("volesti does not support rotation for this representation currently.");
         }
     }
 

R-proj/src/rounding.cpp

@@ -97,24 +97,7 @@ Rcpp::List rounding (Rcpp::Reference P){
         }
         case 4: {
             throw Rcpp::exception("volesti does not support rounding for this representation currently.");
-            /*
-            Vpolytope VP1;
-            Vpolytope VP2;
-            VP1.init(n, Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("V1")),
-                     VT::Ones(Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("V1")).rows()));
-            VP2.init(n, Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("V2")),
-                     VT::Ones(Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("V2")).rows()));
-            VPcVP.init(VP1, VP2);
-
-            if (!VPcVP.is_feasible()) throw Rcpp::exception("Empty set!");
-            InnerBall = VPcVP.ComputeInnerBall();
-
-            diam = 2.0 * std::sqrt(NT(n)) * InnerBall.second;
-            VPcVP.comp_diam(diam);
-            delta = 4.0 * InnerBall.second / std::sqrt(NT(n));*/
-
         }
-        //default: throw Rcpp::exception("Wrong polytope input");
     }
 
     unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
@@ -126,26 +109,28 @@ Rcpp::List rounding (Rcpp::Reference P){
     // initialization
     vars<NT, RNGType> var(rnum,n,walkL,1,0.0,0.0,0,0.0,0,InnerBall.second,diam,rng,urdist,urdist1,
                           delta,verbose,rand_only,false,NN,birk,ball_walk,cdhr,rdhr,billiard);
-    std::pair <NT, NT> round_res;
+    std::pair< std::pair<MT, VT>, NT > round_res;
 
     switch (type) {
         case 1: {
-            round_res = rounding_min_ellipsoid(HP, InnerBall, var);
+            round_res = rounding_min_ellipsoid<MT, VT>(HP, InnerBall, var);
             Mat = extractMatPoly(HP);
             break;
         }
         case 2: {
-            round_res = rounding_min_ellipsoid(VP, InnerBall, var);
+            round_res = rounding_min_ellipsoid<MT, VT>(VP, InnerBall, var);
             Mat = extractMatPoly(VP);
             break;
         }
         case 3: {
-            round_res = rounding_min_ellipsoid(ZP, InnerBall, var);
+            round_res = rounding_min_ellipsoid<MT, VT>(ZP, InnerBall, var);
             Mat = extractMatPoly(ZP);
             break;
         }
     }
 
-    return Rcpp::List::create(Rcpp::Named("Mat") = Mat , Rcpp::Named("round_value") = round_res.first);
+    return Rcpp::List::create(Rcpp::Named("Mat") = Mat, Rcpp::Named("T") = Rcpp::wrap(round_res.first.first),
+                              Rcpp::Named("shift") = Rcpp::wrap(round_res.first.second),
+                              Rcpp::Named("round_value") = round_res.second);
 
 }

R-proj/src/volume.cpp

@@ -337,8 +337,8 @@ double volume (Rcpp::Reference P,
         alpha = Rcpp::as<NT>(Rcpp::as<Rcpp::List>(settings)["alpha"]);
         cb_params++;
     }
-    if (Rcpp::as<Rcpp::List>(algo).containsElementNamed("L")) {
-        diam = Rcpp::as<NT>(Rcpp::as<Rcpp::List>(algo)["L"]);
+    if (Rcpp::as<Rcpp::List>(settings).containsElementNamed("L")) {
+        diam = Rcpp::as<NT>(Rcpp::as<Rcpp::List>(settings)["L"]);
     }
 
     if ((CB && cg_params > 0) || (CG && cb_params > 0) || (!CB & !CG && (cg_params > 0 || cb_params > 0))){

include/convex_bodies/ball.h

@@ -23,11 +23,17 @@ struct Ball{
 
     Ball(Point cc, NT RR) : c(cc),	 R(RR) {}
 
-    NT ComputeDiameter() const
+    NT ComputeDiameter()
     {
         return std::sqrt(R) * 2;
     }
 
+    void set_diameter(const NT &diam) {}
+
+    NT get_diameter() const {
+        return std::sqrt(R) * 2;
+    }
+
     std::pair<Point,NT> InnerBall() const
     {
         return std::pair<Point,NT>(c, R);

include/convex_bodies/ballintersectconvex.h

@@ -21,6 +21,7 @@ class BallIntersectPolytope {
     typedef typename CBall::NT NT;
     typedef typename CBall::PointType PointType;
     typedef Eigen::Matrix<NT,Eigen::Dynamic,1> VT;
+    NT diameter;
 
     BallIntersectPolytope() {}
 
@@ -34,9 +35,22 @@ class BallIntersectPolytope {
         return P.InnerBall();
     }
 
+    //NT ComputeDiameter() const
+   // {
+     //   return NT(2) * B.radius();
+    //}
+
     NT ComputeDiameter() const
     {
-        return NT(2) * B.radius();
+        return 2.0 * B.radius();
+    }
+
+    void set_diameter(const NT &diam) {
+        diameter = diam;
+    }
+
+    NT get_diameter() const {
+        return diameter;
     }
 
     int is_in(PointType &p) const
@@ -54,7 +68,7 @@ class BallIntersectPolytope {
         return P.dimension();
     }
 
-    void comp_diam(NT &diam) {
+    void comp_diam(NT &diam) const {
         diam = 2.0 * B.radius();
     }
 

include/convex_bodies/hpolytope.h

@@ -36,6 +36,7 @@ class HPolytope{
     VT b; // vector b, s.t.: Ax<=b
     unsigned int            _d; //dimension
     std::pair<Point,NT> _inner_ball;
+    NT diameter;
     //NT maxNT = 1.79769e+308;
     //NT minNT = -1.79769e+308;
     NT maxNT = std::numeric_limits<NT>::max();
@@ -67,6 +68,14 @@ class HPolytope{
         return NT(4) * std::sqrt(NT(_d)) * _inner_ball.second;
     }
 
+    void set_diameter(const NT &diam) {
+        diameter = diam;
+    }
+
+    NT get_diameter() const {
+        return diameter;
+    }
+
     // return dimension
     unsigned int dimension() const {
         return _d;
@@ -111,11 +120,11 @@ class HPolytope{
         return 0.0;
     }
 
-    void comp_diam(NT &diam) {
+    void comp_diam(NT &diam) const {
         diam = 4.0 * std::sqrt(NT(_d)) * ComputeInnerBall().second;
     }
 
-    void comp_diam(NT &diam, const NT &cheb_rad) {
+    void comp_diam(NT &diam, const NT &cheb_rad) const {
         diam = 4.0 * std::sqrt(NT(_d)) * cheb_rad;
     }