Use billiard walk in SOB algorithm (#75)

* fix seed in rand_hpoly_generator * enable SOB algorithm with billiard walk * enable SOB algorithm with billiard walk in c++ interface
GeomScale · Mar 28, 2020 · d3bda40 · d3bda40
1 parent f10edad
commit d3bda40
Show file tree

Hide file tree

Showing 3 changed files with 28 additions and 39 deletions.
diff --git a/R-proj/src/volume.cpp b/R-proj/src/volume.cpp
@@ -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;
@@ -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))){

diff --git a/include/volume/volume.h b/include/volume/volume.h
@@ -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;
@@ -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
@@ -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
 
@@ -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());
 
@@ -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;
@@ -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
@@ -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;

diff --git a/test/vol.cpp b/test/vol.cpp
@@ -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) {