Skip to content

Commit

Permalink
Use billiard walk in SOB algorithm (#75)
Browse files Browse the repository at this point in the history 
* fix seed in rand_hpoly_generator

* enable SOB algorithm with billiard walk

* enable SOB algorithm with billiard walk in c++ interface
  • Loading branch information
@TolisChal @vissarion
TolisChal authored and vissarion committed Apr 4, 2020
1 parent fa05bcf commit 69664de
Show file tree
Hide file tree
Showing 10 changed files with 483 additions and 247 deletions.
13 changes: 2 additions & 11 deletions R-proj/src/volume.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -260,14 +260,6 @@ double volume (Rcpp::Reference P,
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;
Expand Down Expand Up @@ -345,9 +337,8 @@ double volume (Rcpp::Reference P,
alpha = Rcpp::as<NT>(Rcpp::as<Rcpp::List>(settings)["alpha"]);
cb_params++;
}
if (Rcpp::as<Rcpp::List>(settings).containsElementNamed("L")) {
diam = Rcpp::as<NT>(Rcpp::as<Rcpp::List>(settings)["L"]);
cb_params++;
if (Rcpp::as<Rcpp::List>(algo).containsElementNamed("L")) {
diam = Rcpp::as<NT>(Rcpp::as<Rcpp::List>(algo)["L"]);
}

if ((CB && cg_params > 0) || (CG && cb_params > 0) || (!CB & !CG && (cg_params > 0 || cb_params > 0))){
Expand Down
4 changes: 0 additions & 4 deletions include/cartesian_geom/point.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -170,13 +170,9 @@ class point

void print() const{
for(unsigned int i=0; i<d; i++){
#ifdef VOLESTI_DEBUG
std::cout<<coeffs(i)<<" ";
#endif
}
#ifdef VOLESTI_DEBUG
std::cout<<"\n";
#endif
}

};
Expand Down
2 changes: 1 addition & 1 deletion include/convex_bodies/ball.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@
template <typename Point>
struct Ball{
public:
typedef Point BallPoint;
typedef Point PointType;
typedef typename Point::FT NT;
typedef typename std::vector<NT>::iterator viterator;
typedef Eigen::Matrix<NT,Eigen::Dynamic,1> VT;
Expand Down
28 changes: 14 additions & 14 deletions include/convex_bodies/ballintersectconvex.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ class BallIntersectPolytope {
CBall B;
public:
typedef typename CBall::NT NT;
typedef typename CBall::BallPoint Point;
typedef typename CBall::PointType PointType;
typedef Eigen::Matrix<NT,Eigen::Dynamic,1> VT;

BallIntersectPolytope() {}
Expand All @@ -29,17 +29,17 @@ class BallIntersectPolytope {
Polytope first() const { return P; }
CBall second() const { return B; }

std::pair<Point,NT> InnerBall()
std::pair<PointType,NT> InnerBall()
{
return P.InnerBall();
}

NT ComputeDiameter()
{
return NT(4) * B.radius();
return NT(2) * B.radius();
}

int is_in(Point &p){
int is_in(PointType &p){
if(B.is_in(p)==-1)
return P.is_in(p);
return 0;
Expand All @@ -61,31 +61,31 @@ class BallIntersectPolytope {
comp_diam(diam);
}

std::pair<NT,NT> line_intersect(Point &r, Point &v) {
std::pair<NT,NT> line_intersect(PointType &r, PointType &v) {

std::pair <NT, NT> polypair = P.line_intersect(r, v);
std::pair <NT, NT> ballpair = B.line_intersect(r, v);
return std::pair<NT, NT>(std::min(polypair.first, ballpair.first),
std::max(polypair.second, ballpair.second));
}

std::pair<NT,NT> line_intersect(Point &r, Point &v, VT &Ar,
std::pair<NT,NT> line_intersect(PointType &r, PointType &v, VT &Ar,
VT &Av) {
std::pair <NT, NT> polypair = P.line_intersect(r, v, Ar, Av);
std::pair <NT, NT> ballpair = B.line_intersect(r, v);
return std::pair<NT, NT>(std::min(polypair.first, ballpair.first),
std::max(polypair.second, ballpair.second));
}

std::pair<NT,NT> line_intersect(Point &r, Point &v, VT &Ar, VT &Av, NT &lambda_prev) {
std::pair<NT,NT> line_intersect(PointType &r, PointType &v, VT &Ar, VT &Av, NT &lambda_prev) {

std::pair <NT, NT> polypair = P.line_intersect(r, v, Ar, Av, lambda_prev);
std::pair <NT, NT> ballpair = B.line_intersect(r, v);
return std::pair<NT, NT>(std::min(polypair.first, ballpair.first),
std::max(polypair.second, ballpair.second));
}

std::pair<NT,int> line_positive_intersect(Point &r, Point &v, VT &Ar, VT &Av) {
std::pair<NT,int> line_positive_intersect(PointType &r, PointType &v, VT &Ar, VT &Av) {

std::pair <NT, int> polypair = P.line_positive_intersect(r, v, Ar, Av);
std::pair <NT, int> ball_lambda = B.line_positive_intersect(r, v);
Expand All @@ -97,7 +97,7 @@ class BallIntersectPolytope {
}


std::pair<NT,int> line_positive_intersect(Point &r, Point &v, VT &Ar, VT &Av,
std::pair<NT,int> line_positive_intersect(PointType &r, PointType &v, VT &Ar, VT &Av,
NT &lambda_prev) {

std::pair <NT, int> polypair = P.line_positive_intersect(r, v, Ar, Av, lambda_prev);
Expand All @@ -110,7 +110,7 @@ class BallIntersectPolytope {
}

//First coordinate ray shooting intersecting convex body
std::pair<NT,NT> line_intersect_coord(Point &r,
std::pair<NT,NT> line_intersect_coord(PointType &r,
const unsigned int &rand_coord,
VT &lamdas) {

Expand All @@ -121,8 +121,8 @@ class BallIntersectPolytope {
}

//Not the first coordinate ray shooting intersecting convex body
std::pair<NT,NT> line_intersect_coord(Point &r,
const Point &r_prev,
std::pair<NT,NT> line_intersect_coord(PointType &r,
const PointType &r_prev,
const unsigned int &rand_coord,
const unsigned int &rand_coord_prev,
VT &lamdas) {
Expand All @@ -133,14 +133,14 @@ class BallIntersectPolytope {
std::max(polypair.second, ballpair.second));
}

std::pair<NT,NT> query_dual(const Point &p, const unsigned int &rand_coord) {
std::pair<NT,NT> query_dual(const PointType &p, const unsigned int &rand_coord) {
std::pair <NT, NT> polypair = P.query_dual(p, rand_coord);
std::pair <NT, NT> ballpair = B.line_intersect_coord(p, rand_coord);
return std::pair<NT, NT>(std::min(polypair.first, ballpair.first),
std::max(polypair.second, ballpair.second));
}

void compute_reflection (Point &v, Point &p, int &facet) {
void compute_reflection (PointType &v, PointType &p, int &facet) {

if (facet == P.num_of_hyperplanes()) {
B.compute_reflection(v, p, facet);
Expand Down
10 changes: 2 additions & 8 deletions include/convex_bodies/hpolytope.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@
template <typename Point>
class HPolytope{
public:
typedef Point PolytopePoint;
typedef Point PointType;
typedef typename Point::FT NT;
typedef typename std::vector<NT>::iterator viterator;
//using RowMatrixXd = Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
Expand All @@ -44,7 +44,7 @@ class HPolytope{

HPolytope()
{
typedef Point PolytopePoint;
typedef Point PointType;
typedef typename Point::FT NT;
inner_ball = ComputeChebychevBall<NT, Point>(A, b);
}
Expand Down Expand Up @@ -141,18 +141,12 @@ class HPolytope{

// print polytope in input format
void print() {
#ifdef VOLESTI_DEBUG
std::cout << " " << A.rows() << " " << _d << " float" << std::endl;
#endif
for (unsigned int i = 0; i < A.rows(); i++) {
for (unsigned int j = 0; j < _d; j++) {
#ifdef VOLESTI_DEBUG
std::cout << A(i, j) << " ";
#endif
}
#ifdef VOLESTI_DEBUG
std::cout << "<= " << b(i) << std::endl;
#endif
}
}

Expand Down
Loading

0 comments on commit 69664de

Please sign in to comment.