From 3ce3ec51b766c2d3e5e16a2bda7ec3ad3f7ec068 Mon Sep 17 00:00:00 2001 From: Tolis Date: Wed, 8 Apr 2020 23:38:59 +0300 Subject: [PATCH 1/3] add seed in rotate_polytope and update Rd files --- R-proj/R/RcppExports.R | 7 ++++--- R-proj/R/rotate_polytope.R | 5 +++-- R-proj/man/poly_gen.Rd | 2 +- R-proj/man/rotate_polytope.Rd | 4 +++- R-proj/man/rotating.Rd | 4 +++- R-proj/src/poly_gen.cpp | 2 +- R-proj/src/rotating.cpp | 14 +++++++++----- include/rotating.h | 13 +++++++++---- 8 files changed, 33 insertions(+), 18 deletions(-) diff --git a/R-proj/R/RcppExports.R b/R-proj/R/RcppExports.R index 25d8845b1..d15f28082 100644 --- a/R-proj/R/RcppExports.R +++ b/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 diff --git a/R-proj/R/rotate_polytope.R b/R-proj/R/rotate_polytope.R index 3a20276c0..1d50bde5f 100644 --- a/R-proj/R/rotate_polytope.R +++ b/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 diff --git a/R-proj/man/poly_gen.Rd b/R-proj/man/poly_gen.Rd index e16aae1ee..30f922032 100644 --- a/R-proj/man/poly_gen.Rd +++ b/R-proj/man/poly_gen.Rd @@ -17,7 +17,7 @@ poly_gen(kind_gen, Vpoly_gen, Zono_gen, dim_gen, m_gen, seed = NULL) \item{m_gen}{An integer to declare the number of generators for the requested random zonotope or the number of vertices for a V-polytope.} -\item{seed}{This variable can be used to set a seed for the random polytope generator.} +\item{seed}{Optional. A fixed seed for the random polytope generator.} } \value{ A numerical matrix describing the requested polytope diff --git a/R-proj/man/rotate_polytope.Rd b/R-proj/man/rotate_polytope.Rd index c53583536..70c9c3adf 100644 --- a/R-proj/man/rotate_polytope.Rd +++ b/R-proj/man/rotate_polytope.Rd @@ -4,12 +4,14 @@ \alias{rotate_polytope} \title{Apply a random rotation to a convex polytope (H-polytope, V-polytope, zonotope or intersection of two V-polytopes)} \usage{ -rotate_polytope(P, T = NULL) +rotate_polytope(P, T = NULL, seed = NULL) } \arguments{ \item{P}{A convex polytope. It is an object from class (a) Hpolytope, (b) Vpolytope, (c) Zonotope, (d) intersection of two V-polytopes.} \item{T}{Optional. A \eqn{d\times d} rotation matrix.} + +\item{seed}{Optional. A fixed seed for the random linear map generator.} } \value{ A list that contains the rotated polytope and the matrix \eqn{T} of the linear transformation. diff --git a/R-proj/man/rotating.Rd b/R-proj/man/rotating.Rd index f3514fb09..4e96293f5 100644 --- a/R-proj/man/rotating.Rd +++ b/R-proj/man/rotating.Rd @@ -4,12 +4,14 @@ \alias{rotating} \title{An internal Rccp function for the random rotation of a convex polytope} \usage{ -rotating(P, T = NULL) +rotating(P, T = NULL, seed = NULL) } \arguments{ \item{P}{A convex polytope (H-, V-polytope or a zonotope).} \item{T}{Optional. A rotation matrix.} + +\item{seed}{Optional. A fixed seed for the random linear map generator.} } \value{ A matrix that describes the rotated polytope diff --git a/R-proj/src/poly_gen.cpp b/R-proj/src/poly_gen.cpp index 6268abe3c..0a0ff9150 100644 --- a/R-proj/src/poly_gen.cpp +++ b/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. diff --git a/R-proj/src/rotating.cpp b/R-proj/src/rotating.cpp index 002be1129..84946a497 100644 --- a/R-proj/src/rotating.cpp +++ b/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 T = R_NilValue){ +Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable T = R_NilValue, + Rcpp::Nullable seed = R_NilValue){ typedef double NT; typedef Cartesian Kernel; @@ -51,6 +53,8 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable::signaling_NaN() : Rcpp::as(seed); + switch (type) { case 1: { // Hpolytope @@ -60,7 +64,7 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable(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(T); VP.linear_transformIt(TransorfMat.inverse()); } else { - TransorfMat = rotating < MT > (VP); + TransorfMat = rotating < MT > (VP, seed2); } Mat = extractMatPoly(VP); break; @@ -86,7 +90,7 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable(T); ZP.linear_transformIt(TransorfMat.inverse()); } else { - TransorfMat = rotating < MT > (ZP); + TransorfMat = rotating < MT > (ZP, seed2); } Mat = extractMatPoly(ZP); break; @@ -102,7 +106,7 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable(T); VPcVP.linear_transformIt(TransorfMat.inverse()); } else { - TransorfMat = rotating < MT > (VPcVP); + TransorfMat = rotating < MT > (VPcVP, seed2); } } } diff --git a/include/rotating.h b/include/rotating.h index 366b8c5f1..641a8d934 100644 --- a/include/rotating.h +++ b/include/rotating.h @@ -12,14 +12,19 @@ template -MT rotating(Polytope &P){ +MT rotating(Polytope &P, double seed = std::numeric_limits::signaling_NaN()){ typedef boost::mt19937 RNGType; - //typedef typename Polytope::MT MT; + unsigned rng_seed = std::chrono::system_clock::now().time_since_epoch().count(); + RNGType rng(rng_seed); + if (!std::isnan(seed)) { + unsigned rng_seed = seed; + rng.seed(rng_seed); + } boost::random::uniform_real_distribution<> urdist(-1.0, 1.0); - unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); - RNGType rng(seed); + //unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + //RNGType rng(seed); unsigned int n = P.dimension(); // pick a random rotation From 7892267cd71be7e8f34f052ba06904ef85647410 Mon Sep 17 00:00:00 2001 From: Tolis Date: Thu, 9 Apr 2020 01:22:25 +0300 Subject: [PATCH 2/3] update R rounding function to return the linear map --- R-proj/R/rotate_polytope.R | 2 +- R-proj/R/round_polytope.R | 7 ++-- R-proj/src/rotating.cpp | 5 ++- R-proj/src/rounding.cpp | 12 +++--- include/rounding.h | 79 ++++++++++++++++++++++++++++++++++++++ 5 files changed, 94 insertions(+), 11 deletions(-) diff --git a/R-proj/R/rotate_polytope.R b/R-proj/R/rotate_polytope.R index 1d50bde5f..12ce246b8 100644 --- a/R-proj/R/rotate_polytope.R +++ b/R-proj/R/rotate_polytope.R @@ -44,7 +44,7 @@ rotate_polytope <- function(P, T = NULL, seed = NULL){ # remove first column A = Mat[,-c(1)] - #A = Mat[,-c(1)] + type = P$type if (type == 2) { PP = Vpolytope$new(A) diff --git a/R-proj/R/round_polytope.R b/R-proj/R/round_polytope.R index 524aab22b..562f91a83 100644 --- a/R-proj/R/round_polytope.R +++ b/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) } diff --git a/R-proj/src/rotating.cpp b/R-proj/src/rotating.cpp index 84946a497..786ce09e6 100644 --- a/R-proj/src/rotating.cpp +++ b/R-proj/src/rotating.cpp @@ -96,7 +96,8 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable(P.field("V1")), VT::Ones(Rcpp::as(P.field("V1")).rows())); @@ -107,7 +108,7 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable (VPcVP, seed2); - } + }*/ } } diff --git a/R-proj/src/rounding.cpp b/R-proj/src/rounding.cpp index 5e16d8481..903ae64bf 100644 --- a/R-proj/src/rounding.cpp +++ b/R-proj/src/rounding.cpp @@ -126,26 +126,28 @@ Rcpp::List rounding (Rcpp::Reference P){ // initialization vars 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 round_res; + std::pair< std::pair, NT > round_res; switch (type) { case 1: { - round_res = rounding_min_ellipsoid(HP, InnerBall, var); + round_res = rounding_min_ellipsoid(HP, InnerBall, var); Mat = extractMatPoly(HP); break; } case 2: { - round_res = rounding_min_ellipsoid(VP, InnerBall, var); + round_res = rounding_min_ellipsoid(VP, InnerBall, var); Mat = extractMatPoly(VP); break; } case 3: { - round_res = rounding_min_ellipsoid(ZP, InnerBall, var); + round_res = rounding_min_ellipsoid(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); } diff --git a/include/rounding.h b/include/rounding.h index 63cd369f6..266517b5a 100644 --- a/include/rounding.h +++ b/include/rounding.h @@ -99,6 +99,85 @@ std::pair rounding_min_ellipsoid(Polytope &P , const std::pair +std::pair< std::pair, NT > rounding_min_ellipsoid(Polytope &P , const std::pair &InnerBall, const Parameters &var) { + + //typedef typename Polytope::VT VT; + typedef typename Parameters::RNGType RNGType; + + unsigned int n=var.n, walk_len=var.walk_steps, i, j = 0; + Point c = InnerBall.first; + NT radius = InnerBall.second; + std::list randPoints; //ds for storing rand points + if (!P.get_points_for_rounding(randPoints)) { // If P is a V-polytope then it will store its vertices in randPoints + // If P is not a V-Polytope or number_of_vertices>20*domension + // 2. Generate the first random point in P + // Perform random walk on random point in the Chebychev ball + Point p = get_point_in_Dsphere(n, radius); + p = p + c; + + //use a large walk length e.g. 1000 + rand_point_generator(P, p, 1, 10*n, randPoints, var); + // 3. Sample points from P + unsigned int num_of_samples = 10*n;//this is the number of sample points will used to compute min_ellipoid + randPoints.clear(); + if (var.bill_walk) { + rand_point_generator(P, p, num_of_samples, 5, randPoints, var); + } else { + rand_point_generator(P, p, num_of_samples, 10 + n * 5, randPoints, var); + } + } + + // Store points in a matrix to call Khachiyan algorithm for the minimum volume enclosing ellipsoid + boost::numeric::ublas::matrix Ap(n,randPoints.size()); + typename std::list::iterator rpit=randPoints.begin(); + + for ( ; rpit!=randPoints.end(); rpit++, j++) { + for (i=0 ; idimension(); i++){ + Ap(i,j)=double((*rpit)[i]); + } + } + boost::numeric::ublas::matrix Q(n,n); //TODO: remove dependence on ublas and copy to eigen + boost::numeric::ublas::vector c2(n); + size_t w=1000; + KhachiyanAlgo(Ap,0.01,w,Q,c2); // call Khachiyan algorithm + + MT E(n,n); + VT e(n); + + //Get ellipsoid matrix and center as Eigen objects + for(unsigned int i=0; i eigensolver(E); + NT rel = std::real(eigensolver.eigenvalues()[0]); + NT Rel = std::real(eigensolver.eigenvalues()[0]); + for(unsigned int i=1; iRel) Rel=std::real(eigensolver.eigenvalues()[i]); + } + + Eigen::LLT lltOfA(E); // compute the Cholesky decomposition of E + MT L = lltOfA.matrixL(); // retrieve factor L in the decomposition + + //Shift polytope in order to contain the origin (center of the ellipsoid) + P.shift(e); + + MT L_1 = L.inverse(); + P.linear_transformIt(L_1.transpose()); + + return std::pair< std::pair, NT > (std::pair(L_1, e), L_1.determinant()); +} + + template void get_vpoly_center(Polytope &P) { From 9068396f5a2f3498709b156fa6819110703da830 Mon Sep 17 00:00:00 2001 From: Tolis Date: Fri, 10 Apr 2020 00:15:59 +0300 Subject: [PATCH 3/3] fix rounding in V-poly --- include/convex_bodies/vpolytope.h | 35 +++++++++++----------- include/rounding.h | 48 ++++++++++++++++--------------- 2 files changed, 44 insertions(+), 39 deletions(-) diff --git a/include/convex_bodies/vpolytope.h b/include/convex_bodies/vpolytope.h index 3182246bd..804bf0166 100644 --- a/include/convex_bodies/vpolytope.h +++ b/include/convex_bodies/vpolytope.h @@ -272,27 +272,30 @@ class VPolytope{ std::list randPoints; - get_points_for_rounding(randPoints); + if (!get_points_for_rounding(randPoints)) { + center = get_mean_of_vertices(); + } else { - boost::numeric::ublas::matrix Ap(_d,randPoints.size()); - typename std::list::iterator rpit=randPoints.begin(); + boost::numeric::ublas::matrix Ap(_d,randPoints.size()); + typename std::list::iterator rpit=randPoints.begin(); - unsigned int i, j = 0; - for ( ; rpit!=randPoints.end(); rpit++, j++) { - const NT* point_data = rpit->getCoefficients().data(); + unsigned int i, j = 0; + for ( ; rpit!=randPoints.end(); rpit++, j++) { + const NT* point_data = rpit->getCoefficients().data(); - for ( i=0; i < rpit->dimension(); i++){ - Ap(i,j)=double(*point_data); - point_data++; + for ( i=0; i < rpit->dimension(); i++){ + Ap(i,j)=double(*point_data); + point_data++; + } } - } - boost::numeric::ublas::matrix Q(_d, _d); - boost::numeric::ublas::vector c2(_d); - size_t w=1000; - KhachiyanAlgo(Ap,0.01,w,Q,c2); // call Khachiyan algorithm + boost::numeric::ublas::matrix Q(_d, _d); + boost::numeric::ublas::vector c2(_d); + size_t w=1000; + KhachiyanAlgo(Ap,0.01,w,Q,c2); // call Khachiyan algorithm - //Get ellipsoid matrix and center as Eigen objects - for(unsigned int i=0; i<_d; i++) center.set_coord(i, NT(c2(i))); + //Get ellipsoid matrix and center as Eigen objects + for(unsigned int i=0; i<_d; i++) center.set_coord(i, NT(c2(i))); + } std::pair res; for (unsigned int i = 0; i < _d; ++i) { diff --git a/include/rounding.h b/include/rounding.h index 266517b5a..fc0dbfcf9 100644 --- a/include/rounding.h +++ b/include/rounding.h @@ -189,33 +189,35 @@ void get_vpoly_center(Polytope &P) { unsigned int n = P.dimension(); std::list randPoints; //ds for storing rand points - P.get_points_for_rounding(randPoints); - - boost::numeric::ublas::matrix Ap(n,randPoints.size()); - typename std::list::iterator rpit=randPoints.begin(); - for (int j=0 ; rpit!=randPoints.end(); rpit++, j++) { - for (int i=0 ; idimension(); i++){ - Ap(i,j)=double((*rpit)[i]); + if (!P.get_points_for_rounding(randPoints)) { + P.shift(P.get_mean_of_vertices().getCoefficients()); + } else { + + boost::numeric::ublas::matrix Ap(n,randPoints.size()); + typename std::list::iterator rpit=randPoints.begin(); + for (int j=0 ; rpit!=randPoints.end(); rpit++, j++) { + for (int i=0 ; idimension(); i++){ + Ap(i,j)=double((*rpit)[i]); + } } - } - boost::numeric::ublas::matrix Q(n,n); - boost::numeric::ublas::vector c2(n); - size_t w=1000; - KhachiyanAlgo(Ap,0.01,w,Q,c2); // call Khachiyan algorithm - - MT E(n,n); - VT e(n); - - //Get ellipsoid matrix and center as Eigen objects - for(unsigned int i=0; i Q(n,n); + boost::numeric::ublas::vector c2(n); + size_t w=1000; + KhachiyanAlgo(Ap,0.01,w,Q,c2); // call Khachiyan algorithm + + //MT E(n,n); + VT e(n); + + //Get ellipsoid matrix and center as Eigen objects + for(unsigned int i=0; i