Exact Hamiltonian Monte Carlo for spherical Gaussian sampling (#144)

* implmentation of exact hmc walk

* implementation of boundary oracles for the hpolytope

* add fake boundary oracles for the other polytope representations

* implement c++ interface for exponential sampling

* improve boundary oracles and fix bugs

* add root finder for the degree two polynomial

* complete the root computation and the failure check

* initial implementation of the random walk

* improve R documentation and resolve gcc compile error

* improve R documentation

* initial implementation of hmc-leapfrog

* add the new walk to the R interface

* implement the C++ interface for the new walk

* fix compile and algorithmic errors

* improve R documentation

* improve comments and R interface

* initial implementation of the Exact HMC for gaussian dist.

* add the new walk to C++ and R interfaces

* implement boundary oracle

* improve boundary oracle for exact hmc for gaussian sampling

* remove test R script

* resolve PR reviews

* fix compile errors

* merge develop branch

* add burn in methods in exact hmc

* delete white space in new gaussian hmc hpp file

* resolve the reviews on the root computation

* resolve review on the header names

* improve definition of macro RVOLESTI

* add an example to sample from exponential exact HMC

* add a c++ example from the gaussian exact HMC

* remove RVOLESTI macro

* declare tol variable as static

* fix bug in variable declaration

* change tolerance variable declaration

* set upper bound for the number of relfections

* set upper bound for the number of relfections

* change the name of TOL to IN_INSIDE_BODY_TOLLERANCE

* remove VOLESTIPY macro

* resolve PR comments

* fix c++ bugs

* resolve PR's comments

R-proj/R/RcppExports.R

@@ -313,24 +313,25 @@ 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{'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 x) \code{'HMC'} for Hamiltonian Monte Carlo (logconcave) xi) \code{'ULD'} for Underdamped Langevin Dynamics using the Randomized Midpoint Method. 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} }{ 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 x) \code{'HMC'} for Hamiltonian Monte Carlo (logconcave densities) xi) \code{'ULD'} for Underdamped Langevin Dynamics using the Randomized Midpoint Method xii) \code{'ExactHMC'} for exact Hamiltonian Monte Carlo with reflections (spherical Gaussian or exponential distribution). 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 or the radius for the step of dikin, vaidya or john walk.}
-#' \item{\code{solver}} {Specify ODE solver for logconcave sampling. Options are i) leapfrog, ii) euler iii) runge-kutta iv) richardson}
-#' \item{\code{step_size} {Optionally chosen step size for logconcave sampling. Defaults to a theoretical value if not provided.}}
+#' \item{\code{solver} }{ Specify ODE solver for logconcave sampling. Options are i) leapfrog, ii) euler iii) runge-kutta iv) richardson}
+#' \item{\code{step_size }{ Optionally chosen step size for logconcave sampling. Defaults to a theoretical value if not provided.}
 #' }
 #' @param distribution Optional. A list that declares the target density and some related parameters as follows:
 #' \itemize{
-#' \item{\code{density} }{ A string: (a) \code{'uniform'} for the uniform distribution or b) \code{'gaussian'} for the multidimensional spherical distribution. The default target distribution is uniform. c) Logconcave with form proportional to exp(-f(x)) where f(x) is L-smooth and m-strongly-convex. }
-#' \item{\code{variance} }{ The variance of the multidimensional spherical gaussian. The default value is 1.}
+#' \item{\code{density} }{ A string: (a) \code{'uniform'} for the uniform distribution or b) \code{'gaussian'} for the multidimensional spherical distribution c) \code{logconcave} with form proportional to exp(-f(x)) where f(x) is L-smooth and m-strongly-convex d) \code{'exponential'} for the exponential distribution. The default target distribution is the uniform distribution.}
+#' \item{\code{variance} }{ The variance of the multidimensional spherical gaussian or the exponential distribution. The default value is 1.}
 #' \item{\code{mode} }{ A \eqn{d}-dimensional numerical vector that declares the mode of the Gaussian distribution. The default choice is the center of the as that one computed by the function \code{inner_ball()}.}
-#' \item{\code{L_}} { Smoothness constant (for logconcave). }
-#' \item{\code{m}} { Strong-convexity constant (for logconcave). }
-#' \item{\code{negative_logprob}} { Negative log-probability (for logconcave). }
-#' \item{\code{negative_logprob_gradient}} { Negative log-probability gradient (for logconcave). }
+#' \item{\code{bias} }{ The bias vector for the exponential distribution. The default vector is \eqn{c_1 = 1} and \eqn{c_i = 0} for \eqn{i \neq 1}.}
+#' \item{\code{L_} }{ Smoothness constant (for logconcave). }
+#' \item{\code{m} }{ Strong-convexity constant (for logconcave). }
+#' \item{\code{negative_logprob} }{ Negative log-probability (for logconcave). }
+#' \item{\code{negative_logprob_gradient} }{ Negative log-probability gradient (for logconcave). }
 #' }
 #' @param seed Optional. A fixed seed for the number generator.
 #'
@@ -349,6 +350,9 @@ rounding <- function(P, method = NULL, seed = NULL) {
 #' @references \cite{Shen, Ruoqi, and Yin Tat Lee,
 #' \dQuote{"The randomized midpoint method for log-concave sampling.",} \emph{Advances in Neural Information Processing Systems}, 2019.}
 #'
+#' @references \cite{Augustin Chevallier, Sylvain Pion, Frederic Cazals,
+#' \dQuote{"Hamiltonian Monte Carlo with boundary reflections, and application to polytope volume calculations,"} \emph{Research Report preprint hal-01919855}, 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

R-proj/src/copula.cpp

@@ -7,6 +7,7 @@
 
 //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program.
 
+
 #include <Rcpp.h>
 #include <RcppEigen.h>
 #include <chrono>

R-proj/src/direct_sampling.cpp

@@ -7,6 +7,7 @@
 
 //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 and 2019 program.
 
+
 #include <Rcpp.h>
 #include <RcppEigen.h>
 #include <chrono>

R-proj/src/exact_vol.cpp

@@ -5,6 +5,7 @@
 // Copyright (c) 2012-2019 Vissarion Fisikopoulos
 // Copyright (c) 2018-2019 Apostolos Chalkis
 
+
 #include <Rcpp.h>
 #include <RcppEigen.h>
 #include <chrono>

R-proj/src/get_full_dimensional_polytope.cpp

@@ -7,6 +7,7 @@
 
 //Contributed and/or modified by Alexandros Manochis, as part of Google Summer of Code 2020 program.
 
+
 #include <Rcpp.h>
 #include <RcppEigen.h>
 #include <chrono>

R-proj/src/geweke.cpp

@@ -7,6 +7,7 @@
 
 //Contributed and/or modified by Alexandros Manochis, as part of Google Summer of Code 2020 program.
 
+
 #include <Rcpp.h>
 #include <RcppEigen.h>
 #include <chrono>

R-proj/src/ode_solve.cpp

@@ -8,6 +8,7 @@
 
 //Contributed and/or modified by Marios Papachristou, as part of Google Summer of Code 2018 and 2019 program.
 
+
 #include <Rcpp.h>
 #include <RcppEigen.h>
 #include <chrono>

R-proj/src/poly_gen.cpp

@@ -7,6 +7,7 @@
 
 //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program.
 
+
 #include <Rcpp.h>
 #include <RcppEigen.h>
 #include <chrono>

R-proj/src/polytopes_modules.cpp

@@ -3,6 +3,7 @@
 // Copyright (c) 2012-2019 Vissarion Fisikopoulos
 // Copyright (c) 2018-2019 Apostolos Chalkis
 
+
 #include <Rcpp.h>
 
 

R-proj/src/psrf_multivariate.cpp

@@ -7,6 +7,7 @@
 
 //Contributed and/or modified by Alexandros Manochis, as part of Google Summer of Code 2020 program.
 
+
 #include <Rcpp.h>
 #include <RcppEigen.h>
 #include <chrono>

R-proj/src/psrf_univariate.cpp

@@ -7,6 +7,7 @@
 
 //Contributed and/or modified by Alexandros Manochis, as part of Google Summer of Code 2020 program.
 
+
 #include <Rcpp.h>
 #include <RcppEigen.h>
 #include <chrono>

R-proj/src/raftery.cpp

@@ -7,6 +7,7 @@
 
 //Contributed and/or modified by Alexandros Manochis, as part of Google Summer of Code 2020 program.
 
+
 #include <Rcpp.h>
 #include <RcppEigen.h>
 #include <chrono>

R-proj/src/rotating.cpp

@@ -7,6 +7,7 @@
 
 //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program.
 
+
 #include <Rcpp.h>
 #include <RcppEigen.h>
 #include <chrono>

R-proj/src/rounding.cpp

@@ -8,6 +8,7 @@
 //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 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>