Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improvements in rounding and rotating R functions #76

Merged
merged 3 commits into from
Apr 10, 2020
Merged

Improvements in rounding and rotating R functions #76

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
3ce3ec5
add seed in rotate_polytope and update Rd files
TolisChal Apr 8, 2020
7892267
update R rounding function to return the linear map
TolisChal Apr 8, 2020
9068396
fix rounding in V-poly
TolisChal Apr 9, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 4 additions & 3 deletions R-proj/R/RcppExports.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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
Expand Down
7 changes: 4 additions & 3 deletions R-proj/R/rotate_polytope.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
#'
Expand All @@ -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
Expand All @@ -43,7 +44,7 @@ rotate_polytope <- function(P, T = NULL){

# remove first column
A = Mat[,-c(1)]
#A = Mat[,-c(1)]

type = P$type
if (type == 2) {
PP = Vpolytope$new(A)
Expand Down
7 changes: 4 additions & 3 deletions R-proj/R/round_polytope.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
}
2 changes: 1 addition & 1 deletion R-proj/man/poly_gen.Rd
View file
Open in desktop

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

4 changes: 3 additions & 1 deletion R-proj/man/rotate_polytope.Rd
View file
Open in desktop

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

4 changes: 3 additions & 1 deletion R-proj/man/rotating.Rd
View file
Open in desktop

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

2 changes: 1 addition & 1 deletion R-proj/src/poly_gen.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
19 changes: 12 additions & 7 deletions R-proj/src/rotating.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,13 +27,15 @@
//'
//' @param P A convex polytope (H-, V-polytope or a zonotope).
//' @param T Optional. A rotation matrix.
//' @param seed Optional. A fixed seed for the random linear map generator.
//'
//' @section warning:
//' Do not use this function.
//'
//' @return A matrix that describes the rotated polytope
// [[Rcpp::export]]
Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable<Rcpp::NumericMatrix> T = R_NilValue){
Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable<Rcpp::NumericMatrix> T = R_NilValue,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why not creating two overloads of rotating, one with fixed seed and one without ?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We use the same code structure for polytope generators as well.

Rcpp::Nullable<double> seed = R_NilValue){

typedef double NT;
typedef Cartesian<NT> Kernel;
Expand All @@ -51,6 +53,8 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable<Rcpp::NumericMat
unsigned int n = P.field("dimension");
int type = P.field("type");

double seed2 = (!seed.isNotNull()) ? std::numeric_limits<double>::signaling_NaN() : Rcpp::as<double>(seed);

switch (type) {
case 1: {
// Hpolytope
Expand All @@ -60,7 +64,7 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable<Rcpp::NumericMat
TransorfMat = Rcpp::as<MT>(T);
HP.linear_transformIt(TransorfMat.inverse());
} else {
TransorfMat = rotating < MT > (HP);
TransorfMat = rotating < MT > (HP, seed2);
}
Mat = extractMatPoly(HP);
break;
Expand All @@ -73,7 +77,7 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable<Rcpp::NumericMat
TransorfMat = Rcpp::as<MT>(T);
VP.linear_transformIt(TransorfMat.inverse());
} else {
TransorfMat = rotating < MT > (VP);
TransorfMat = rotating < MT > (VP, seed2);
}
Mat = extractMatPoly(VP);
break;
Expand All @@ -86,13 +90,14 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable<Rcpp::NumericMat
TransorfMat = Rcpp::as<MT>(T);
ZP.linear_transformIt(TransorfMat.inverse());
} else {
TransorfMat = rotating < MT > (ZP);
TransorfMat = rotating < MT > (ZP, seed2);
}
Mat = extractMatPoly(ZP);
break;
}
case 4: {
Vpolytope VP1;
throw Rcpp::exception("volesti does not support rotation for this representation currently.");
/*Vpolytope VP1;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is this piece of code going to be used in the future?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I removed this part.

Vpolytope VP2;
InterVP VPcVP;
VP1.init(n, Rcpp::as<MT>(P.field("V1")), VT::Ones(Rcpp::as<MT>(P.field("V1")).rows()));
Expand All @@ -102,8 +107,8 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable<Rcpp::NumericMat
TransorfMat = Rcpp::as<MT>(T);
VPcVP.linear_transformIt(TransorfMat.inverse());
} else {
TransorfMat = rotating < MT > (VPcVP);
}
TransorfMat = rotating < MT > (VPcVP, seed2);
}*/
}
}

Expand Down
12 changes: 7 additions & 5 deletions R-proj/src/rounding.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -126,26 +126,28 @@ Rcpp::List rounding (Rcpp::Reference P){
// initialization
vars<NT, RNGType> var(rnum,n,walkL,1,0.0,0.0,0,0.0,0,InnerBall.second,diam,rng,urdist,urdist1,
delta,verbose,rand_only,false,NN,birk,ball_walk,cdhr,rdhr,billiard);
std::pair <NT, NT> round_res;
std::pair< std::pair<MT, VT>, NT > round_res;

switch (type) {
case 1: {
round_res = rounding_min_ellipsoid(HP, InnerBall, var);
round_res = rounding_min_ellipsoid<MT, VT>(HP, InnerBall, var);
Mat = extractMatPoly(HP);
break;
}
case 2: {
round_res = rounding_min_ellipsoid(VP, InnerBall, var);
round_res = rounding_min_ellipsoid<MT, VT>(VP, InnerBall, var);
Mat = extractMatPoly(VP);
break;
}
case 3: {
round_res = rounding_min_ellipsoid(ZP, InnerBall, var);
round_res = rounding_min_ellipsoid<MT, VT>(ZP, InnerBall, var);
Mat = extractMatPoly(ZP);
break;
}
}

return Rcpp::List::create(Rcpp::Named("Mat") = Mat , Rcpp::Named("round_value") = round_res.first);
return Rcpp::List::create(Rcpp::Named("Mat") = Mat, Rcpp::Named("T") = Rcpp::wrap(round_res.first.first),
Rcpp::Named("shift") = Rcpp::wrap(round_res.first.second),
Rcpp::Named("round_value") = round_res.second);

}
35 changes: 19 additions & 16 deletions include/convex_bodies/vpolytope.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -272,27 +272,30 @@ class VPolytope{


std::list<Point> randPoints;
get_points_for_rounding(randPoints);
if (!get_points_for_rounding(randPoints)) {
center = get_mean_of_vertices();
} else {

boost::numeric::ublas::matrix<double> Ap(_d,randPoints.size());
typename std::list<Point>::iterator rpit=randPoints.begin();
boost::numeric::ublas::matrix<double> Ap(_d,randPoints.size());
typename std::list<Point>::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<double> Q(_d, _d);
boost::numeric::ublas::vector<double> c2(_d);
size_t w=1000;
KhachiyanAlgo(Ap,0.01,w,Q,c2); // call Khachiyan algorithm
boost::numeric::ublas::matrix<double> Q(_d, _d);
boost::numeric::ublas::vector<double> 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<NT,NT> res;
for (unsigned int i = 0; i < _d; ++i) {
Expand Down
13 changes: 9 additions & 4 deletions include/rotating.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,14 +12,19 @@


template <typename MT, typename Polytope>
MT rotating(Polytope &P){
MT rotating(Polytope &P, double seed = std::numeric_limits<double>::signaling_NaN()){
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same here, why not two overloads?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

overloaded


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();
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

needed?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

removed

//RNGType rng(seed);
unsigned int n = P.dimension();

// pick a random rotation
Expand Down
Loading