Refactoring polytopes (#82)

* Eliminating "init" and refactoring of the constructors and introducing std::make_pair * Add implementation for the Big3 * Eliminate calls to Hpolytope.init * Eliminate free_them_all function call * fix R tests and improve c++ polytope generators * fix errors in R interface and cross polytope generator * Polishing according to the review in the PR. * fix bugs on zonotope destructor * remove Chebychev ball computation from the constructor Co-authored-by: TSIGARIDAS Elias <[email protected]> Co-authored-by: Tolis <[email protected]>
GeomScale · Sep 29, 2020 · e547120 · e547120
1 parent 6b07cf5
commit e547120
Show file tree

Hide file tree

Showing 54 changed files with 843 additions and 760 deletions.
diff --git a/R-proj/R/RcppExports.R b/R-proj/R/RcppExports.R
@@ -298,7 +298,7 @@ rounding <- function(P, method = NULL, seed = NULL) {
 #' \dQuote{Billiard walk - a new sampling algorithm for control and optimization,} \emph{IFAC Proceedings Volumes,} 2014.},
 #'
 #' @references \cite{Y. Chen, R. Dwivedi, M. J. Wainwright and B. Yu,
-#' \dQuote{Vaidya walk: A sampling algorithm based on the volumetric barrier,} \emph{55th Annual Allerton Conference on Communication, Control, and Computing,} 2017.}
+#' \dQuote{Fast MCMC Sampling Algorithms on Polytopes,} \emph{Journal of Machine Learning Research,} 2018.}
 #'
 #' @return A \eqn{d\times n} matrix that contains, column-wise, the sampled points from the convex polytope P.
 #' @examples

diff --git a/R-proj/man/compute_indicators.Rd b/R-proj/man/compute_indicators.Rd
diff --git a/R-proj/man/copula.Rd b/R-proj/man/copula.Rd
diff --git a/R-proj/man/sample_points.Rd b/R-proj/man/sample_points.Rd
diff --git a/R-proj/man/writeSdpaFormatFile.Rd b/R-proj/man/writeSdpaFormatFile.Rd
diff --git a/R-proj/man/zonotope_approximation.Rd b/R-proj/man/zonotope_approximation.Rd
diff --git a/R-proj/src/Makevars b/R-proj/src/Makevars
@@ -1,4 +1,4 @@
-PKG_CPPFLAGS=-I../../external/boost -I../../external/LPsolve_src/run_headers -I../../external/minimum_ellipsoid -I../../include
+PKG_CPPFLAGS= -I../../external/boost -I../../external/LPsolve_src/run_headers -I../../external/minimum_ellipsoid -I../../include
 
 PKG_CXXFLAGS= -lm -ldl -Wno-ignored-attributes -DBOOST_NO_AUTO_PTR -DDISABLE_NLP_ORACLES
 

diff --git a/R-proj/src/exact_vol.cpp b/R-proj/src/exact_vol.cpp
@@ -88,11 +88,10 @@ double exact_vol(Rcpp::Nullable<Rcpp::Reference> P) {
 
     } else if (type == 3) {
 
-        typedef Zonotope <Point> zonotope;
-        zonotope ZP;
+        typedef Zonotope<Point> zonotope;
         dim = Rcpp::as<Rcpp::Reference>(P).field("dimension");
 
-        ZP.init(dim, Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("G")),
+        zonotope ZP(dim, Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("G")),
                 VT::Ones(Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("G")).rows()));
         vol = exact_zonotope_vol<NT>(ZP);
 

diff --git a/R-proj/src/inner_ball.cpp b/R-proj/src/inner_ball.cpp
@@ -55,8 +55,7 @@ Rcpp::NumericVector inner_ball(Rcpp::Reference P,
     switch (type) {
         case 1: {
             // Hpolytope
-            Hpolytope HP;
-            HP.init(n, Rcpp::as<MT>(P.field("A")), Rcpp::as<VT>(P.field("b")));
+            Hpolytope HP(n, Rcpp::as<MT>(P.field("A")), Rcpp::as<VT>(P.field("b")));
             if (lp_solve) {
                 InnerBall = ComputeChebychevBall<NT, Point>(HP.get_mat(), HP.get_vec());
             } else {
@@ -67,29 +66,23 @@ Rcpp::NumericVector inner_ball(Rcpp::Reference P,
         }
         case 2: {
             // Vpolytope
-            Vpolytope VP;
-            VP.init(n, Rcpp::as<MT>(P.field("V")), VT::Ones(Rcpp::as<MT>(P.field("V")).rows()));
+            Vpolytope VP(n, Rcpp::as<MT>(P.field("V")), VT::Ones(Rcpp::as<MT>(P.field("V")).rows()));
             InnerBall = VP.ComputeInnerBall();
             if (InnerBall.second < 0.0) throw Rcpp::exception("Unable to compute a feasible point.");
             break;
         }
         case 3: {
             // Zonotope
-            zonotope ZP;
-            InnerBall = ZP.ComputeInnerBall();
-            ZP.init(n, Rcpp::as<MT>(P.field("G")), VT::Ones(Rcpp::as<MT>(P.field("G")).rows()));
+            zonotope ZP(n, Rcpp::as<MT>(P.field("G")), VT::Ones(Rcpp::as<MT>(P.field("G")).rows()));
             InnerBall = ZP.ComputeInnerBall();
             if (InnerBall.second < 0.0) throw Rcpp::exception("Unable to compute a feasible point.");
             break;
         }
         case 4: {
             // Intersection of two V-polytopes
-            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);
+            Vpolytope VP1(n, Rcpp::as<MT>(P.field("V1")), VT::Ones(Rcpp::as<MT>(P.field("V1")).rows()));
+            Vpolytope VP2(n, Rcpp::as<MT>(P.field("V2")), VT::Ones(Rcpp::as<MT>(P.field("V2")).rows()));
+            InterVP VPcVP(VP1, VP2);
             if (!VPcVP.is_feasible()) throw Rcpp::exception("Empty set!");
             InnerBall = VPcVP.ComputeInnerBall();
             if (InnerBall.second < 0.0) throw Rcpp::exception("Unable to compute a feasible point.");

diff --git a/R-proj/src/poly_gen.cpp b/R-proj/src/poly_gen.cpp
@@ -66,13 +66,13 @@ Rcpp::NumericMatrix poly_gen (int kind_gen, bool Vpoly_gen, bool Zono_gen, int d
         switch (kind_gen) {
 
             case 1:
-                return extractMatPoly(gen_cube<Vpolytope>(dim_gen, true));
+                return extractMatPoly(generate_cube<Vpolytope>(dim_gen, true));
 
             case 2:
-                return extractMatPoly(gen_cross<Vpolytope>(dim_gen, true));
+                return extractMatPoly(generate_cross<Vpolytope>(dim_gen, true));
 
             case 3:
-                return extractMatPoly(gen_simplex<Vpolytope>(dim_gen, true));
+                return extractMatPoly(generate_simplex<Vpolytope>(dim_gen, true));
 
             case 4:
                 return extractMatPoly(random_vpoly<Vpolytope, RNGType>(dim_gen, m_gen, seed_rcpp));
@@ -85,25 +85,25 @@ Rcpp::NumericMatrix poly_gen (int kind_gen, bool Vpoly_gen, bool Zono_gen, int d
         switch (kind_gen) {
 
             case 1:
-                return extractMatPoly(gen_cube<Hpolytope>(dim_gen, false));
+                return extractMatPoly(generate_cube<Hpolytope>(dim_gen, false));
 
             case 2:
-                return extractMatPoly(gen_cross<Hpolytope>(dim_gen, false));
+                return extractMatPoly(generate_cross<Hpolytope>(dim_gen, false));
 
             case 3:
-                return extractMatPoly(gen_simplex<Hpolytope>(dim_gen, false));
+                return extractMatPoly(generate_simplex<Hpolytope>(dim_gen, false));
 
             case 4:
-                return extractMatPoly(gen_prod_simplex<Hpolytope>(dim_gen));
+                return extractMatPoly(generate_prod_simplex<Hpolytope>(dim_gen));
 
             case 5:
-                return extractMatPoly(gen_skinny_cube<Hpolytope>(dim_gen));
+                return extractMatPoly(generate_skinny_cube<Hpolytope>(dim_gen));
 
             case 6:
                 return extractMatPoly(random_hpoly<Hpolytope, RNGType>(dim_gen, m_gen, seed_rcpp));
 
             case 7:
-                return extractMatPoly(gen_birk<Hpolytope>(dim_gen));
+                return extractMatPoly(generate_birkhoff<Hpolytope>(dim_gen));
 
         }
     }

diff --git a/R-proj/src/rotating.cpp b/R-proj/src/rotating.cpp
@@ -54,8 +54,7 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable<Rcpp::NumericMat
     switch (type) {
         case 1: {
             // Hpolytope
-            Hpolytope HP;
-            HP.init(n, Rcpp::as<MT>(P.field("A")), Rcpp::as<VT>(P.field("b")));
+            Hpolytope HP(n, Rcpp::as<MT>(P.field("A")), Rcpp::as<VT>(P.field("b")));
             if (T.isNotNull()) {
                 TransorfMat = Rcpp::as<MT>(T);
                 HP.linear_transformIt(TransorfMat.inverse());
@@ -67,8 +66,7 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable<Rcpp::NumericMat
         }
         case 2: {
             // Vpolytope
-            Vpolytope VP;
-            VP.init(n, Rcpp::as<MT>(P.field("V")), VT::Ones(Rcpp::as<MT>(P.field("V")).rows()));
+            Vpolytope VP(n, Rcpp::as<MT>(P.field("V")), VT::Ones(Rcpp::as<MT>(P.field("V")).rows()));
             if (T.isNotNull()) {
                 TransorfMat = Rcpp::as<MT>(T);
                 VP.linear_transformIt(TransorfMat.inverse());
@@ -80,8 +78,7 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable<Rcpp::NumericMat
         }
         case 3: {
             // Zonotope
-            zonotope ZP;
-            ZP.init(n, Rcpp::as<MT>(P.field("G")), VT::Ones(Rcpp::as<MT>(P.field("G")).rows()));
+            zonotope ZP(n, Rcpp::as<MT>(P.field("G")), VT::Ones(Rcpp::as<MT>(P.field("G")).rows()));
             if (T.isNotNull()) {
                 TransorfMat = Rcpp::as<MT>(T);
                 ZP.linear_transformIt(TransorfMat.inverse());
@@ -96,8 +93,6 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable<Rcpp::NumericMat
         }
     }
 
-
-
     TransorfMat.conservativeResize(n+1, n);
     TransorfMat.row(n) = VT::Ones(n);
     MT res(TransorfMat.rows(), Rcpp::as<MT>(Mat).rows()+n);

diff --git a/R-proj/src/rounding.cpp b/R-proj/src/rounding.cpp
@@ -90,23 +90,18 @@ Rcpp::List rounding (Rcpp::Reference P,
         rng.set_seed(seed_rcpp);
     }
 
-    Hpolytope HP;
-    Vpolytope VP;
-    zonotope ZP;
-
     std::pair <Point, NT> InnerBall;
     Rcpp::NumericMatrix Mat;
+
     std::tuple<MT, VT, NT> round_res;
     switch (type) {
         case 1: {
             // Hpolytope
-            if (Rcpp::as<MT>(P.field("Aeq")).rows() == 0) {
-                // full dimensional
-                HP.init(n, Rcpp::as<MT>(P.field("A")), Rcpp::as<VT>(P.field("b")));
-            } else {
-                // low dimensional
+
+            if (Rcpp::as<MT>(P.field("Aeq")).rows() > 0) {
                 throw Rcpp::exception("volesti supports rounding for full dimensional polytopes. Maybe call function get_full_dimensional_polytope()");
-            }
+            } 
+            Hpolytope HP(n, Rcpp::as<MT>(P.field("A")), Rcpp::as<VT>(P.field("b")));
             InnerBall = HP.ComputeInnerBall();
             if (InnerBall.second < 0.0) throw Rcpp::exception("Unable to compute a feasible point.");
             if (method_rcpp.compare(std::string("max_ellipsoid")) == 0) {
@@ -119,7 +114,7 @@ Rcpp::List rounding (Rcpp::Reference P,
         }
         case 2: {
             // Vpolytope
-            VP.init(n, Rcpp::as<MT>(P.field("V")), VT::Ones(Rcpp::as<MT>(P.field("V")).rows()));
+            Vpolytope VP(n, Rcpp::as<MT>(P.field("V")), VT::Ones(Rcpp::as<MT>(P.field("V")).rows()));
             InnerBall = VP.ComputeInnerBall();
             if (InnerBall.second < 0.0) throw Rcpp::exception("Unable to compute a feasible point.");
             round_res = apply_rounding<MT, VT, BilliardWalk>(VP, method_rcpp, walkL, InnerBall, rng);
@@ -128,7 +123,7 @@ Rcpp::List rounding (Rcpp::Reference P,
         }
         case 3: {
             // Zonotope
-            ZP.init(n, Rcpp::as<MT>(P.field("G")), VT::Ones(Rcpp::as<MT>(P.field("G")).rows()));
+            zonotope ZP(n, Rcpp::as<MT>(P.field("G")), VT::Ones(Rcpp::as<MT>(P.field("G")).rows()));
             InnerBall = ZP.ComputeInnerBall();
             if (InnerBall.second < 0.0) throw Rcpp::exception("Unable to compute a feasible point.");
             round_res = apply_rounding<MT, VT, BilliardWalk>(ZP, method_rcpp, walkL, InnerBall, rng);
@@ -143,4 +138,5 @@ Rcpp::List rounding (Rcpp::Reference P,
     return Rcpp::List::create(Rcpp::Named("Mat") = Mat, Rcpp::Named("T") = Rcpp::wrap(std::get<0>(round_res)),
                               Rcpp::Named("shift") = Rcpp::wrap(std::get<1>(round_res)),
                               Rcpp::Named("round_value") = std::get<2>(round_res));
+
 }
diff --git a/R-proj/src/sample_points.cpp b/R-proj/src/sample_points.cpp
@@ -198,7 +198,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable<Rcpp::Reference> P,
     typedef Eigen::Matrix<NT,Eigen::Dynamic,1> VT;
     typedef Eigen::Matrix<NT,Eigen::Dynamic,Eigen::Dynamic> MT;
 
-    unsigned int type = Rcpp::as<Rcpp::Reference>(P).field("type"),dim = Rcpp::as<Rcpp::Reference>(P).field("dimension"),
+    unsigned int type = Rcpp::as<Rcpp::Reference>(P).field("type"), dim = Rcpp::as<Rcpp::Reference>(P).field("dimension"),
           walkL = 1, numpoints, nburns = 0;
 
     RNGType rng(dim);
@@ -207,11 +207,6 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable<Rcpp::Reference> P,
         rng.set_seed(seed_rcpp);
     }
 
-    Hpolytope HP;
-    Vpolytope VP;
-    zonotope ZP;
-    InterVP VPcVP;
-
     NT radius = 1.0, L;
     bool set_mode = false, gaussian = false, set_starting_point = false, set_L = false;
 
@@ -357,7 +352,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable<Rcpp::Reference> P,
     switch(type) {
         case 1: {
             // Hpolytope
-            HP.init(dim, Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("A")),
+            Hpolytope HP(dim, Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("A")),
                     Rcpp::as<VT>(Rcpp::as<Rcpp::Reference>(P).field("b")));
 
             InnerBall = HP.ComputeInnerBall();
@@ -379,7 +374,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable<Rcpp::Reference> P,
         }
         case 2: {
             // Vpolytope
-            VP.init(dim, Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("V")),
+            Vpolytope VP(dim, Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("V")),
                     VT::Ones(Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("V")).rows()));
 
             InnerBall = VP.ComputeInnerBall();
@@ -400,7 +395,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable<Rcpp::Reference> P,
         }
         case 3: {
             // Zonotope
-            ZP.init(dim, Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("G")),
+            zonotope ZP(dim, Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("G")),
                     VT::Ones(Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("G")).rows()));
 
             InnerBall = ZP.ComputeInnerBall();
@@ -415,19 +410,17 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable<Rcpp::Reference> P,
                 StartingPoint = StartingPoint - mode;
                 ZP.shift(mode.getCoefficients());
             }
-            sample_from_polytope(ZP, type, rng, randPoints, walkL, numpoints, gaussian, a, L, 
+            sample_from_polytope(ZP, type, rng, randPoints, walkL, numpoints, gaussian, a, L,
                                  StartingPoint, nburns, set_L, walk);
             break;
         }
         case 4: {
             // Intersection of two V-polytopes
-            Vpolytope VP1;
-            Vpolytope VP2;
-            VP1.init(dim, Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("V1")),
+            Vpolytope VP1(dim, 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(dim, Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("V2")),
+            Vpolytope VP2(dim, 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);
+            InterVP VPcVP(VP1, VP2);
 
             if (!VPcVP.is_feasible()) throw Rcpp::exception("Empty set!");
             InnerBall = VPcVP.ComputeInnerBall();