Improve interfaces, generators and bug fixes (#58)

* use cdhr in rounding, improve t-test iterations, change diameter of H-polytopes, minor modifications in both interfaces.

* improve c++ interface, add c++ test for cooling bodies with billiard walk

* improve new c++ test

* add random generators

* update R random generators

* update Rd files

* update R volume interface

* fix generators in both c++ and R interfaces, improve exact_volume check in c++ interface

* fix bug in hpoly zonotope volume approximation

* improve c++ documentation (help command)

* fix c++ tests

* fix c++ tests

* improve cpp generator interface

* improve cpp generator interface

R-proj/R/RcppExports.R

@@ -122,8 +122,8 @@ inner_ball <- function(P) {
 #' Do not use this function.
 #'
 #' @return A numerical matrix describing the requested polytope
-poly_gen <- function(kind_gen, Vpoly_gen, dim_gen, m_gen) {
-    .Call(`_volesti_poly_gen`, kind_gen, Vpoly_gen, dim_gen, m_gen)
+poly_gen <- function(kind_gen, Vpoly_gen, Zono_gen, dim_gen, m_gen) {
+    .Call(`_volesti_poly_gen`, kind_gen, Vpoly_gen, Zono_gen, dim_gen, m_gen)
 }
 
 #'  An internal Rccp function for the random rotation of a convex polytope

R-proj/R/gen_cross.R

@@ -25,7 +25,7 @@ gen_cross <- function(dimension, repr) {
     stop('Not a known representation.')
   }
 
-  Mat = poly_gen(kind_gen, Vpoly_gen, dimension, m_gen)
+  Mat = poly_gen(kind_gen, Vpoly_gen, FALSE, dimension, m_gen)
 
   # first column is the vector b
   b = Mat[,1]

R-proj/R/gen_cube.R

@@ -25,7 +25,7 @@ gen_cube <- function(dimension, repr) {
     stop('Not a known representation.')
   }
 
-  Mat = poly_gen(kind_gen, Vpoly_gen, dimension, m_gen)
+  Mat = poly_gen(kind_gen, Vpoly_gen, FALSE, dimension, m_gen)
 
   # first column is the vector b
   b = Mat[,1]

R-proj/R/gen_prod_simplex.R

@@ -16,7 +16,7 @@ gen_prod_simplex <- function(dimension) {
   m_gen = 0
   Vpoly_gen = FALSE
 
-  Mat = poly_gen(kind_gen, Vpoly_gen, dimension, m_gen)
+  Mat = poly_gen(kind_gen, Vpoly_gen, FALSE, dimension, m_gen)
 
   # first column is the vector b
   b = Mat[,1]

R-proj/R/gen_rand_hpoly.R

@@ -15,7 +15,7 @@ gen_rand_hpoly <- function(dimension, m) {
   kind_gen = 6
   Vpoly_gen = FALSE
 
-  Mat = poly_gen(kind_gen, Vpoly_gen, dimension, m)
+  Mat = poly_gen(kind_gen, Vpoly_gen, FALSE, dimension, m)
 
   # first column is the vector b
   b = Mat[,1]

R-proj/R/gen_rand_vpoly.R

@@ -3,19 +3,27 @@
 #' This function can be used to generate a \eqn{d}-dimensional polytope in V-representation with \eqn{m} vertices. We pick \eqn{m} random points from the boundary of the \eqn{d}-dimensional unit hypersphere as vertices.
 #' 
 #' @param dimension The dimension of the convex polytope.
-#' @param m The number of the vertices.
+#' @param n_vertices The number of the vertices.
+#' @param generator The body that the generator samples uniformly the vertices from: (a) 'cube' or (b) 'sphere'.
 #' 
 #' @return A polytope class representing a V-polytope.
 #' @examples 
 #' # generate a 10-dimensional polytope defined as the convex hull of 25 random vertices
 #' P = gen_rand_vpoly(10, 25)
 #' @export
-gen_rand_vpoly <- function(dimension, m) {
+gen_rand_vpoly <- function(dimension, n_vertices, generator = NULL) {
 
   kind_gen = 4
-  Vpoly_gen = TRUE
 
-  Mat = poly_gen(kind_gen, Vpoly_gen, dimension, m)
+  if(!is.null(generator)){
+    if (generator == 'cube'){
+      kind_gen = 5
+    } else if (generator != 'sphere') {
+      stop("Wrong generator!")
+    }
+  }
+
+  Mat = poly_gen(kind_gen, TRUE, FALSE, dimension, n_vertices)
 
   # first column is the vector b
   b = Mat[,1]

R-proj/R/gen_rand_zonotope.R

@@ -3,20 +3,30 @@
 #' This function can be used to generate a random \eqn{d}-dimensional zonotope defined by the Minkowski sum of \eqn{m} \eqn{d}-dimensional segments. We consider \eqn{m} random directions in \eqn{R^d} and for each direction we pick a random length in \eqn{[(,\sqrt{d}]} in order to define \eqn{m} segments.
 #' 
 #' @param dimension The dimension of the zonotope.
-#' @param NumGen The number of segments that generate the zonotope.
+#' @param n_segments The number of segments that generate the zonotope.
+#' @param generator The distribution to pick the length of each segment from \eqn{[0,100]}: (a) 'uniform', (b) 'gaussian' or (c) 'exponential'.
 #' 
 #' @return A polytope class representing a zonotope.
 #'
 #' @examples 
 #' # generate a 10-dimensional zonotope defined by the Minkowski sum of 20 segments
 #' P = gen_rand_zonotope(10, 20)
 #' @export
-gen_rand_zonotope <- function(dimension, NumGen) {
+gen_rand_zonotope <- function(dimension, n_segments, generator = NULL) {
 
-  kind_gen = 0
-  Vpoly_gen = FALSE
+  kind_gen = 1
 
-  Mat = poly_gen(kind_gen, Vpoly_gen, dimension, NumGen)
+  if (!is.null(generator)) {
+    if (generator == 'gaussian') {
+      kind_gen = 2
+    } else if (generator == 'exponential') {
+      kind_gen = 3
+    } else if (generator != 'uniform'){
+      stop("Wrong generator!")
+    }
+  }
+
+  Mat = poly_gen(kind_gen, FALSE, TRUE, dimension, n_segments)
 
   # first column is the vector b
   b = Mat[,1]

R-proj/R/gen_simplex.R

@@ -25,7 +25,7 @@ gen_simplex <- function(dimension, repr) {
     stop('Not a known representation.')
   }
 
-  Mat = poly_gen(kind_gen, Vpoly_gen, dimension, m_gen)
+  Mat = poly_gen(kind_gen, Vpoly_gen, FALSE, dimension, m_gen)
 
   # first column is the vector b
   b = Mat[,1]

R-proj/R/gen_skinny_cube.R

@@ -16,7 +16,7 @@ gen_skinny_cube <- function(dimension) {
   m_gen = 0
   Vpoly_gen = FALSE
 
-  Mat = poly_gen(kind_gen, Vpoly_gen, dimension, m_gen)
+  Mat = poly_gen(kind_gen, Vpoly_gen, FALSE, dimension, m_gen)
 
   # first column is the vector b
   b = Mat[,1]

R-proj/R/round_polytope.R

@@ -45,10 +45,8 @@ round_polytope <- function(P, random_walk = NULL, walk_length = NULL, parameters
     PP = list("P" = Vpolytope$new(A), "round_value" = ret_list$round_value)
   }else if (type == 3) {
     PP = list("P" = Zonotope$new(A), "round_value" = ret_list$round_value)
-  } else if(type == 1) {
-    PP = list("P" = Hpolytope$new(A,b), "round_value" = ret_list$round_value)
   } else {
-    PP = list("P" = VPolyintersectVPoly$new(V1 = t(Mat[,dim(P$V1)[1]]), V2 = t(Mat[,dim(P$V2)[1]])), "round_value" = ret_list$round_value)
+    PP = list("P" = Hpolytope$new(A,b), "round_value" = ret_list$round_value)
   }
   return(PP)
 }

R-proj/src/poly_gen.cpp

@@ -18,7 +18,10 @@
 #include "hpolytope.h"
 #include "vpolytope.h"
 #include "zpolytope.h"
-#include "polytope_generators.h"
+#include "known_polytope_generators.h"
+#include "h_polytopes_gen.h"
+#include "v_polytopes_gen.h"
+#include "z_polytopes_gen.h"
 #include "extractMatPoly.h"
 
 //' An internal Rccp function as a polytope generator
@@ -33,7 +36,7 @@
 //'
 //' @return A numerical matrix describing the requested polytope
 // [[Rcpp::export]]
-Rcpp::NumericMatrix poly_gen (int kind_gen, bool Vpoly_gen, int dim_gen, int m_gen) {
+Rcpp::NumericMatrix poly_gen (int kind_gen, bool Vpoly_gen, bool Zono_gen, int dim_gen, int m_gen) {
 
     typedef double NT;
     typedef Cartesian <NT> Kernel;
@@ -43,9 +46,17 @@ Rcpp::NumericMatrix poly_gen (int kind_gen, bool Vpoly_gen, int dim_gen, int m_g
     typedef VPolytope <Point, RNGType> Vpolytope;
     typedef Zonotope <Point> zonotope;
 
-    if (kind_gen == 0) {
+    if (Zono_gen) {
+        switch (kind_gen) {
 
-        return extractMatPoly(gen_zonotope<zonotope, RNGType>(dim_gen, m_gen));
+            case 1:
+                return extractMatPoly(gen_zonotope_uniform<zonotope, RNGType>(dim_gen, m_gen));
+            case 2:
+                return extractMatPoly(gen_zonotope_gaussian<zonotope, RNGType>(dim_gen, m_gen));
+            case 3:
+                return extractMatPoly(gen_zonotope_exponential<zonotope, RNGType>(dim_gen, m_gen));
+
+        }
 
     } else if (Vpoly_gen) {
         switch (kind_gen) {
@@ -62,6 +73,9 @@ Rcpp::NumericMatrix poly_gen (int kind_gen, bool Vpoly_gen, int dim_gen, int m_g
             case 4:
                 return extractMatPoly(random_vpoly<Vpolytope, RNGType>(dim_gen, m_gen));
 
+            case 5:
+                return extractMatPoly(random_vpoly_incube<Vpolytope, RNGType>(dim_gen, m_gen));
+
         }
     } else {
         switch (kind_gen) {
@@ -87,6 +101,6 @@ Rcpp::NumericMatrix poly_gen (int kind_gen, bool Vpoly_gen, int dim_gen, int m_g
         }
     }
 
-    return Rcpp::NumericMatrix(0, 0);
+    throw Rcpp::exception("Wrong inputs!");
 
 }

R-proj/src/volume.cpp

@@ -200,35 +200,6 @@ double volume (Rcpp::Reference P,  Rcpp::Nullable<unsigned int> walk_length = R_
     NT C = 2.0, ratio = 1.0-1.0/(NT(n)), frac = 0.1, e, delta = -1.0, lb = 0.1, ub = 0.15, p = 0.75, rmax = 0.0,
             alpha = 0.2, diam = -1.0;
 
-    if (!random_walk.isNotNull()) {
-        if ( type == 1 ){
-            cdhr = true;
-        } else {
-            billiard = true;
-            win_len = 150;
-            NN = 125;
-        }
-    }else if (Rcpp::as<std::string>(random_walk).compare(std::string("CDHR")) == 0) {
-        cdhr = true;
-        win_len = 3*n*n+400;
-    } else if (Rcpp::as<std::string>(random_walk).compare(std::string("RDHR")) == 0) {
-        rdhr = true;
-    } else if (Rcpp::as<std::string>(random_walk).compare(std::string("BaW"))==0) {
-        ball_walk = true;
-    } else if (Rcpp::as<std::string>(random_walk).compare(std::string("BiW"))==0) {
-        billiard = true;
-        win_len = 150;
-        NN = 125;
-    }else {
-        throw Rcpp::exception("Unknown walk type!");
-    }
-
-    if (!rounding.isNotNull() && type == 2){
-        round = true;
-    } else {
-        round = (!rounding.isNotNull()) ? false : Rcpp::as<bool>(rounding);
-    }
-
     if(!algo.isNotNull()){
 
         if (type == 2 || type == 3) {
@@ -257,14 +228,61 @@ double volume (Rcpp::Reference P,  Rcpp::Nullable<unsigned int> walk_length = R_
         CB = true;
         e = (!error.isNotNull()) ? 0.1 : Rcpp::as<NT>(error);
         walkL = (!walk_length.isNotNull()) ? 1 : Rcpp::as<int>(walk_length);
-        win_len = (cdhr) ? 3*n*n+400 : 2*n*n+250;
-        if (billiard){
-            win_len = 150;
+
+    } else {
+        throw Rcpp::exception("Unknown method!");
+    }
+
+    if (!random_walk.isNotNull()) {
+        if ( type == 1 ){
+            cdhr = true;
+            if (CB) win_len = 3*n*n+400;
+        } else {
+            if (CB) {
+                billiard = true;
+                win_len = 170;
+                NN = 125;
+            } else {
+                rdhr = true;
+            }
+        }
+    }else if (Rcpp::as<std::string>(random_walk).compare(std::string("CDHR")) == 0) {
+        cdhr = true;
+        if (CB) win_len = 3*n*n+400;
+    } else if (Rcpp::as<std::string>(random_walk).compare(std::string("RDHR")) == 0) {
+        rdhr = true;
+    } else if (Rcpp::as<std::string>(random_walk).compare(std::string("BaW"))==0) {
+        ball_walk = true;
+    } else if (Rcpp::as<std::string>(random_walk).compare(std::string("BiW"))==0) {
+        if (CG) {
+            if (type !=1){
+                Rcpp::Rcout << "Billiard walk is not supported for CG algorithm. RDHR is used."<<std::endl;
+                rdhr = true;
+            } else {
+                Rcpp::Rcout << "Billiard walk is not supported for CG algorithm. CDHR is used."<<std::endl;
+                cdhr = true;
+            }
+        } else if (!CB) {
+            if (type !=1){
+                Rcpp::Rcout << "Billiard walk is not supported for SOB algorithm. RDHR is used."<<std::endl;
+                rdhr = true;
+            } else {
+                Rcpp::Rcout << "Billiard walk is not supported for SOB algorithm. CDHR is used."<<std::endl;
+                cdhr = true;
+            }
+        } else {
+            billiard = true;
+            win_len = 170;
             NN = 125;
         }
+    }else {
+        throw Rcpp::exception("Unknown walk type!");
+    }
 
+    if (!rounding.isNotNull() && type == 2){
+        round = true;
     } else {
-        throw Rcpp::exception("Unknown method!");
+        round = (!rounding.isNotNull()) ? false : Rcpp::as<bool>(rounding);
     }
 
     if(parameters.isNotNull()) {
@@ -319,8 +337,6 @@ double volume (Rcpp::Reference P,  Rcpp::Nullable<unsigned int> walk_length = R_
         }
     }
 
-    //std::cout<<"Wlen = "<<win_len<<", NN = "<<NN<<", nu = "<<nu<<", lb = "<<lb<<", ub = "<<ub<<", alpha = "<<alpha<<", p = "<<p<<", W = "<<walkL<<", billiard = "<<billiard<<std::endl;
-
     switch(type) {
         case 1: {
             // Hpolytope

include/annealing/ball_annealing.h

@@ -16,14 +16,14 @@ bool check_convergence(ConvexBody &P, PointList &randPoints, const NT &lb, const
 
     std::vector<NT> ratios;
     std::pair<NT,NT> mv;
-    int m = randPoints.size()/nu;
+    int m = randPoints.size()/nu, i = 1;
     NT T, rs, alpha_check = 0.01;
     size_t countsIn = 0;
 
-    for(typename PointList::iterator pit=randPoints.begin(); pit!=randPoints.end(); ++pit){
+    for(typename std::list<Point>::iterator pit=randPoints.begin(); pit!=randPoints.end(); ++pit, i++){
 
         if (P.is_in(*pit)==-1) countsIn++;
-        if ((std::distance(randPoints.begin(), pit)+1) % m == 0) {
+        if (i % m == 0) {
             ratios.push_back(NT(countsIn)/m);
             countsIn = 0;
             if (ratios.size()>1 && precheck) {
@@ -47,7 +47,8 @@ bool check_convergence(ConvexBody &P, PointList &randPoints, const NT &lb, const
     ratio = mv.first;
     rs = std::sqrt(mv.second);
     boost::math::students_t dist(nu - 1);
-    T = rs*(boost::math::quantile(boost::math::complement(dist, alpha)) / std::sqrt(NT(nu)));
+    T = rs*(boost::math::quantile(boost::math::complement(dist, alpha))
+            / std::sqrt(NT(nu)));
     if (ratio > lb + T) {
         if (lastball) return true;
         if ((precheck && ratio < ub - T) || (!precheck && ratio < ub + T)) return true;

include/annealing/hpoly_annealing.h

@@ -147,7 +147,8 @@ bool get_sequence_of_zonopolys(Zonotope &Z, const HPolytope &HP, std::vector<HPo
     typedef typename Zonotope::MT MT;
 
     int n = var.n;
-    MT G = Z.get_mat().transpose(), AG = HP.get_mat()*G;
+    MT G = Z.get_mat().transpose();
+    MT AG = HP.get_mat()*G;
     NT ratio;
     std::list<Point> randPoints;
     Point q(n);