diff --git a/R-proj/R/RcppExports.R b/R-proj/R/RcppExports.R index 58db215af..7ebdb5e6f 100644 --- a/R-proj/R/RcppExports.R +++ b/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 diff --git a/R-proj/R/gen_cross.R b/R-proj/R/gen_cross.R index 026735ae0..6101f5ee1 100644 --- a/R-proj/R/gen_cross.R +++ b/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] diff --git a/R-proj/R/gen_cube.R b/R-proj/R/gen_cube.R index cce62489c..5ecf1bcc0 100644 --- a/R-proj/R/gen_cube.R +++ b/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] diff --git a/R-proj/R/gen_prod_simplex.R b/R-proj/R/gen_prod_simplex.R index 95e506948..2eda50bd9 100644 --- a/R-proj/R/gen_prod_simplex.R +++ b/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] diff --git a/R-proj/R/gen_rand_hpoly.R b/R-proj/R/gen_rand_hpoly.R index 0093e09e5..9a5059119 100644 --- a/R-proj/R/gen_rand_hpoly.R +++ b/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] diff --git a/R-proj/R/gen_rand_vpoly.R b/R-proj/R/gen_rand_vpoly.R index 8b8d7e5da..47001dbb5 100644 --- a/R-proj/R/gen_rand_vpoly.R +++ b/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] diff --git a/R-proj/R/gen_rand_zonotope.R b/R-proj/R/gen_rand_zonotope.R index e13a079be..f35798be9 100644 --- a/R-proj/R/gen_rand_zonotope.R +++ b/R-proj/R/gen_rand_zonotope.R @@ -3,7 +3,8 @@ #' 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. #' @@ -11,12 +12,21 @@ #' # 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] diff --git a/R-proj/R/gen_simplex.R b/R-proj/R/gen_simplex.R index b1f04d917..6d930cdd4 100644 --- a/R-proj/R/gen_simplex.R +++ b/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] diff --git a/R-proj/R/gen_skinny_cube.R b/R-proj/R/gen_skinny_cube.R index 2c2dbc67f..bcf3dda40 100644 --- a/R-proj/R/gen_skinny_cube.R +++ b/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] diff --git a/R-proj/R/round_polytope.R b/R-proj/R/round_polytope.R index f980b8639..16b8367ae 100644 --- a/R-proj/R/round_polytope.R +++ b/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) } diff --git a/R-proj/man/gen_rand_vpoly.Rd b/R-proj/man/gen_rand_vpoly.Rd index 14d904e64..a994ce0ef 100644 --- a/R-proj/man/gen_rand_vpoly.Rd +++ b/R-proj/man/gen_rand_vpoly.Rd @@ -4,12 +4,14 @@ \alias{gen_rand_vpoly} \title{Generator function for random V-polytopes} \usage{ -gen_rand_vpoly(dimension, m) +gen_rand_vpoly(dimension, n_vertices, generator = NULL) } \arguments{ \item{dimension}{The dimension of the convex polytope.} -\item{m}{The number of the vertices.} +\item{n_vertices}{The number of the vertices.} + +\item{generator}{The body that the generator samples uniformly the vertices from: (a) 'cube' or (b) 'sphere'.} } \value{ A polytope class representing a V-polytope. diff --git a/R-proj/man/gen_rand_zonotope.Rd b/R-proj/man/gen_rand_zonotope.Rd index 67b3f7f45..8278490d4 100644 --- a/R-proj/man/gen_rand_zonotope.Rd +++ b/R-proj/man/gen_rand_zonotope.Rd @@ -4,12 +4,14 @@ \alias{gen_rand_zonotope} \title{Generator function for zonotopes} \usage{ -gen_rand_zonotope(dimension, NumGen) +gen_rand_zonotope(dimension, n_segments, generator = NULL) } \arguments{ \item{dimension}{The dimension of the zonotope.} -\item{NumGen}{The number of segments that generate the zonotope.} +\item{n_segments}{The number of segments that generate the zonotope.} + +\item{generator}{The distribution to pick the length of each segment from \eqn{[0,100]}: (a) 'uniform', (b) 'gaussian' or (c) 'exponential'.} } \value{ A polytope class representing a zonotope. diff --git a/R-proj/man/poly_gen.Rd b/R-proj/man/poly_gen.Rd index d199997d5..2e7ce5e02 100644 --- a/R-proj/man/poly_gen.Rd +++ b/R-proj/man/poly_gen.Rd @@ -4,7 +4,7 @@ \alias{poly_gen} \title{An internal Rccp function as a polytope generator} \usage{ -poly_gen(kind_gen, Vpoly_gen, dim_gen, m_gen) +poly_gen(kind_gen, Vpoly_gen, Zono_gen, dim_gen, m_gen) } \arguments{ \item{kind_gen}{An integer to declare the type of the polytope.} diff --git a/R-proj/src/poly_gen.cpp b/R-proj/src/poly_gen.cpp index f9a9e618f..03cea06fb 100644 --- a/R-proj/src/poly_gen.cpp +++ b/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 Kernel; @@ -43,9 +46,17 @@ Rcpp::NumericMatrix poly_gen (int kind_gen, bool Vpoly_gen, int dim_gen, int m_g typedef VPolytope Vpolytope; typedef Zonotope zonotope; - if (kind_gen == 0) { + if (Zono_gen) { + switch (kind_gen) { - return extractMatPoly(gen_zonotope(dim_gen, m_gen)); + case 1: + return extractMatPoly(gen_zonotope_uniform(dim_gen, m_gen)); + case 2: + return extractMatPoly(gen_zonotope_gaussian(dim_gen, m_gen)); + case 3: + return extractMatPoly(gen_zonotope_exponential(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(dim_gen, m_gen)); + case 5: + return extractMatPoly(random_vpoly_incube(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!"); } diff --git a/R-proj/src/volume.cpp b/R-proj/src/volume.cpp index 0ef136cc5..52bb99cdc 100644 --- a/R-proj/src/volume.cpp +++ b/R-proj/src/volume.cpp @@ -200,35 +200,6 @@ double volume (Rcpp::Reference P, Rcpp::Nullable 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(random_walk).compare(std::string("CDHR")) == 0) { - cdhr = true; - win_len = 3*n*n+400; - } else if (Rcpp::as(random_walk).compare(std::string("RDHR")) == 0) { - rdhr = true; - } else if (Rcpp::as(random_walk).compare(std::string("BaW"))==0) { - ball_walk = true; - } else if (Rcpp::as(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(rounding); - } - if(!algo.isNotNull()){ if (type == 2 || type == 3) { @@ -257,14 +228,61 @@ double volume (Rcpp::Reference P, Rcpp::Nullable walk_length = R_ CB = true; e = (!error.isNotNull()) ? 0.1 : Rcpp::as(error); walkL = (!walk_length.isNotNull()) ? 1 : Rcpp::as(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(random_walk).compare(std::string("CDHR")) == 0) { + cdhr = true; + if (CB) win_len = 3*n*n+400; + } else if (Rcpp::as(random_walk).compare(std::string("RDHR")) == 0) { + rdhr = true; + } else if (Rcpp::as(random_walk).compare(std::string("BaW"))==0) { + ball_walk = true; + } else if (Rcpp::as(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."<(rounding); } if(parameters.isNotNull()) { @@ -319,8 +337,6 @@ double volume (Rcpp::Reference P, Rcpp::Nullable walk_length = R_ } } - //std::cout<<"Wlen = "< +Polytope random_hpoly(unsigned int dim, unsigned int m) { + + typedef typename Polytope::MT MT; + typedef typename Polytope::VT VT; + typedef typename Polytope::NT NT; + typedef typename Polytope::PolytopePoint Point; + + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + RNGType rng(seed); + boost::random::uniform_real_distribution<> urdist1(-10, 10); + Point p(dim); + typename std::vector::iterator pit; + MT A(m, dim); + VT b(m); + unsigned int j; + + for(unsigned int i=0; i(dim); + pit = p.iter_begin(); + j = 0; + for ( ; pit!=p.iter_end(); ++pit, ++j) { + A(i,j) = *pit; + } + b(i) = 10.0; + + } + Polytope HP; + HP.init(dim, A, b); + + return HP; +} + +#endif diff --git a/include/generators/polytope_generators.h b/include/generators/known_polytope_generators.h similarity index 76% rename from include/generators/polytope_generators.h rename to include/generators/known_polytope_generators.h index 9653f95c0..df52215ba 100644 --- a/include/generators/polytope_generators.h +++ b/include/generators/known_polytope_generators.h @@ -5,8 +5,8 @@ // Licensed under GNU LGPL.3, see LICENCE file -#ifndef POLYTOPE_GENERATORS_H -#define POLYTOPE_GENERATORS_H +#ifndef KNOWN_POLYTOPE_GENERATORS_H +#define KNOWN_POLYTOPE_GENERATORS_H #include #include "samplers.h" @@ -306,103 +306,6 @@ Polytope gen_skinny_cube(const unsigned int &dim, bool Vpoly = false) { return P; } - -template -Polytope gen_zonotope(const unsigned int &dim, const unsigned int &m) { - - typedef typename Polytope::MT MT; - typedef typename Polytope::VT VT; - typedef typename Polytope::NT NT; - - unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); - RNGType rng(seed); - boost::normal_distribution<> rdist(0, 1); - - MT A; - VT b; - A.resize(m, dim); - b.resize(m); - Polytope P; - - for (unsigned int i = 0; i < m; ++i) { - b(i) = 1.0; - for (unsigned int j = 0; j < dim; ++j) { - A(i,j) = rdist(rng); - } - } - - P.init(dim, A, b); - return P; -} - -template -Polytope random_vpoly(const unsigned int &dim, const unsigned int &k) { - - typedef typename Polytope::MT MT; - typedef typename Polytope::VT VT; - typedef typename Polytope::NT NT; - typedef typename Polytope::PolytopePoint Point; - - unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); - RNGType rng(seed); - - Point p; - typename std::vector::iterator pit; - MT V(k, dim); - unsigned int j; - - for (unsigned int i = 0; i < k; ++i) { - p = get_direction(dim); - pit = p.iter_begin(); - j = 0; - for ( ; pit!=p.iter_end(); ++pit, ++j) { - V(i,j) = *pit; - } - } - - Polytope VP; - VT b = VT::Ones(k); - VP.init(dim, V, b); - - return VP; - -} - -template -Polytope random_hpoly(const unsigned int &dim, const unsigned int &m) { - - typedef typename Polytope::MT MT; - typedef typename Polytope::VT VT; - typedef typename Polytope::NT NT; - typedef typename Polytope::PolytopePoint Point; - - unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); - RNGType rng(seed); - boost::random::uniform_real_distribution<> urdist1(-10, 10); - Point p(dim); - typename std::vector::iterator pit; - MT A(m, dim); - VT b(m); - unsigned int j; - - for(unsigned int i=0; i(dim); - pit = p.iter_begin(); - j = 0; - for ( ; pit!=p.iter_end(); ++pit, ++j) { - A(i,j) = *pit; - } - b(i) = 10.0; - - } - Polytope HP; - HP.init(dim, A, b); - - return HP; -} - - - /* * ToDo: brkhoff polytope generator template diff --git a/include/generators/v_polytopes_gen.h b/include/generators/v_polytopes_gen.h new file mode 100644 index 000000000..73e148bbf --- /dev/null +++ b/include/generators/v_polytopes_gen.h @@ -0,0 +1,149 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 20012-2018 Vissarion Fisikopoulos +// Copyright (c) 2018 Apostolos Chalkis + +// Licensed under GNU LGPL.3, see LICENCE file + +#ifndef V_POLYTOPES_GEN_H +#define V_POLYTOPES_GEN_H + +#include +#include "samplers.h" + + +template +void removeRow(MT &matrix, unsigned int rowToRemove) +{ + unsigned int numRows = matrix.rows()-1; + unsigned int numCols = matrix.cols(); + + if( rowToRemove < numRows ) + matrix.block(rowToRemove,0,numRows-rowToRemove,numCols) = matrix.bottomRows(numRows-rowToRemove); + + matrix.conservativeResize(numRows,numCols); +} + +template +Polytope random_vpoly(unsigned int dim, unsigned int k) { + + typedef typename Polytope::MT MT; + typedef typename Polytope::VT VT; + typedef typename Polytope::NT NT; + typedef typename Polytope::PolytopePoint Point; + + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + RNGType rng(seed); + + Point p; + typename std::vector::iterator pit; + MT V(k, dim); + unsigned int j; + + for (unsigned int i = 0; i < k; ++i) { + p = get_direction(dim); + pit = p.iter_begin(); + j = 0; + for ( ; pit!=p.iter_end(); ++pit, ++j) { + V(i,j) = *pit; + } + } + + Polytope VP; + VT b = VT::Ones(k); + VP.init(dim, V, b); + + return VP; + +} + + +template +Polytope random_vpoly_incube(unsigned int d, unsigned int k) { + + typedef typename Polytope::MT MT; + typedef typename Polytope::VT VT; + typedef typename Polytope::NT NT; + typedef typename Polytope::PolytopePoint Point; + + REAL *conv_mem; + int *colno_mem; + + conv_mem = (REAL *) malloc(k * sizeof(*conv_mem)); + colno_mem = (int *) malloc(k * sizeof(*colno_mem)); + + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + RNGType rng(seed); + boost::random::uniform_real_distribution<> urdist1(-1, 1); + + Point p(d); + typename std::vector::iterator pit; + MT V(k, d); + unsigned int j, count_row,it=0; + std::vector indices; + Polytope VP; + VT b = VT::Ones(k); + + for (unsigned int i = 0; i < k; ++i) { + for (int j = 0; j < d; ++j) { + V(i, j) = urdist1(rng); + } + } + if(k==d+1){ + VP.init(d, V, b); + return VP; + } + + MT V2(k,d); + V2 = V; + indices.clear(); + while(it<20) { + V.resize(V2.rows(), d); + V = V2; + for (int i = 0; i < indices.size(); ++i) { + V.conservativeResize(V.rows()+1, d); + for (int j = 0; j < d; ++j) { + V(V.rows()-1, j) = urdist1(rng); + } + } + indices.clear(); + V2.resize(k, d); + V2 = V; + + for (int i = 0; i < k; ++i) { + for (int j = 0; j < d; ++j) { + p.set_coord(j, V(i, j)); + } + removeRow(V2, i); + if (memLP_Vpoly(V2, p, conv_mem, colno_mem)){ + indices.push_back(i); + } + V2.resize(k, d); + V2 = V; + } + if (indices.size()==0) { + VP.init(d, V, b); + return VP; + } + V2.resize(k - indices.size(), d); + count_row =0; + for (int i = 0; i < k; ++i) { + if(std::find(indices.begin(), indices.end(), i) != indices.end()) { + continue; + } else { + for (int j = 0; j < d; ++j) V2(count_row, j) = V(i,j); + count_row++; + } + } + it++; + } + + VP.init(d, V2, VT::Ones(V2.rows())); + free(colno_mem); + free(conv_mem); + + return VP; + +} + +#endif diff --git a/include/generators/z_polytopes_gen.h b/include/generators/z_polytopes_gen.h new file mode 100644 index 000000000..3aa818df0 --- /dev/null +++ b/include/generators/z_polytopes_gen.h @@ -0,0 +1,123 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 20012-2018 Vissarion Fisikopoulos +// Copyright (c) 2018 Apostolos Chalkis + +// Licensed under GNU LGPL.3, see LICENCE file + +#ifndef Z_POLYTOPES_GEN_H +#define Z_POLYTOPES_GEN_H + +#include +#include "samplers.h" + +template +Polytope gen_zonotope_gaussian(int dim, int m) { + + typedef typename Polytope::MT MT; + typedef typename Polytope::VT VT; + typedef typename Polytope::NT NT; + + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + RNGType rng(seed); + boost::normal_distribution<> rdist(0, 1); + boost::normal_distribution<> rdist2(50, 33.3); + + MT A; + VT b; + A.resize(m, dim); + b.resize(m); + Polytope P; + NT rand_gaus; + + for (unsigned int i = 0; i < m; ++i) { + b(i) = 1.0; + for (unsigned int j = 0; j < dim; ++j) { + A(i,j) = rdist(rng); + } + A.row(i)=A.row(i)/A.row(i).norm(); + while(true){ + rand_gaus = rdist2(rng); + if (rand_gaus > 0.0 && rand_gaus<100.0){ + A.row(i) = A.row(i) * rand_gaus; + break; + } + } + } + + P.init(dim, A, b); + return P; +} + + +template +Polytope gen_zonotope_uniform(int dim, int m) { + + typedef typename Polytope::MT MT; + typedef typename Polytope::VT VT; + typedef typename Polytope::NT NT; + + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + RNGType rng(seed); + boost::normal_distribution<> rdist(0, 1); + boost::random::uniform_real_distribution<> urdist1(0, 100); + + MT A; + VT b; + A.resize(m, dim); + b.resize(m); + Polytope P; + + for (unsigned int i = 0; i < m; ++i) { + b(i) = 1.0; + for (unsigned int j = 0; j < dim; ++j) { + A(i,j) = rdist(rng); + } + A.row(i)=A.row(i)/A.row(i).norm(); + A.row(i) = A.row(i) * urdist1(rng); + } + + P.init(dim, A, b); + return P; +} + + +template +Polytope gen_zonotope_exponential(int dim, int m) { + + typedef typename Polytope::MT MT; + typedef typename Polytope::VT VT; + typedef typename Polytope::NT NT; + + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + RNGType rng(seed); + boost::normal_distribution<> rdist(0, 1); + boost::normal_distribution<> expdist(1.0/30.0); + + MT A; + VT b; + A.resize(m, dim); + b.resize(m); + Polytope P; + NT rand_exp; + + for (unsigned int i = 0; i < m; ++i) { + b(i) = 1.0; + for (unsigned int j = 0; j < dim; ++j) { + A(i,j) = rdist(rng); + } + A.row(i)=A.row(i)/A.row(i).norm(); + while(true){ + rand_exp = expdist(rng); + if (rand_exp > 0.0 && rand_exp<100.0){ + A.row(i) = A.row(i) * rand_exp; + break; + } + } + } + + P.init(dim, A, b); + return P; +} + +#endif diff --git a/include/rounding.h b/include/rounding.h index d79fd82a4..65d6a6389 100644 --- a/include/rounding.h +++ b/include/rounding.h @@ -34,11 +34,14 @@ std::pair rounding_min_ellipsoid(Polytope &P , const std::pair::iterator pit=randPoints.begin(); pit!=randPoints.end(); ++pit){ current_dist=(*pit-c).squared_length(); diff --git a/include/volume/cooling_balls.h b/include/volume/cooling_balls.h index cd1c0b7e3..417d0d98c 100644 --- a/include/volume/cooling_balls.h +++ b/include/volume/cooling_balls.h @@ -37,6 +37,7 @@ NT vol_cooling_balls(Polytope &P, UParameters &var, AParameters &var_ban, std::p std::vector BallSet; std::vector ratios; Point c = InnerBall.first; + P.normalize(); if (round) { #ifdef VOLESTI_DEBUG @@ -52,7 +53,10 @@ NT vol_cooling_balls(Polytope &P, UParameters &var, AParameters &var_ban, std::p std::pair res = P.ComputeInnerBall(); c = res.first; radius = res.second; - P.comp_diam(var.diameter, radius); + P.normalize(); + if (var.bill_walk) { + P.comp_diam(var.diameter, radius); + } if (var.ball_walk){ var.delta = 4.0 * radius / NT(n); } @@ -63,7 +67,6 @@ NT vol_cooling_balls(Polytope &P, UParameters &var, AParameters &var_ban, std::p // Move the chebychev center to the origin and apply the same shifting to the polytope VT c_e = Eigen::Map(&c.get_coeffs()[0], c.dimension()); P.shift(c_e); - P.normalize(); if ( !get_sequence_of_polyballs(P, BallSet, ratios, N * nu, nu, lb, ub, radius, alpha, var, rmax) ){ return -1.0; diff --git a/include/volume/cooling_hpoly.h b/include/volume/cooling_hpoly.h index 20372f832..57777bd5a 100644 --- a/include/volume/cooling_hpoly.h +++ b/include/volume/cooling_hpoly.h @@ -54,7 +54,7 @@ NT vol_cooling_hpoly (Zonotope &ZP, UParameters &var, AParameters &var_ban, GPar VT Zs_max(2*m); UParameters var3 = var; var3.cdhr_walk = true; - var3.ball_walk = var3.rdhr_walk = false; + var3.ball_walk = var3.rdhr_walk = var3.bill_walk = false; if ( !get_first_poly(ZP, HP, Zs_max, lb, ub, ratio, var3) ) { return -1.0; } @@ -88,7 +88,7 @@ NT vol_cooling_hpoly (Zonotope &ZP, UParameters &var, AParameters &var_ban, GPar if (!window2) { UParameters var2 = var; var2.cdhr_walk = true; - var2.ball_walk = var2.rdhr_walk = false; + var2.ball_walk = var2.rdhr_walk = var2.bill_walk = false; var2.walk_steps = 10+2*n; vol *= esti_ratio_interval(HP, ZP, ratio, er0, win_len, N*nu, prob, var2); } else { diff --git a/include/volume/volume.h b/include/volume/volume.h index 906f711e2..f5aeac3a3 100644 --- a/include/volume/volume.h +++ b/include/volume/volume.h @@ -96,7 +96,7 @@ NT volume(Polytope &P, Point p = get_point_on_Dsphere(n, radius); std::list randPoints; //ds for storing rand points //use a large walk length e.g. 1000 - + rand_point_generator(P, p, 1, 50*n, randPoints, var); double tstart2 = (double)clock()/(double)CLOCKS_PER_SEC; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 8fcb7a815..50c8e0eca 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -73,12 +73,13 @@ else () add_executable (VpolyCV_test VpolyCV_test.cpp $) add_executable (VpolyVol_test VpolyVol_test.cpp $) add_executable (ZonotopeVol_test ZonotopeVol_test.cpp $) + add_executable (cool_bodies_bill_test cooling_bodies_bill_test.cpp $) #add_executable (ZonotopeVolCV_test ZonotopeVolCV_test.cpp $) add_test(NAME volume_cube COMMAND volume_test -tc=cube) add_test(NAME volume_cross COMMAND volume_test -tc=cross) add_test(NAME volume_birkhoff COMMAND volume_test -tc=birk) - add_test(NAME volume_prod_simplex COMMAND volume_test -tc=prod_simplex) + #add_test(NAME volume_prod_simplex COMMAND volume_test -tc=prod_simplex) add_test(NAME volume_simplex COMMAND volume_test -tc=simplex) add_test(NAME volume_skinny_cube COMMAND volume_test -tc=skinny_cube) @@ -108,6 +109,13 @@ else () add_test(NAME cheb_simplex COMMAND cheb_test -tc=cheb_simplex) add_test(NAME cheb_skinny_cube COMMAND cheb_test -tc=cheb_skinny_cube) + add_test(NAME cool_bodies_cube COMMAND cool_bodies_bill_test -tc=cube) + add_test(NAME cool_bodies_cross COMMAND cool_bodies_bill_test -tc=cross) + add_test(NAME cool_bodies_birkhoff COMMAND cool_bodies_bill_test -tc=birk) + #add_test(NAME cool_bodies_prod_simplex COMMAND cool_bodies_bill_test -tc=prod_simplex) + add_test(NAME cool_bodies_simplex COMMAND cool_bodies_bill_test -tc=simplex) + add_test(NAME cool_bodies_skinny_cube COMMAND cool_bodies_bill_test -tc=skinny_cube) + #add_test(NAME round_skinny_cube COMMAND rounding_test -tc=round_skinny_cube) #add_test(NAME round_rot_skinny_cube COMMAND rounding_test -tc=round_rot_skinny_cube) @@ -121,6 +129,7 @@ else () TARGET_LINK_LIBRARIES(VpolyCV_test ${LP_SOLVE}) TARGET_LINK_LIBRARIES(VpolyVol_test ${LP_SOLVE}) TARGET_LINK_LIBRARIES(ZonotopeVol_test ${LP_SOLVE}) + TARGET_LINK_LIBRARIES(cool_bodies_bill_test ${LP_SOLVE}) #TARGET_LINK_LIBRARIES(ZonotopeVolCV_test ${LP_SOLVE}) endif() diff --git a/test/VpolyCV_test.cpp b/test/VpolyCV_test.cpp index 87fa8e490..6ee673af3 100644 --- a/test/VpolyCV_test.cpp +++ b/test/VpolyCV_test.cpp @@ -13,7 +13,7 @@ #include "random/normal_distribution.hpp" #include "random/uniform_real_distribution.hpp" #include "volume.h" -#include "polytope_generators.h" +#include "known_polytope_generators.h" #include #include diff --git a/test/VpolyVol_test.cpp b/test/VpolyVol_test.cpp index 8f5ba3c9d..34e8a6c3a 100644 --- a/test/VpolyVol_test.cpp +++ b/test/VpolyVol_test.cpp @@ -13,7 +13,7 @@ #include "random/normal_distribution.hpp" #include "random/uniform_real_distribution.hpp" #include "volume.h" -#include "polytope_generators.h" +#include "known_polytope_generators.h" #include template diff --git a/test/ZonotopeVolCV_test.cpp b/test/ZonotopeVolCV_test.cpp index 0f9eff163..b65604265 100644 --- a/test/ZonotopeVolCV_test.cpp +++ b/test/ZonotopeVolCV_test.cpp @@ -9,7 +9,7 @@ #include #include "Eigen/Eigen" #include "volume.h" -#include "polytope_generators.h" +#include "z_polytopes_gen.h" #include "exact_vols.h" #include @@ -26,7 +26,7 @@ void test_zono_volume(int n, int m, NT tolerance = 0.3) typedef typename Kernel::Point Point; typedef boost::mt19937 RNGType; typedef Zonotope Zonotope; - Zonotope ZP = gen_zonotope(n, m); + Zonotope ZP = gen_zonotope_uniform(n, m); // Setup the parameters int walk_len=1; diff --git a/test/ZonotopeVol_test.cpp b/test/ZonotopeVol_test.cpp index af35577f8..19fce904a 100644 --- a/test/ZonotopeVol_test.cpp +++ b/test/ZonotopeVol_test.cpp @@ -13,7 +13,7 @@ #include "random/normal_distribution.hpp" #include "random/uniform_real_distribution.hpp" #include "volume.h" -#include "polytope_generators.h" +#include "z_polytopes_gen.h" #include "exact_vols.h" #include @@ -30,7 +30,7 @@ void test_zono_volume(int n, int m, NT tolerance = 0.15) typedef typename Kernel::Point Point; typedef boost::mt19937 RNGType; typedef Zonotope Zonotope; - Zonotope ZP = gen_zonotope(n, m); + Zonotope ZP = gen_zonotope_uniform(n, m); // Setup the parameters int walk_len=10 + n/10; diff --git a/test/chebychev_test.cpp b/test/chebychev_test.cpp index 942581f0b..27a03e52d 100644 --- a/test/chebychev_test.cpp +++ b/test/chebychev_test.cpp @@ -15,7 +15,7 @@ #include "random/uniform_real_distribution.hpp" #include "volume.h" #include "misc.h" -#include "polytope_generators.h" +#include "known_polytope_generators.h" template NT factorial(NT n) diff --git a/test/cooling_bodies_bill_test.cpp b/test/cooling_bodies_bill_test.cpp new file mode 100644 index 000000000..fef5faffe --- /dev/null +++ b/test/cooling_bodies_bill_test.cpp @@ -0,0 +1,283 @@ +#include "doctest.h" +#include +#include "Eigen/Eigen" +#define VOLESTI_DEBUG +#include +#include "random.hpp" +#include "random/uniform_int.hpp" +#include "random/normal_distribution.hpp" +#include "random/uniform_real_distribution.hpp" +#include +#include +#include +#include +#include +#include +#include "cartesian_geom/cartesian_kernel.h" +#include "vars.h" +#include "hpolytope.h" +#include "vpolytope.h" +#include "zpolytope.h" +#include "ball.h" +#include "ballintersectconvex.h" +#include "vpolyintersectvpoly.h" +#include "samplers.h" +#include "rounding.h" +#include "rotating.h" +#include "linear_extensions.h" +#include "gaussian_annealing.h" +#include "cooling_balls.h" +#include "cooling_hpoly.h" +#include "misc.h" +#include "known_polytope_generators.h" +#include + +template +NT factorial(NT n) +{ + return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n; +} + +template +void test_cool_bodies(Polytope &HP, NT expected, NT tolerance=0.1, bool round = false, NT diam = -1.0) +{ + + typedef typename Polytope::PolytopePoint Point; + + // Setup the parameters + int n = HP.dimension(); + int walk_len=1; + int nexp=1, n_threads=1; + NT e=0.1, err=0.0000000001, diameter = diam, round_val = 1.0; + int rnum = std::pow(e,-2) * 400 * n * std::log(n); + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + RNGType rng(seed); + boost::normal_distribution<> rdist(0,1); + boost::random::uniform_real_distribution<>(urdist); + boost::random::uniform_real_distribution<> urdist1(-1,1); + + + std::pair InnerBall; + if (round) { + InnerBall = HP.ComputeInnerBall(); + vars var2(rnum,n,walk_len,n_threads,err,e,0,0,0,0,0.0,rng, + urdist,urdist1,-1.0,false,false,false,false,false,false,false,false,true); + std::pair res_round = rounding_min_ellipsoid(HP , InnerBall, var2); + round_val = res_round.first; + } + + InnerBall = HP.ComputeInnerBall(); + if(diameter < 0.0) { + HP.comp_diam(diameter, InnerBall.second); + if (round) diameter*=2.0; + } + + vars var(rnum,n,walk_len,n_threads,err,e,0,0,0,InnerBall.second,diameter,rng, + urdist,urdist1,-1.0,false,false,false,false,false,false,false,false,true); + + NT lb = 0.1, ub = 0.15, p = 0.75, rmax = 0.0, alpha = 0.2; + int W = 250, NNu = 140, nu =10; + bool win2 = false; + vars_ban var_ban(lb, ub, p, rmax, alpha, W, NNu, nu, win2); + + //Compute chebychev ball// + std::pair CheBall; + + // Estimate the volume + std::cout << "Number type: " << typeid(NT).name() << std::endl; + NT vol = 0; + unsigned int const num_of_exp = 10; + for (unsigned int i=0; i(10, false); + test_cool_bodies(P, 1024.0); + + std::cout << "--- Testing volume of H-cube20" << std::endl; + P = gen_cube(20, false); + test_cool_bodies(P, 1048576.0); + + std::cout << "--- Testing volume of H-cube30" << std::endl; + P = gen_cube(30, false); + test_cool_bodies(P, 1073742000.0); +} + +template +void call_test_cross(){ + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef boost::mt19937 RNGType; + typedef HPolytope Hpolytope; + + std::cout << "--- Testing volume of H-cross10" << std::endl; + Hpolytope P = gen_cross(10, false); + test_cool_bodies(P, 0.0002821869); +} + +template +void call_test_birk() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef boost::mt19937 RNGType; + typedef HPolytope Hpolytope; + Hpolytope P; + + std::cout << "--- Testing volume of H-birk3" << std::endl; + std::ifstream inp; + std::vector > Pin; + inp.open("../R-proj/inst/extdata/birk3.ine",std::ifstream::in); + read_pointset(inp,Pin); + P.init(Pin); + test_cool_bodies(P, 0.125); + + std::cout << "--- Testing volume of H-birk4" << std::endl; + std::ifstream inp2; + std::vector > Pin2; + inp2.open("../R-proj/inst/extdata/birk4.ine",std::ifstream::in); + read_pointset(inp2,Pin2); + P.init(Pin2); + test_cool_bodies(P, 0.000970018, 0.2); + + std::cout << "--- Testing volume of H-birk5" << std::endl; + std::ifstream inp3; + std::vector > Pin3; + inp3.open("../R-proj/inst/extdata/birk5.ine",std::ifstream::in); + read_pointset(inp3,Pin3); + P.init(Pin3); + test_cool_bodies(P, 0.000000225, 0.2); + + //std::cout << "--- Testing volume of H-birk6" << std::endl; + //std::ifstream inp4; + //std::vector > Pin4; + //inp4.open("../R-proj/inst/extdata/birk6.ine",std::ifstream::in); + //read_pointset(inp4,Pin4); + //P.init(Pin4); + //test_volume(P, 0.0000000000009455459196, 0.5); +} + +template +void call_test_prod_simplex() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef boost::mt19937 RNGType; + typedef HPolytope Hpolytope; + Hpolytope P; + + std::cout << "--- Testing volume of H-prod_simplex5" << std::endl; + P = gen_prod_simplex(5); + test_cool_bodies(P, std::pow(1.0 / factorial(5.0), 2.0),0.15, true); + + std::cout << "--- Testing volume of H-prod_simplex10" << std::endl; + P = gen_prod_simplex(10); + test_cool_bodies(P, std::pow(1.0 / factorial(10.0), 2.0), 0.15, true); + + //std::cout << "--- Testing volume of H-prod_simplex15" << std::endl; + //P = gen_prod_simplex(15); + //test_volume(P, std::pow(1.0 / factorial(15.0), 2)); + + std::cout << "--- Testing volume of H-prod_simplex20" << std::endl; + P = gen_prod_simplex(20); + test_cool_bodies(P, std::pow(1.0 / factorial(20.0), 2.0), 0.2, true); +} + +template +void call_test_simplex() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef boost::mt19937 RNGType; + typedef HPolytope Hpolytope; + Hpolytope P; + + std::cout << "--- Testing volume of H-simplex10" << std::endl; + P = gen_simplex(10, false); + test_cool_bodies(P, 1.0 / factorial(10.0), 0.1, false, 1.5); + + std::cout << "--- Testing volume of H-simplex20" << std::endl; + P = gen_simplex(20, false); + test_cool_bodies(P, 1.0 / factorial(20.0), 0.1, false, 1.5); + + std::cout << "--- Testing volume of H-simplex30" << std::endl; + P = gen_simplex(30, false); + test_cool_bodies(P, 1.0 / factorial(30.0), 0.1, false, 1.5); + + //std::cout << "--- Testing volume of H-simplex40" << std::endl; + //P = gen_simplex(40, false); + //test_volume(P, 1.0 / factorial(40.0)); + + //std::cout << "--- Testing volume of H-simplex50" << std::endl; + //P = gen_simplex(50, false); + //test_volume(P, 1.0 / factorial(50.0)); +} + +template +void call_test_skinny_cube() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef boost::mt19937 RNGType; + typedef HPolytope Hpolytope; + Hpolytope P; + + std::cout << "--- Testing volume of H-skinny_cube10" << std::endl; + P = gen_skinny_cube(10); + test_cool_bodies(P, 102400.0, 0.1, true); + + //std::cout << "--- Testing volume of H-skinny_cube20" << std::endl; + //P = gen_skinny_cube(20); + //test_volume(P, 104857600.0); +} + + +TEST_CASE("cube") { +call_test_cube(); +//call_test_cube(); +//call_test_cube(); +} + +TEST_CASE("cross") { +call_test_cross(); +//call_test_cross(); +//call_test_cross(); +} + +TEST_CASE("birk") { +call_test_birk(); +//call_test_birk(); +//call_test_birk(); +} + +TEST_CASE("prod_simplex") { +call_test_prod_simplex(); +//call_test_prod_simplex(); +//call_test_prod_simplex(); +} + +TEST_CASE("simplex") { +call_test_simplex(); +//call_test_simplex(); +//call_test_simplex(); +} + +TEST_CASE("skinny_cube") { +call_test_skinny_cube(); +//call_test_skinny_cube(); +//call_test_skinny_cube(); +} + diff --git a/test/generator.cpp b/test/generator.cpp index 0a7b63c57..799e9a4eb 100644 --- a/test/generator.cpp +++ b/test/generator.cpp @@ -12,84 +12,81 @@ #include "random/uniform_int.hpp" #include "random/normal_distribution.hpp" #include "random/uniform_real_distribution.hpp" -#include "convex_bodies/hpolytope.h" -#include "convex_bodies/vpolytope.h" -#include "convex_bodies/zpolytope.h" -#include "polytope_generators.h" +#include "vpolyoracles.h" +#include "hpolytope.h" +#include "vpolytope.h" +#include "zpolytope.h" +#include "known_polytope_generators.h" +#include "h_polytopes_gen.h" +#include "v_polytopes_gen.h" +#include "z_polytopes_gen.h" #include #include template -void create_txt(MT A, VT b, int kind, bool Vpoly) { +void create_txt(MT A, VT b, int kind, bool Vpoly, bool Zono) { int d = A.cols(), m = A.rows(); - std::string bar = "_"; - std::string ext = ".ext"; - std::string ine = ".ine"; - std::string name; + std::string bar = "_", ext = ".ext", ine = ".ine", name; std::ofstream outputFile; - if(Vpoly) { - if (kind == 0) { - std::string poly = "zonotope"; - name = poly + bar + std::to_string(d) + bar + std::to_string(m) + ext; - outputFile.open(name); - outputFile< 3 || kind <1) std::cout << "Wrong declaration of zonotope's generator, try -help" << std::endl; if (m > 0) { - Zonotope ZP = gen_zonotope(d, m); - create_txt(ZP.get_mat(), ZP.get_vec(), kind, true); + Zonotope ZP; + switch (kind) { + case 1: + ZP = gen_zonotope_uniform(d, m); + break; + case 2: + ZP = gen_zonotope_gaussian(d, m); + break; + case 3: + ZP = gen_zonotope_exponential(d, m); + break; + default: + kind = 1; + ZP = gen_zonotope_uniform(d, m); + break; + } + create_txt(ZP.get_mat(), ZP.get_vec(), kind, false, true); } else { std::cout << "Wrong inputs, try -help" << std::endl; exit(-1); } } else if (Hpoly) { Hpolytope HP; - if (cube) { - HP = gen_cube(d, false); - } else if (cross) { - HP = gen_cross(d, false); - } else if (simplex) { - HP = gen_simplex(d, false); - } else if (prod_simplex) { - HP = gen_prod_simplex(d); - } else if (skinny_cube) { - HP = gen_skinny_cube(d); - } else { - std::cout << "Wrong inputs, try -help" << std::endl; - exit(-1); + switch (kind) { + case 1: + HP = gen_cube(d, false); + break; + case 2: + HP = gen_cross(d, false); + break; + case 3: + HP = gen_simplex(d, false); + break; + case 4: + HP = gen_prod_simplex(d); + break; + case 5: + HP = gen_skinny_cube(d); + break; + case 6: + if (m > d + 1) { + HP = random_hpoly(d, m); + } else { + std::cout << "The number of facets has to be >= d+1" << std::endl; + exit(-1); + } + break; + default: + std::cout << "Wrong inputs, try -help" << std::endl; + exit(-1); } - create_txt(HP.get_mat(), HP.get_vec(), kind, false); + create_txt(HP.get_mat(), HP.get_vec(), kind, false, false); } else if (Vpoly) { Vpolytope VP; - if (cube) { - VP = gen_cube(d, true); - } else if (cross) { - VP = gen_cross(d, true); - } else if (simplex) { - VP = gen_simplex(d, true); - } else if (prod_simplex) { - std::cout<<"No prod_simplex in V-representation can be generated, try -help"<(d, true); + break; + case 2: + VP = gen_cross(d, true); + break; + case 3: + VP = gen_simplex(d, true); + break; + case 4: + if (!randv) { + std::cout<<"No prod_simplex in V-representation can be generated, try -help"< d) { + VP = random_vpoly_incube(d, m); + } else { + std::cout << "The number of vertices has to be >=d+1" << std::endl; + exit(-1); + } + break; + case 5: + if (!randv) { + std::cout<<"No skinny_cube in V-representation can be generated, try -help"< d) { + VP = random_vpoly(d, m); + } else { + std::cout << "The number of vertices has to be >=d+1" << std::endl; + exit(-1); + } + break; + default: + std::cout << "Wrong inputs, try -help" << std::endl; + exit(-1); } - create_txt(VP.get_mat(), VP.get_vec(), kind, true); + create_txt(VP.get_mat(), VP.get_vec(), kind, true, false); } else { - std::cout<<"Wrong inputs, try -help"<200\n"<< + "-cb : use Cooling Bodies algorithm for volume computation. This is the default choice for V-polytopes, Zonotopes and H-polytopes in dimension <=200\n"<< + "-sob : use Sequence of Balls algorithm for volume computation\n"<< "-w, --walk_len [walk_len] : the random walk length (default 1)\n"<< - "-rdhr : use random directions HnR, default is coordinate directions HnR\n" - "-e, --error epsilon : the goal error of approximation\n"<< - "-bw : use ball walk for sampling\n"<< + "-rdhr : use random directions HnR\n"<< + "-cdhr : use coordinate directions HnR. This is the default choice for H-polytopes\n"<< + "-biw : use Billiard walk. This is the default choice for V-polytopes and zonotopes\n"<< + "-L : the maximum length of the billiard trajectory\n"<< + "-baw : use ball walk\n"<< "-bwr : the radius of the ball walk (default r*chebychev_radius/sqrt(max(1.0, a_i)*dimension\n"<< - "-Win : the size of the open window for the ratios convergence\n"<< - "-C : a constant for the upper boud of variance/mean^2 in schedule annealing\n" - "-N : the number of points to sample in each step of schedule annealing. Default value N = 500*C + dimension^2/2\n"<< - "-frac : the fraction of the total error to spend in the first gaussian (default frac=0.1)\n"<< - "-ratio : parameter of schedule annealing, larger ratio means larger steps in schedule annealing (default 1-1/dimension)\n"<< + "-e, --error epsilon : the goal error of approximation\n"<< + "-win_len : the size of the open window for the ratios convergence (for CB and CG algorithms)\n"<< + "-C : a constant for the upper boud of variance/mean^2 in schedule annealing. The default values is 2 (for CG algorithm)\n" + "-N : the number of points to sample in each step of schedule annealing. The default value N = 500*C + dimension^2/2 (for CG algorithm)\n"<< + "-frac : the fraction of the total error to spend in the first gaussian. The default frac=0.1 (for CG algo)\n"<< + "-ratio : parameter of schedule annealing, larger ratio means larger steps in schedule annealing. The default 1-1/dimension (for CG algorithm)\n"<< + "-lb : lower bound for the volume ratios in CB algorithm. The default values is 0.1\n"<< + "-ub : upper bound for the volume ratios in CB algorithm. The default values is 0.15\n"<< + "-alpha : the significance level for the statistical test in CB algorithm\n"<< + "-nu : the degrees of freedom of t-student to use in the t-tests in CB algorithm. The default value is 10\n"<< + "-nuN : the degrees of freedom of t-student to use in the t-tests and the number of samples to perform the statistical tests in CB algorithm. The default values is 10 and 125 (when -biw is used) or 120 + d*d/10 (when other random walks are used)\n"<< std::endl; return 0; } @@ -157,7 +162,7 @@ int main(const int argc, const char** argv) exactvol = atof(argv[++i]); correct = true; } - if(!strcmp(argv[i],"-exact_zono")){ + if(!strcmp(argv[i],"-exact_vol")){ exact_zono = true; correct = true; } @@ -181,12 +186,12 @@ int main(const int argc, const char** argv) cdhr = true; correct = true; } - if(!strcmp(argv[i],"-BiW")){ + if(!strcmp(argv[i],"-biw")){ user_randwalk = true; billiard = true; correct = true; } - if(!strcmp(argv[i],"-BaW")){ + if(!strcmp(argv[i],"-baw")){ user_randwalk = true; ball_walk = true; correct = true; @@ -199,7 +204,7 @@ int main(const int argc, const char** argv) hpoly = true; correct = true; } - if(!strcmp(argv[i],"-WinLen")){ + if(!strcmp(argv[i],"-win_len")){ W = atof(argv[++i]); user_W = true; correct = true; @@ -335,7 +340,7 @@ int main(const int argc, const char** argv) set_error = true; correct = true; } - if(!strcmp(argv[i],"-w")||!strcmp(argv[i],"--walk_len")){ + if(!strcmp(argv[i],"-w")||!strcmp(argv[i],"-walk_len")){ walk_len = atof(argv[++i]); user_walk_len = true; correct = true; @@ -344,7 +349,7 @@ int main(const int argc, const char** argv) nexp = atof(argv[++i]); correct = true; } - if(!strcmp(argv[i],"-diameter")){ + if(!strcmp(argv[i],"-L")){ diameter = atof(argv[++i]); correct = true; } @@ -413,11 +418,59 @@ int main(const int argc, const char** argv) } if (exact_zono) { + if (!Zono) { + std::cout<<"Exact volume computation is available only for zonotopes."<(ZP); std::cout<<"Zonotope's exact volume = "< #include diff --git a/test/volume_test.cpp b/test/volume_test.cpp index d6d21fa23..cae626a3f 100644 --- a/test/volume_test.cpp +++ b/test/volume_test.cpp @@ -15,7 +15,7 @@ #include "random/uniform_real_distribution.hpp" #include "volume.h" #include "misc.h" -#include "polytope_generators.h" +#include "known_polytope_generators.h" #include template