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

Use billiard walk in SOB algorithm #75

Merged
merged 6 commits into from
Mar 28, 2020
Merged

Use billiard walk in SOB algorithm #75

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
cce53a4
fix seed in rand_hpoly_generator
TolisChal Mar 25, 2020
f2d7349
Merge branch 'improve_interface' into fix_seed
TolisChal Mar 25, 2020
74c5f71
enable SOB algorithm with billiard walk
TolisChal Mar 27, 2020
aa96a62
enable SOB algorithm with billiard walk in c++ interface
TolisChal Mar 27, 2020
2862cdd
merge fix_seed branch
TolisChal Mar 27, 2020
cf4bc82
Merge branch 'develop' into sob_billiard
TolisChal Mar 28, 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
13 changes: 2 additions & 11 deletions R-proj/src/volume.cpp
View file
Open in desktop
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
50 changes: 26 additions & 24 deletions include/volume/volume.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,8 +54,9 @@ NT volume(Polytope &P,
RNGType &rng = var.rng;

//0. Get the Chebychev ball (largest inscribed ball) with center and radius
Point c=InnerBall.first;
NT radius=InnerBall.second;
Point c = InnerBall.first;
NT radius = InnerBall.second;
P.normalize();

//1. Rounding of the polytope if round=true
NT round_value=1;
Expand All @@ -73,16 +74,17 @@ NT volume(Polytope &P,
std::pair<Point,NT> res=P.ComputeInnerBall();
c=res.first; radius=res.second;
P.comp_diam(var.diameter, radius);
P.normalize();
if (var.ball_walk){
var.delta = 4.0 * radius / NT(n);
}
}

// Move the chebychev center to the origin and apply the same shifting to the polytope
P.shift(c.getCoefficients());
c=Point(n);
c = Point(n);

rnum=rnum/n_threads;
rnum = rnum / n_threads;
NT vol=0;

// Perform the procedure for a number of threads and then take the average
Expand All @@ -93,7 +95,7 @@ NT volume(Polytope &P,
if(print) std::cout<<"\nGenerate the first random point in P"<<std::endl;
#endif

Point p = get_point_on_Dsphere<RNGType , Point>(n, radius);
Point p = get_point_in_Dsphere<RNGType , Point>(n, radius);
std::list<Point> randPoints; //ds for storing rand points
//use a large walk length e.g. 1000

Expand All @@ -115,37 +117,36 @@ NT volume(Polytope &P,
// 4. Construct the sequence of balls
// 4a. compute the radius of the largest ball
NT current_dist, max_dist=NT(0);
for(typename std::list<Point>::iterator pit=randPoints.begin(); pit!=randPoints.end(); ++pit){
current_dist=(*pit).squared_length();
if(current_dist>max_dist){
max_dist=current_dist;
for(typename std::list<Point>::iterator pit=randPoints.begin(); pit!=randPoints.end(); ++pit) {
current_dist = (*pit).squared_length();
if (current_dist > max_dist) {
max_dist = current_dist;
}
}
max_dist=std::sqrt(max_dist);
max_dist = std::sqrt(max_dist);
#ifdef VOLESTI_DEBUG
if(print) std::cout<<"\nFurthest distance from Chebychev point= "<<max_dist<<std::endl;
#endif

//
// 4b. Number of balls
int nb1 = n * (std::log(radius)/std::log(2.0));
int nb2 = std::ceil(n * (std::log(max_dist)/std::log(2.0)));
int nb2 = std::ceil(n * (std::log(max_dist) / std::log(2.0)));

#ifdef VOLESTI_DEBUG
if(print) std::cout<<"\nConstructing the sequence of balls"<<std::endl;
#endif

std::vector<Ball> balls;

for(int i=nb1; i<=nb2; ++i){
for(int i=nb1; i<=nb2; ++i) {

if(i==nb1){
balls.push_back(Ball(c,radius*radius));
vol = (std::pow(M_PI,n/2.0)*(std::pow(balls[0].radius(), n) ) ) / (tgamma(n/2.0+1));
}else{
balls.push_back(Ball(c,std::pow(std::pow(2.0,NT(i)/NT(n)),2)));
if (i == nb1) {
balls.push_back(Ball(c, radius * radius));
vol = (std::pow(M_PI, n / 2.0) * (std::pow(balls[0].radius(), n))) / (tgamma(n / 2.0 + 1));
} else {
balls.push_back(Ball(c, std::pow(std::pow(2.0, NT(i) / NT(n)), 2)));
}

}
assert(!balls.empty());

Expand Down Expand Up @@ -187,9 +188,9 @@ NT volume(Polytope &P,

//keep the points in randPoints that fall in PBSmall
typename std::list<Point>::iterator rpit=randPoints.begin();
while(rpit!=randPoints.end()){
if (PBSmall.second().is_in(*rpit) == 0){//not in
rpit=randPoints.erase(rpit);
while(rpit!=randPoints.end()) {
if (PBSmall.second().is_in(*rpit) == 0) {//not in
rpit = randPoints.erase(rpit);
} else {
++nump_PBSmall;
++rpit;
Expand All @@ -207,10 +208,11 @@ NT volume(Polytope &P,
<<std::endl;
#endif

PBLarge.comp_diam(var.diameter);
//generate more random points in PBLarge to have "rnum" in total
rand_point_generator(PBLarge,p_gen,rnum-nump_PBLarge,walk_len,randPoints,PBSmall,nump_PBSmall,var);
rand_point_generator(PBLarge, p_gen, rnum-nump_PBLarge, walk_len, randPoints, PBSmall, nump_PBSmall, var);

vol *= NT(rnum)/NT(nump_PBSmall);
vol *= NT(rnum) / NT(nump_PBSmall);

#ifdef VOLESTI_DEBUG
if(print) std::cout<<nump_PBSmall<<"/"<<rnum<<" = "<<NT(rnum)/nump_PBSmall
Expand All @@ -221,7 +223,7 @@ NT volume(Polytope &P,
//don't continue in pairs of balls that are almost inside P, i.e. ratio ~= 2
}
}
vol=round_value*vol;
vol = round_value * vol;
#ifdef VOLESTI_DEBUG
if(print) std::cout<<"rand points = "<<rnum<<std::endl;
if(print) std::cout<<"walk len = "<<walk_len<<std::endl;
Expand Down
4 changes: 0 additions & 4 deletions test/vol.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -462,10 +462,6 @@ int main(const int argc, const char** argv)
} else {
cdhr = true;
}
} else if (!CB && !CG && billiard) {
std::cout<<"Billiard is not supported for SOB algorithm. CDHR is used."<<std::endl;
billiard = false;
cdhr = true;
} else if (CG && billiard) {
billiard = false;
if (Zono || Vpoly) {
Expand Down