Dikin, Vaidya, John walk (#88)

* implement dikin walk * implement john walk * implement vaidya walk * include random walk * add c++ wrappers for the ellipsoid random walks * update R interface for sampling and parameterized ellipsoid walks * update Rd files and roxygen comments * implement PR comments * use constructors in dikin,vaidya and john walks * improve implementations * improve computation of inner ball for H-polytopes * improve inner ball computation and improve tests * update termination criterions in rounding methods * update Rd file of Birkhoff R generator * improve univariate psrf implemetations * fix c++ tests * update parameters in rounding and minor improvements * fix gcc tests and improve birkhoff generator Rd file * fix c++ tests * update citations
GeomScale · Sep 18, 2020 · 0fad1ce · 0fad1ce
1 parent 375fa1f
commit 0fad1ce
Show file tree

Hide file tree

Showing 19 changed files with 1,005 additions and 25 deletions.
diff --git a/R-proj/R/RcppExports.R b/R-proj/R/RcppExports.R
@@ -168,7 +168,7 @@ geweke <- function(samples, frac_first = NULL, frac_last = NULL) {
 #'
 #' # compute an inscribed ball of the 3-dimensional unit cube in V-representation
 #' P = gen_cube(3, 'V')
-#' ball_vec = inner_ball(P)
+#' ball_vec = inner_ball(P, lpsolve = TRUE)
 #' @export
 inner_ball <- function(P, lpsolve = NULL) {
     .Call(`_volesti_inner_ball`, P, lpsolve)
@@ -276,12 +276,12 @@ rounding <- function(P, method = NULL, seed = NULL) {
 #' @param n The number of points that the function is going to sample from the convex polytope.
 #' @param random_walk Optional. A list that declares the random walk and some related parameters as follows:
 #' \itemize{
-#' \item{\code{walk} }{ A string to declare the random walk: i) \code{'CDHR'} for Coordinate Directions Hit-and-Run, ii) \code{'RDHR'} for Random Directions Hit-and-Run, iii) \code{'BaW'} for Ball Walk, iv) \code{'BiW'} for Billiard walk, v) \code{'BCDHR'} boundary sampling by keeping the extreme points of CDHR or vi) \code{'BRDHR'} boundary sampling by keeping the extreme points of RDHR. The default walk is \code{'BiW'} for the uniform distribution or \code{'CDHR'} for the Gaussian distribution.}
+#' \item{\code{walk} }{ A string to declare the random walk: i) \code{'CDHR'} for Coordinate Directions Hit-and-Run, ii) \code{'RDHR'} for Random Directions Hit-and-Run, iii) \code{'BaW'} for Ball Walk, iv) \code{'BiW'} for Billiard walk, v) \code{'dikin'} for dikin walk, vi) \code{'vaidya'} for vaidya walk, vii) \code{'john'} for john walk, viii) \code{'BCDHR'} boundary sampling by keeping the extreme points of CDHR or ix) \code{'BRDHR'} boundary sampling by keeping the extreme points of RDHR. The default walk is \code{'aBiW'} for the uniform distribution or \code{'CDHR'} for the Gaussian distribution and H-polytopes and \code{'BiW'} or \code{'RDHR'} for the same distributions and V-polytopes and zonotopes.}
 #' \item{\code{walk_length} }{ The number of the steps per generated point for the random walk. The default value is \eqn{1}.}
 #' \item{\code{nburns} }{ The number of points to burn before start sampling. The default value is \eqn{1}.}
 #' \item{\code{starting_point} }{ A \eqn{d}-dimensional numerical vector that declares a starting point in the interior of the polytope for the random walk. The default choice is the center of the ball as that one computed by the function \code{inner_ball()}.}
 #' \item{\code{BaW_rad} }{ The radius for the ball walk.}
-#' \item{\code{L} }{ The maximum length of the billiard trajectory.}
+#' \item{\code{L} }{ The maximum length of the billiard trajectory or the radius for the step of dikin, vaidya or john walk.}
 #' }
 #' @param distribution Optional. A list that declares the target density and some related parameters as follows:
 #' \itemize{
@@ -291,6 +291,15 @@ rounding <- function(P, method = NULL, seed = NULL) {
 #' }
 #' @param seed Optional. A fixed seed for the number generator.
 #'
+#' @references \cite{Robert L. Smith,
+#' \dQuote{Efficient Monte Carlo Procedures for Generating Points Uniformly Distributed Over Bounded Regions,} \emph{Operations Research,} 1984.},
+#'
+#' @references \cite{B.T. Polyak, E.N. Gryazina,
+#' \dQuote{Billiard walk - a new sampling algorithm for control and optimization,} \emph{IFAC Proceedings Volumes,} 2014.},
+#'
+#' @references \cite{Y. Chen, R. Dwivedi, M. J. Wainwright and B. Yu,
+#' \dQuote{Vaidya walk: A sampling algorithm based on the volumetric barrier,} \emph{55th Annual Allerton Conference on Communication, Control, and Computing,} 2017.}
+#'
 #' @return A \eqn{d\times n} matrix that contains, column-wise, the sampled points from the convex polytope P.
 #' @examples
 #' # uniform distribution from the 3d unit cube in H-representation using ball walk
@@ -377,12 +386,12 @@ loadSdpaFormatFile <- function(inputFile = NULL) {
 #' @return The approximation of the volume of a convex polytope.
 #' @examples
 #'
-#' # calling SOB algorithm for a H-polytope (3d unit simplex)
-#' HP = gen_cube(3,'H')
+#' # calling SOB algorithm for a H-polytope (5d unit simplex)
+#' HP = gen_cube(5,'H')
 #' vol = volume(HP)
 #'
-#' # calling CG algorithm for a V-polytope (2d simplex)
-#' VP = gen_simplex(2,'V')
+#' # calling CG algorithm for a V-polytope (3d simplex)
+#' VP = gen_simplex(3,'V')
 #' vol = volume(VP, settings = list("algorithm" = "CG"))
 #'
 #' # calling CG algorithm for a 2-dimensional zonotope defined as the Minkowski sum of 4 segments

diff --git a/R-proj/man/inner_ball.Rd b/R-proj/man/inner_ball.Rd
diff --git a/R-proj/man/sample_points.Rd b/R-proj/man/sample_points.Rd
diff --git a/R-proj/man/volume.Rd b/R-proj/man/volume.Rd
diff --git a/R-proj/src/sample_points.cpp b/R-proj/src/sample_points.cpp
@@ -7,6 +7,8 @@
 
 //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 and 2019 program.
 
+// Contributed and/or modified by Alexandros Manochis, as part of Google Summer of Code 2020 program.
+
 #include <Rcpp.h>
 #include <RcppEigen.h>
 #include <chrono>
@@ -19,7 +21,9 @@
 #include "volume/volume_cooling_gaussians.hpp"
 #include "sampling/sampling.hpp"
 
-enum random_walks {ball_walk, rdhr, cdhr, billiard, accelarated_billiard, brdhr, bcdhr};
+
+enum random_walks {ball_walk, rdhr, cdhr, billiard, accelarated_billiard, 
+                   dikin_walk, vaidya_walk, john_walk, brdhr, bcdhr};
 
 template <typename Polytope, typename RNGType, typename PointList, typename NT, typename Point>
 void sample_from_polytope(Polytope &P, int type, RNGType &rng, PointList &randPoints,
@@ -56,6 +60,33 @@ void sample_from_polytope(Polytope &P, int type, RNGType &rng, PointList &randPo
                                              StartingPoint, nburns);
         }
         break;
+    case vaidya_walk:
+        if (set_L) {
+            VaidyaWalk WalkType(L);
+            uniform_sampling(randPoints, P, rng, WalkType, walkL, numpoints, StartingPoint, nburns);
+        } else {
+            uniform_sampling<VaidyaWalk>(randPoints, P, rng, walkL, numpoints,
+                                     StartingPoint, nburns);
+        }
+        break;
+    case dikin_walk:
+        if (set_L) {
+            DikinWalk WalkType(L);
+            uniform_sampling(randPoints, P, rng, WalkType, walkL, numpoints, StartingPoint, nburns);
+        } else {
+            uniform_sampling<DikinWalk>(randPoints, P, rng, walkL, numpoints,
+                                     StartingPoint, nburns);
+        }
+        break;
+    case john_walk:
+        if (set_L) {
+            JohnWalk WalkType(L);
+            uniform_sampling(randPoints, P, rng, WalkType, walkL, numpoints, StartingPoint, nburns);
+        } else {
+            uniform_sampling<JohnWalk>(randPoints, P, rng, walkL, numpoints,
+                                    StartingPoint, nburns);
+        }
+        break;
     case billiard:
         if(set_L) {
             BilliardWalk WalkType(L);
@@ -94,7 +125,9 @@ void sample_from_polytope(Polytope &P, int type, RNGType &rng, PointList &randPo
                                                  StartingPoint, nburns);
             }
         }
-    break;
+        break;
+    default:
+        throw Rcpp::exception("Unknown random walk!");
     }
 }
 
@@ -106,12 +139,12 @@ void sample_from_polytope(Polytope &P, int type, RNGType &rng, PointList &randPo
 //' @param n The number of points that the function is going to sample from the convex polytope.
 //' @param random_walk Optional. A list that declares the random walk and some related parameters as follows:
 //' \itemize{
-//' \item{\code{walk} }{ A string to declare the random walk: i) \code{'CDHR'} for Coordinate Directions Hit-and-Run, ii) \code{'RDHR'} for Random Directions Hit-and-Run, iii) \code{'BaW'} for Ball Walk, iv) \code{'BiW'} for Billiard walk, v) \code{'BCDHR'} boundary sampling by keeping the extreme points of CDHR or vi) \code{'BRDHR'} boundary sampling by keeping the extreme points of RDHR. The default walk is \code{'BiW'} for the uniform distribution or \code{'CDHR'} for the Gaussian distribution.}
+//' \item{\code{walk} }{ A string to declare the random walk: i) \code{'CDHR'} for Coordinate Directions Hit-and-Run, ii) \code{'RDHR'} for Random Directions Hit-and-Run, iii) \code{'BaW'} for Ball Walk, iv) \code{'BiW'} for Billiard walk, v) \code{'dikin'} for dikin walk, vi) \code{'vaidya'} for vaidya walk, vii) \code{'john'} for john walk, viii) \code{'BCDHR'} boundary sampling by keeping the extreme points of CDHR or ix) \code{'BRDHR'} boundary sampling by keeping the extreme points of RDHR. The default walk is \code{'aBiW'} for the uniform distribution or \code{'CDHR'} for the Gaussian distribution and H-polytopes and \code{'BiW'} or \code{'RDHR'} for the same distributions and V-polytopes and zonotopes.}
 //' \item{\code{walk_length} }{ The number of the steps per generated point for the random walk. The default value is \eqn{1}.}
 //' \item{\code{nburns} }{ The number of points to burn before start sampling. The default value is \eqn{1}.}
 //' \item{\code{starting_point} }{ A \eqn{d}-dimensional numerical vector that declares a starting point in the interior of the polytope for the random walk. The default choice is the center of the ball as that one computed by the function \code{inner_ball()}.}
 //' \item{\code{BaW_rad} }{ The radius for the ball walk.}
-//' \item{\code{L} }{ The maximum length of the billiard trajectory.}
+//' \item{\code{L} }{ The maximum length of the billiard trajectory or the radius for the step of dikin, vaidya or john walk.}
 //' }
 //' @param distribution Optional. A list that declares the target density and some related parameters as follows:
 //' \itemize{
@@ -121,6 +154,15 @@ void sample_from_polytope(Polytope &P, int type, RNGType &rng, PointList &randPo
 //' }
 //' @param seed Optional. A fixed seed for the number generator.
 //'
+//' @references \cite{Robert L. Smith,
+//' \dQuote{Efficient Monte Carlo Procedures for Generating Points Uniformly Distributed Over Bounded Regions,} \emph{Operations Research,} 1984.},
+//'
+//' @references \cite{B.T. Polyak, E.N. Gryazina,
+//' \dQuote{Billiard walk - a new sampling algorithm for control and optimization,} \emph{IFAC Proceedings Volumes,} 2014.},
+//'
+//' @references \cite{Y. Chen, R. Dwivedi, M. J. Wainwright and B. Yu,
+//' \dQuote{Fast MCMC Sampling Algorithms on Polytopes,} \emph{Journal of Machine Learning Research,} 2018.}
+//'
 //' @return A \eqn{d\times n} matrix that contains, column-wise, the sampled points from the convex polytope P.
 //' @examples
 //' # uniform distribution from the 3d unit cube in H-representation using ball walk
@@ -172,6 +214,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable<Rcpp::Reference> P,
 
     NT radius = 1.0, L;
     bool set_mode = false, gaussian = false, set_starting_point = false, set_L = false;
+
     random_walks walk;
 
     std::list<Point> randPoints;
@@ -230,38 +273,59 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable<Rcpp::Reference> P,
         }
     } else if (Rcpp::as<std::string>(Rcpp::as<Rcpp::List>(random_walk)["walk"]).compare(std::string("CDHR")) == 0) {
         walk = cdhr;
+    } else if (Rcpp::as<std::string>(Rcpp::as<Rcpp::List>(random_walk)["walk"]).compare(std::string("dikin")) == 0) {
+        walk = dikin_walk;
+        if (Rcpp::as<Rcpp::List>(random_walk).containsElementNamed("L")) {
+            L = Rcpp::as<NT>(Rcpp::as<Rcpp::List>(random_walk)["L"]);
+            set_L = true;
+            if (L <= 0.0) throw Rcpp::exception("L must be a postitive number!");
+        }
+    } else if (Rcpp::as<std::string>(Rcpp::as<Rcpp::List>(random_walk)["walk"]).compare(std::string("vaidya")) == 0) {
+        walk = vaidya_walk;
+        if (Rcpp::as<Rcpp::List>(random_walk).containsElementNamed("L")) {
+            L = Rcpp::as<NT>(Rcpp::as<Rcpp::List>(random_walk)["L"]);
+            set_L = true;
+            if (L <= 0.0) throw Rcpp::exception("L must be a postitive number!");
+        }
+    } else if (Rcpp::as<std::string>(Rcpp::as<Rcpp::List>(random_walk)["walk"]).compare(std::string("john")) == 0) {
+        walk = john_walk;
+        if (Rcpp::as<Rcpp::List>(random_walk).containsElementNamed("L")) {
+            L = Rcpp::as<NT>(Rcpp::as<Rcpp::List>(random_walk)["L"]);
+            set_L = true;
+            if (L <= 0.0) throw Rcpp::exception("L must be a postitive number!");
+        }
     } else if (Rcpp::as<std::string>(Rcpp::as<Rcpp::List>(random_walk)["walk"]).compare(std::string("RDHR")) == 0) {
         walk = rdhr;
     } else if (Rcpp::as<std::string>(Rcpp::as<Rcpp::List>(random_walk)["walk"]).compare(std::string("BaW")) == 0) {
         walk = ball_walk;
         if (Rcpp::as<Rcpp::List>(random_walk).containsElementNamed("BaW_rad")) {
             L = Rcpp::as<NT>(Rcpp::as<Rcpp::List>(random_walk)["BaW_rad"]);
             set_L = true;
-            if (L<=0.0) throw Rcpp::exception("BaW diameter must be a postitive number!");
+            if (L <= 0.0) throw Rcpp::exception("BaW diameter must be a postitive number!");
         }
     } else if (Rcpp::as<std::string>(Rcpp::as<Rcpp::List>(random_walk)["walk"]).compare(std::string("BiW")) == 0) {
         if (gaussian) throw Rcpp::exception("Billiard walk can be used only for uniform sampling!");
         walk = billiard;
         if (Rcpp::as<Rcpp::List>(random_walk).containsElementNamed("L")) {
             L = Rcpp::as<NT>(Rcpp::as<Rcpp::List>(random_walk)["L"]);
             set_L = true;
-            if (L<=0.0) throw Rcpp::exception("L must be a postitive number!");
+            if (L <= 0.0) throw Rcpp::exception("L must be a postitive number!");
         }
     } else if (Rcpp::as<std::string>(Rcpp::as<Rcpp::List>(random_walk)["walk"]).compare(std::string("aBiW")) == 0) {
         if (gaussian) throw Rcpp::exception("Billiard walk can be used only for uniform sampling!");
         walk = accelarated_billiard;
         if (Rcpp::as<Rcpp::List>(random_walk).containsElementNamed("L")) {
             L = Rcpp::as<NT>(Rcpp::as<Rcpp::List>(random_walk)["L"]);
             set_L = true;
-            if (L<=0.0) throw Rcpp::exception("L must be a postitive number!");
+            if (L <= 0.0) throw Rcpp::exception("L must be a postitive number!");
         }
     } else if (Rcpp::as<std::string>(Rcpp::as<Rcpp::List>(random_walk)["walk"]).compare(std::string("BRDHR")) == 0) {
         if (gaussian) throw Rcpp::exception("Gaussian sampling from the boundary is not supported!");
         walk = brdhr;
-    }else if (Rcpp::as<std::string>(Rcpp::as<Rcpp::List>(random_walk)["walk"]).compare(std::string("BCDHR")) == 0) {
+    } else if (Rcpp::as<std::string>(Rcpp::as<Rcpp::List>(random_walk)["walk"]).compare(std::string("BCDHR")) == 0) {
         if (gaussian) throw Rcpp::exception("Gaussian sampling from the boundary is not supported!");
         walk = bcdhr;
-    }else {
+    } else {
         throw Rcpp::exception("Unknown walk type!");
     }
 

diff --git a/include/cartesian_geom/point.h b/include/cartesian_geom/point.h
@@ -86,6 +86,11 @@ class point
         coeffs(i) = coord;
     }
 
+    void set_coeffs (const Coeff& coeffs2) {
+        d = coeffs2.rows();
+        coeffs = coeffs2;
+    }
+
     void set_to_origin() {
         coeffs.setZero(d);
     }

diff --git a/include/random_walks/ellipsoid_walks/README.md b/include/random_walks/ellipsoid_walks/README.md
@@ -0,0 +1,7 @@
+1) This folder contains the C++ implementations -with several modifications- of the Dikin walk, 
+   the Vaidya walk and the John walk from https://github.com/rzrsk/vaidya-walk by Raaz Dwivedi.
+
+2) The implemented random walks are presented in the paper of
+   Y. Chen, R. Dwivedi, M. J. Wainwright and B. Yu, 
+   "Fast MCMC Sampling Algorithms on Polytopes",
+   Journal of Machine Learning Research, 2018.