diff --git a/R-proj/R/RcppExports.R b/R-proj/R/RcppExports.R index 2af18a4b4..f7837109a 100644 --- a/R-proj/R/RcppExports.R +++ b/R-proj/R/RcppExports.R @@ -296,24 +296,26 @@ 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 or exponential) xi) \code{'ULD'} for Underdamped Langevin Dynamics using the Randomized Midpoint Method xii) \code{'ExactHMC'} for exact Hamiltonian Monte Carlo with reflections (only 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.}} +#' \item{\code{number_of_steps }{ The number of steps for the leapfrog method, only for the exponential distribution.} #' } #' @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 uniform. } +#' \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. #' diff --git a/R-proj/src/sample_points.cpp b/R-proj/src/sample_points.cpp index 585ed6910..61b5bd331 100644 --- a/R-proj/src/sample_points.cpp +++ b/R-proj/src/sample_points.cpp @@ -2,8 +2,8 @@ // VolEsti (volume computation and sampling library) -// Copyright (c) 2012-2020 Vissarion Fisikopoulos -// Copyright (c) 2018-2020 Apostolos Chalkis +// Copyright (c) 2012-2021 Vissarion Fisikopoulos +// Copyright (c) 2018-2021 Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 and 2019 program. @@ -36,6 +36,8 @@ enum random_walks { brdhr, bcdhr, hmc, + exponential_hmc, + exponential_hmc_leapfrog, uld }; @@ -50,10 +52,10 @@ template < > void sample_from_polytope(Polytope &P, int type, RNGType &rng, PointList &randPoints, unsigned int const& walkL, unsigned int const& numpoints, - bool const& gaussian, NT const& a, NT const& L, + bool const& gaussian, NT const& a, NT const& L, Point const& c, Point const& StartingPoint, unsigned int const& nburns, - bool const& set_L, random_walks walk, - NegativeGradientFunctor *F=NULL, NegativeLogprobFunctor *f=NULL, + bool const& set_L, bool const& set_nsteps, unsigned int const& nsteps, NT const& eta, + random_walks walk, NegativeGradientFunctor *F=NULL, NegativeLogprobFunctor *f=NULL, ode_solvers solver_type = no_solver) { switch (walk) @@ -129,6 +131,24 @@ void sample_from_polytope(Polytope &P, int type, RNGType &rng, PointList &randPo StartingPoint, nburns); } break; + case exponential_hmc: + if (set_L) { + HMCExponentialWalk WalkType(L); + exponential_sampling(randPoints, P, rng, WalkType, walkL, numpoints, c, a, StartingPoint, nburns); + } else { + exponential_sampling(randPoints, P, rng, walkL, numpoints, c, a, + StartingPoint, nburns); + } + break; + case exponential_hmc_leapfrog: + if (set_nsteps) { + HMCLeafrogExponential WalkType(nsteps); + exponential_sampling(randPoints, P, rng, WalkType, walkL, numpoints, c, a, eta, StartingPoint, nburns); + } else { + exponential_sampling(randPoints, P, rng, walkL, numpoints, c, a, eta, + StartingPoint, nburns); + } + break; case ball_walk: if (set_L) { if (gaussian) { @@ -222,24 +242,26 @@ 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{'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 or exponential densities) xi) \code{'ULD'} for Underdamped Langevin Dynamics using the Randomized Midpoint Method xii) \code{'ExactHMC'} for exact Hamiltonian Monte Carlo with reflections (only 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. It is required for the exponential distribution.} +//' \item{\code{number_of_steps }{ Optionally the number of steps for the leapfrog method, only for the exponential distribution.} //' } //' @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. //' @@ -296,7 +318,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, typedef Eigen::Matrix MT; unsigned int type = Rcpp::as(P).field("type"), dim = Rcpp::as(P).field("dimension"), - walkL = 1, numpoints, nburns = 0; + walkL = 1, numpoints, nburns = 0, nsteps = 0; RcppFunctor::GradientFunctor *F = NULL; RcppFunctor::FunctionFunctor *f = NULL; @@ -308,13 +330,16 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, rng.set_seed(seed_rcpp); } + Point c(dim); + NT radius = 1.0, L; - bool set_mode = false, gaussian = false, logconcave = false, - set_starting_point = false, set_L = false; + bool set_mode = false, gaussian = false, logconcave = false, exponential = false, set_eta = false, + set_starting_point = false, set_L = false, set_nsteps = false; random_walks walk; ode_solvers solver; // Used only for logconcave sampling + NT eta; std::list randPoints; std::pair InnerBall; Point mode(dim); @@ -330,6 +355,10 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, } else if ( Rcpp::as(Rcpp::as(distribution)["density"]).compare(std::string("gaussian")) == 0) { gaussian = true; + } else if ( + Rcpp::as(Rcpp::as(distribution)["density"]).compare(std::string("exponential")) == 0) { + walk = exponential_hmc_leapfrog; + exponential = true; } else if ( Rcpp::as(Rcpp::as(distribution)["density"]).compare(std::string("logconcave")) == 0) { logconcave = true; @@ -351,13 +380,34 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, NT a = 0.5; if (Rcpp::as(distribution).containsElementNamed("variance")) { a = 1.0 / (2.0 * Rcpp::as(Rcpp::as(distribution)["variance"])); - if (!gaussian) { - Rcpp::warning("The variance can be set only for Gaussian sampling!"); + if (exponential) a = Rcpp::as(Rcpp::as(distribution)["variance"]); + if (!gaussian && !exponential) { + Rcpp::warning("The variance can be set only for Gaussian and exponential sampling!"); } else if (a <= 0.0) { throw Rcpp::exception("The variance has to be positive!"); } } + if (Rcpp::as(distribution).containsElementNamed("bias")) { + c = Point(Rcpp::as(Rcpp::as(distribution)["bias"])); + } + + if (Rcpp::as(random_walk).containsElementNamed("number_of_steps")) { + nsteps = NT(Rcpp::as(Rcpp::as(random_walk)["number_of_steps"])); + if (nsteps <= 0) { + throw Rcpp::exception("Number of steps must be a positive integer"); + } + set_nsteps = true; + } + + if (Rcpp::as(random_walk).containsElementNamed("step_size")) { + eta = NT(Rcpp::as(Rcpp::as(random_walk)["step_size"])); + if (eta <= NT(0)) { + throw Rcpp::exception("Step size must be positive"); + } + set_eta = true; + } + if (Rcpp::as(distribution).containsElementNamed("negative_logprob") && Rcpp::as(distribution).containsElementNamed("negative_logprob_gradient")) { @@ -383,7 +433,6 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, throw Rcpp::exception("The strong-convexity constant is absent"); } - if (Rcpp::as(random_walk).containsElementNamed("step_size")) { eta = NT(Rcpp::as(Rcpp::as(random_walk)["step_size"])); if (eta <= NT(0)) { @@ -414,9 +463,14 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, } - if (!random_walk.isNotNull() || !Rcpp::as(random_walk).containsElementNamed("walk")) { - if (gaussian) { + if (exponential) { + if (type !=1) { + throw Rcpp::exception("Exponential sampling is supported only for H-polytopes"); + } + if (!set_eta) throw Rcpp::exception("You have to set the step size of the leapfrog method."); + walk = exponential_hmc_leapfrog; + } else if (gaussian) { if (type == 1) { walk = cdhr; } else { @@ -478,14 +532,22 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, if (L <= 0.0) throw Rcpp::exception("L must be a postitive number!"); } } else if (Rcpp::as(Rcpp::as(random_walk)["walk"]).compare(std::string("BRDHR")) == 0) { - if (gaussian) throw Rcpp::exception("Gaussian sampling from the boundary is not supported!"); + if (gaussian || exponential) throw Rcpp::exception("Gaussian sampling from the boundary is not supported!"); walk = brdhr; } else if (Rcpp::as(Rcpp::as(random_walk)["walk"]).compare(std::string("BCDHR")) == 0) { if (gaussian) throw Rcpp::exception("Gaussian sampling from the boundary is not supported!"); walk = bcdhr; + } else if(Rcpp::as(Rcpp::as(random_walk)["walk"]).compare(std::string("ExactHMC")) == 0) { + if (!exponential) throw Rcpp::exception("Exact HMC is supported only for exponential sampling."); + walk = exponential_hmc; } else if (Rcpp::as(Rcpp::as(random_walk)["walk"]).compare(std::string("HMC")) == 0) { - if (!logconcave) throw Rcpp::exception("HMC is not supported for non first-order sampling"); - walk = hmc; + if (!logconcave && !exponential) throw Rcpp::exception("HMC is not supported for non first-order sampling"); + if (exponential) { + if (!set_eta) throw Rcpp::exception("You have to set the step size"); + walk = exponential_hmc_leapfrog; + } else { + walk = hmc; + } } else if (Rcpp::as(Rcpp::as(random_walk)["walk"]).compare(std::string("ULD")) == 0) { if (!logconcave) throw Rcpp::exception("ULD is not supported for non first-order sampling"); walk = uld; @@ -537,8 +599,8 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, StartingPoint = StartingPoint - mode; HP.shift(mode.getCoefficients()); } - sample_from_polytope(HP, type, rng, randPoints, walkL, numpoints, gaussian, a, L, - StartingPoint, nburns, set_L, walk, F, f, solver); + sample_from_polytope(HP, type, rng, randPoints, walkL, numpoints, gaussian, a, L, c, + StartingPoint, nburns, set_L, set_nsteps, nsteps, eta, walk, F, f, solver); break; } case 2: { @@ -558,8 +620,8 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, StartingPoint = StartingPoint - mode; VP.shift(mode.getCoefficients()); } - sample_from_polytope(VP, type, rng, randPoints, walkL, numpoints, gaussian, a, L, - StartingPoint, nburns, set_L, walk, F, f, solver); + sample_from_polytope(VP, type, rng, randPoints, walkL, numpoints, gaussian, a, L, c, + StartingPoint, nburns, set_L, set_nsteps, nsteps, eta, walk, F, f, solver); break; } case 3: { @@ -579,8 +641,8 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, StartingPoint = StartingPoint - mode; ZP.shift(mode.getCoefficients()); } - sample_from_polytope(ZP, type, rng, randPoints, walkL, numpoints, gaussian, a, L, - StartingPoint, nburns, set_L, walk, F, f, solver); + sample_from_polytope(ZP, type, rng, randPoints, walkL, numpoints, gaussian, a, L, c, + StartingPoint, nburns, set_L, set_nsteps, nsteps, eta, walk, F, f, solver); break; } case 4: { @@ -602,8 +664,8 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, StartingPoint = StartingPoint - mode; VPcVP.shift(mode.getCoefficients()); } - sample_from_polytope(VPcVP, type, rng, randPoints, walkL, numpoints, gaussian, a, L, - StartingPoint, nburns, set_L, walk, F, f, solver); + sample_from_polytope(VPcVP, type, rng, randPoints, walkL, numpoints, gaussian, a, L, c, + StartingPoint, nburns, set_L, set_nsteps, nsteps, eta, walk, F, f, solver); break; } } diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index 3cf5f41f2..e978e5e0d 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -16,6 +16,7 @@ #include #include #include "preprocess/max_inscribed_ball.hpp" +#include "root_finders/quadratic_polynomial_solvers.hpp" #ifndef VOLESTIPY #include "lp_oracles/solve_lp.h" #endif @@ -638,6 +639,124 @@ class HPolytope { } + //------------------------------oracles for exponential sampling---------------////// + + // compute intersection points of a ray starting from r and pointing to v + // with polytope discribed by A and b + std::pair quadratic_positive_intersect(Point const& r, + Point const& v, + VT& Ac, + NT const& T, + VT& Ar, + VT& Av, + int& facet_prev) const + { + NT lamda = 0, lamda2 =0, lamda1 =0, alpha; + NT min_plus = std::numeric_limits::max(); + NT max_minus = std::numeric_limits::lowest(); + VT sum_nom; + bool real; + int m = num_of_hyperplanes(), facet = -1; + + Ar.noalias() = A * r.getCoefficients(); + sum_nom = Ar - b; + Av.noalias() = A * v.getCoefficients();; + + NT* Av_data = Av.data(); + NT* sum_nom_data = sum_nom.data(); + NT* Ac_data = Ac.data(); + + for (int i = 0; i < m; i++) { + alpha = -((*Ac_data) / (2.0 * T)); + solve_quadratic_polynomial(alpha, (*Av_data), (*sum_nom_data), lamda1, lamda2, real); + if (real) { + lamda = pick_intersection_time(lamda1, lamda2, i, facet_prev); + if (lamda < min_plus && lamda > 0) { + min_plus = lamda; + facet = i; + } + } + Av_data++; + sum_nom_data++; + Ac_data++; + } + facet_prev = facet; + return std::make_pair(min_plus, facet); + } + + std::pair quadratic_positive_intersect(Point const& r, + Point const& v, + VT& Ac, + NT const& T, + VT& Ar, + VT& Av, + NT const& lambda_prev, + int& facet_prev) const + { + + NT lamda = 0, lamda2 =0, lamda1 =0, alpha; + NT min_plus = std::numeric_limits::max(); + NT max_minus = std::numeric_limits::lowest(); + VT sum_nom; + NT mult; + unsigned int j; + int m = num_of_hyperplanes(), facet = -1; + bool real; + + Ar.noalias() += ((lambda_prev * lambda_prev) / (-2.0*T)) * Ac + lambda_prev * Av; + sum_nom = Ar - b; + Av.noalias() = A * v.getCoefficients(); + + NT* sum_nom_data = sum_nom.data(); + NT* Av_data = Av.data(); + NT* Ac_data = Ac.data(); + + for (int i = 0; i < m; i++) { + alpha = -((*Ac_data) / (2.0 * T)); + solve_quadratic_polynomial(alpha, (*Av_data), (*sum_nom_data), lamda1, lamda2, real); + if (real) { + lamda = pick_intersection_time(lamda1, lamda2, i, facet_prev); + if (lamda < min_plus && lamda > 0) { + min_plus = lamda; + facet = i; + } + } + Av_data++; + sum_nom_data++; + Ac_data++; + } + facet_prev = facet; + return std::make_pair(min_plus, facet); + } + + NT pick_intersection_time(NT lamda1, NT lamda2, int current_facet, int previous_facet) const + { + NT lamda; + if (lamda1 * lamda2 < NT(0)) { + if (previous_facet == current_facet) { + if (std::max(lamda1, lamda2) < 1e-10) { + lamda = std::min(lamda1, lamda2); + } else { + lamda = std::max(lamda1, lamda2); + } + } else { + lamda = std::max(lamda1, lamda2); + } + } else { + if (previous_facet == current_facet) { + if (std::min(lamda1, lamda2) >= NT(0) && std::min(lamda1, lamda2) < 1e-10) { + lamda = std::max(lamda1, lamda2); + } else { + lamda = std::min(lamda1, lamda2); + } + } else { + lamda = std::min(lamda1, lamda2); + } + } + return lamda; + } + + // Apply linear transformation, of square matrix T^{-1}, in H-polytope P:= Ax<=b void linear_transformIt(MT const& T) { diff --git a/include/convex_bodies/vpolyintersectvpoly.h b/include/convex_bodies/vpolyintersectvpoly.h index 3924fbe54..19b5f8d60 100644 --- a/include/convex_bodies/vpolyintersectvpoly.h +++ b/include/convex_bodies/vpolyintersectvpoly.h @@ -50,7 +50,7 @@ class IntersectionOfVpoly { return _inner_ball; } - int is_in(const Point &p) const { + int is_in(const Point &p, NT tol=NT(0)) const { if(P1.is_in(p)==-1) return P2.is_in(p); return 0; @@ -322,6 +322,34 @@ class IntersectionOfVpoly { } + //------------------------------oracles for exponential sampling---------------////// + + // compute intersection points of a ray starting from r and pointing to v + // with polytope discribed by A and b + std::pair quadratic_positive_intersect(Point const& r, + Point const& v, + VT const& Ac, + NT const& T, + VT& Ar, + VT& Av, + int& facet_prev) const + { + return std::make_pair(0, 0); + } + + std::pair quadratic_positive_intersect(Point const& r, + Point const& v, + VT const& Ac, + NT const& T, + VT& Ar, + VT& Av, + NT const& lambda_prev, + int& facet_prev) const + { + return std::make_pair(0, 0); + } + + // shift polytope by a point c void shift(const VT &c) { P1.shift(c); diff --git a/include/convex_bodies/vpolytope.h b/include/convex_bodies/vpolytope.h index a661be269..ce6b1ffdf 100644 --- a/include/convex_bodies/vpolytope.h +++ b/include/convex_bodies/vpolytope.h @@ -493,6 +493,34 @@ class VPolytope { } + //------------------------------oracles for exponential sampling---------------////// + + // compute intersection points of a ray starting from r and pointing to v + // with polytope discribed by A and b + std::pair quadratic_positive_intersect(Point const& r, + Point const& v, + VT const& Ac, + NT const& T, + VT& Ar, + VT& Av, + int& facet_prev) const + { + return std::make_pair(0, 0); + } + + std::pair quadratic_positive_intersect(Point const& r, + Point const& v, + VT const& Ac, + NT const& T, + VT& Ar, + VT& Av, + NT const& lambda_prev, + int& facet_prev) const + { + return std::make_pair(0, 0); + } + + // shift polytope by a point c void shift(const VT &c) { MT V2 = V.transpose().colwise() - c; diff --git a/include/convex_bodies/zpolytope.h b/include/convex_bodies/zpolytope.h index 4d70b8a46..2d1f781cc 100644 --- a/include/convex_bodies/zpolytope.h +++ b/include/convex_bodies/zpolytope.h @@ -336,7 +336,7 @@ class Zonotope { // check if point p belongs to the convex hull of V-Polytope P - int is_in(Point const& p) const + int is_in(Point const& p, NT tol=NT(0)) const { if(memLP_Zonotope(V, p, row_mem, colno_mem)) { @@ -477,6 +477,35 @@ class Zonotope { { return line_intersect_coord(r, rand_coord, lamdas); } + + + //------------------------------oracles for exponential sampling---------------////// + + // compute intersection points of a ray starting from r and pointing to v + // with polytope discribed by A and b + std::pair quadratic_positive_intersect(Point const& r, + Point const& v, + VT const& Ac, + NT const& T, + VT& Ar, + VT& Av, + int& facet_prev) const + { + return std::make_pair(0, 0); + } + + std::pair quadratic_positive_intersect(Point const& r, + Point const& v, + VT const& Ac, + NT const& T, + VT& Ar, + VT& Av, + NT const& lambda_prev, + int& facet_prev) const + { + return std::make_pair(0, 0); + } + // shift polytope by a point c diff --git a/include/random_walks/exponential_hamiltonian_monte_carlo_walk.hpp b/include/random_walks/exponential_hamiltonian_monte_carlo_walk.hpp new file mode 100644 index 000000000..81407206b --- /dev/null +++ b/include/random_walks/exponential_hamiltonian_monte_carlo_walk.hpp @@ -0,0 +1,205 @@ +// VolEsti (volume computation and sampling library) +// Copyright (c) 2021 Vissarion Fisikopoulos +// Copyright (c) 2021 Apostolos Chalkis + + +// Licensed under GNU LGPL.3, see LICENCE file + +#ifndef RANDOM_WALKS_EXPONENTIAL_HMC_WALK_HPP +#define RANDOM_WALKS_EXPONENTIAL_HMC_WALK_HPP + +#include "sampling/sphere.hpp" + + +// Exact HMC for sampling from the Exponential distribution restricted to a convex polytope + +struct HMCExponentialWalk +{ + HMCExponentialWalk(double L) + : param(L, true) + {} + + HMCExponentialWalk() + : param(0, false) + {} + + struct parameters + { + parameters(double L, bool set) + : m_L(L), set_L(set) + {} + double m_L; + bool set_L; + }; + + parameters param; + + +template +< + typename Polytope, + typename RandomNumberGenerator +> +struct Walk +{ + typedef typename Polytope::PointType Point; + typedef typename Point::FT NT; + typedef typename Polytope::MT MT; + typedef typename Polytope::VT VT; + + template + Walk(GenericPolytope const& P, Point const& p, Point const& c, NT const& T, RandomNumberGenerator &rng) + { + _Len = compute_diameter + ::template compute(P); + _Ac = P.get_mat() * c.getCoefficients(); + _Temp = T; + _c = c; + initialize(P, p, rng); + } + + template + Walk(GenericPolytope const& P, Point const& p, Point const& c, NT const& T, RandomNumberGenerator &rng, + parameters const& params) + { + _Len = params.set_L ? params.m_L + : compute_diameter + ::template compute(P); + _Ac = P.get_mat() * c.getCoefficients(); + _Temp = T; + _c = c; + initialize(P, p, rng); + } + + template + < + typename GenericPolytope + > + inline bool apply(GenericPolytope const& P, + Point& p, // a point to start + unsigned int const& walk_length, + RandomNumberGenerator &rng) + { + unsigned int n = P.dimension(); + NT T; + int failures = 0, it; + Point p0; + + for (auto j=0u; j::apply(n, rng, false); + p0 = _p; + it = 0; + while (it < 200*n) + { + auto pbpair = P.quadratic_positive_intersect(_p, _v, _Ac, _Temp, _lambdas, + _Av, _lambda_prev, _facet_prev); + if (T <= pbpair.first || pbpair.second < 0) { + _p += ((T * T) / (-2.0*_Temp)) *_c + (T * _v); + _lambda_prev = T; + break; + } + _lambda_prev = pbpair.first; + _p += ((_lambda_prev * _lambda_prev) / (-2.0*_Temp)) *_c + (_lambda_prev * _v); + T -= _lambda_prev; + _v += (-_lambda_prev/_Temp) * _c; + P.compute_reflection(_v, _p, pbpair.second); + it++; + } + + } while (P.is_in(_p, _tol) == 0); + if (it == 200*n){ + _p = p0; + } + } + p = _p; + return true; + } + + inline void update_delta(NT L) + { + _Len = L; + } + +private : + + template + < + typename GenericPolytope + > + inline void initialize(GenericPolytope const& P, + Point const& p, + RandomNumberGenerator &rng) + { + unsigned int n = P.dimension(); + NT T; + _lambdas.setZero(P.num_of_hyperplanes()); + _Av.setZero(P.num_of_hyperplanes()); + _p = p; + _v = GetDirection::apply(n, rng, false); + + do { + T = -std::log(rng.sample_urdist()) * _Len; + Point p0 = _p; + int it = 0; + + std::pair pbpair + = P.quadratic_positive_intersect(_p, _v, _Ac, _Temp, _lambdas, _Av, _facet_prev); + if (T <= pbpair.first || pbpair.second < 0) { + _p += ((T * T) / (-2.0*_Temp)) *_c + (T * _v); + _lambda_prev = T; + return; + } + _lambda_prev = pbpair.first; + _p += ((_lambda_prev * _lambda_prev) / (-2.0*_Temp)) *_c + (_lambda_prev * _v); + _v += (-_lambda_prev/_Temp) * _c; + T -= _lambda_prev; + P.compute_reflection(_v, _p, pbpair.second); + + while (it <= 100*n) + { + std::pair pbpair + = P.quadratic_positive_intersect(_p, _v, _Ac, _Temp, _lambdas, _Av, _lambda_prev, _facet_prev); + if (T <= pbpair.first || pbpair.second < 0) { + _p += ((T * T) / (-2.0*_Temp)) *_c + (T * _v); + _lambda_prev = T; + break; + }else if (it == 100*n) { + _lambda_prev = rng.sample_urdist() * pbpair.first; + _p += ((_lambda_prev * _lambda_prev) / (-2.0*_Temp)) *_c + (_lambda_prev * _v); + break; + } + _lambda_prev = pbpair.first; + _p += ((_lambda_prev * _lambda_prev) / (-2.0*_Temp)) *_c + (_lambda_prev * _v); + _v += (-_lambda_prev/_Temp) * _c; + T -= _lambda_prev; + P.compute_reflection(_v, _p, pbpair.second); + it++; + } + } while (P.is_in(_p, _tol) == 0); + } + + NT _Len; + VT _Ac; + Point _p; + Point _v; + Point _c; + NT _Temp; + const double _tol = 1e-10; + NT _lambda_prev; + int _facet_prev; + typename Point::Coeff _lambdas; + typename Point::Coeff _Av; +}; + +}; + + +#endif // RANDOM_WALKS_EXPONENTIAL_HMC_WALK_HPP + diff --git a/include/random_walks/exponential_hmc_leapfrog.hpp b/include/random_walks/exponential_hmc_leapfrog.hpp new file mode 100644 index 000000000..db1ad6f71 --- /dev/null +++ b/include/random_walks/exponential_hmc_leapfrog.hpp @@ -0,0 +1,248 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2021 Vissarion Fisikopoulos +// Copyright (c) 2021 Apostolos Chalkis + +// Licensed under GNU LGPL.3, see LICENCE file + +#ifndef RANDOM_WALKS_HAMILTONIAN_HMC_LEAPFROG_HPP +#define RANDOM_WALKS_HAMILTONIAN_HMC_LEAPFROG_HPP + +#include "sampling/sphere.hpp" + + +// Hamiltonian Monte Carlo with leapfrog for sampling from the exponential distribution + +struct HMCLeafrogExponential +{ + HMCLeafrogExponential(unsigned int steps) + : param(steps, true) + {} + + HMCLeafrogExponential() + : param(0, false) + {} + + struct parameters + { + parameters(unsigned int steps, bool set) + : nsteps(steps), set_steps(set) + {} + unsigned int nsteps; + bool set_steps; + }; + + struct update_parameters + { + update_parameters() + : facet_prev(0), hit_ball(false), inner_vi_ak(0.0), ball_inner_norm(0.0) + {} + int facet_prev; + bool hit_ball; + double inner_vi_ak; + double ball_inner_norm; + }; + + parameters param; + + + template + < + typename Polytope, + typename RandomNumberGenerator + > + struct Walk + { + typedef typename Polytope::PointType Point; + typedef typename Polytope::MT MT; + typedef typename Point::FT NT; + + template + Walk(GenericPolytope const& P, Point const& p, Point const& c, NT const& T, + NT const& eta, RandomNumberGenerator &rng) + { + _update_parameters = update_parameters(); + NT L = compute_diameter + ::template compute(P); + _steps = int(L / eta); + _eta = eta; + _Temp = T; + _c = c; + _AA.noalias() = P.get_mat() * P.get_mat().transpose(); + initialize(P, p, rng); + } + + template + Walk(GenericPolytope const& P, Point const& p, Point const& c, NT const& T, + NT const& eta, RandomNumberGenerator &rng, parameters const& params) + { + _update_parameters = update_parameters(); + _eta = eta; + _Temp = T; + if (params.set_steps) { + _steps = params.nsteps; + } else { + NT L = compute_diameter + ::template compute(P); + _steps = int(L / _eta); + } + _AA.noalias() = P.get_mat() * P.get_mat().transpose(); + _c = c; + initialize(P, p, rng); + } + + template + < + typename GenericPolytope + > + inline void apply(GenericPolytope const& P, + Point &p, // a point to start + unsigned int const& walk_length, + RandomNumberGenerator &rng) + { + unsigned int n = P.dimension(); + NT T; + const NT dl = 0.995; + int it; + + for (auto j=0u; j::apply(n, rng, false); + _p0 = p; + _H = _c.dot(_p) + 0.5 * _v.dot(_v); + for (auto k=0u; k<_steps; ++k) + { + T = _eta; + _v += (_eta / (-2.0 * _Temp)) * _c; + + it = 0; + std::pair pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, + _lambda_prev, _update_parameters); + if (T <= pbpair.first) { + _p += (T * _v); + _v += (_eta / (-2.0 * _Temp)) * _c; + _lambda_prev = T; + continue; + } + + _lambda_prev = dl * pbpair.first; + _p += (_lambda_prev * _v); + T -= _lambda_prev; + P.compute_reflection(_v, _p, _update_parameters); + it++; + + while (it < 500*n) + { + std::pair pbpair + = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, + _AA, _update_parameters); + if (T <= pbpair.first) { + _p += (T * _v); + _v += (_eta / (-2.0 * _Temp)) * _c; + _lambda_prev = T; + break; + } + _lambda_prev = dl * pbpair.first; + _p += (_lambda_prev * _v); + T -= _lambda_prev; + P.compute_reflection(_v, _p, _update_parameters); + it++; + } + } + _Htilde = _c.dot(_p) + 0.5 * _v.dot(_v); + + NT log_prob = _H - _Htilde < 0 ? _H - _Htilde : 0; + NT u_logprob = log(rng.sample_urdist()); + + if (u_logprob > log_prob) { + _p = _p0; + _Av = _Av_prev; + _lambda_prev = _lambda_prev_0; + _lambdas = _lambdas_prev; + } else { + _Av_prev = _Av; + _lambda_prev_0 = _lambda_prev; + _lambdas_prev = _lambdas; + } + } + p = _p; + } + + private : + + template + < + typename GenericPolytope + > + inline void initialize(GenericPolytope const& P, + Point const& p, + RandomNumberGenerator &rng) + { + unsigned int n = P.dimension(); + const NT dl = 0.995; + NT T = _eta; + _lambdas.setZero(P.num_of_hyperplanes()); + _Av.setZero(P.num_of_hyperplanes()); + _p = p; + _v = GetDirection::apply(n, rng, false); + _v += (_eta / (-2.0 * _Temp)) * _c; + + int it = 0; + + std::pair pbpair + = P.line_first_positive_intersect(_p, _v, _lambdas, _Av, _update_parameters); + if (T <= pbpair.first) { + _p += (T * _v); + _lambda_prev = T; + + _Av_prev = _Av; + _lambda_prev_0 = _lambda_prev; + _lambdas_prev = _lambdas; + return; + } + _lambda_prev = dl * pbpair.first; + _p += (_lambda_prev * _v); + T -= _lambda_prev; + P.compute_reflection(_v, _p, _update_parameters); + + while (it <= 500*n) + { + std::pair pbpair + = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _AA, _update_parameters); + if (T <= pbpair.first) { + _p += (T * _v); + _lambda_prev = T; + break; + } else if (it == 500*n) { + _lambda_prev = rng.sample_urdist() * pbpair.first; + _p += (_lambda_prev * _v); + break; + } + _lambda_prev = dl * pbpair.first; + _p += (_lambda_prev * _v); + T -= _lambda_prev; + P.compute_reflection(_v, _p, _update_parameters); + it++; + } + _Av_prev = _Av; + _lambda_prev_0 = _lambda_prev; + _lambdas_prev = _lambdas; + } + + + Point _p, _v, _c, _p0; + NT _lambda_prev, _lambda_prev_0, _eta, _H, _Htilde, _Temp; + unsigned int _steps; + MT _AA; + update_parameters _update_parameters; + typename Point::Coeff _lambdas, _lambdas_prev; + typename Point::Coeff _Av, _Av_prev; + }; + +}; + + +#endif + + + diff --git a/include/random_walks/random_walks.hpp b/include/random_walks/random_walks.hpp index 4c3f5082c..2bb8768e5 100644 --- a/include/random_walks/random_walks.hpp +++ b/include/random_walks/random_walks.hpp @@ -20,6 +20,8 @@ #include "random_walks/uniform_john_walk.hpp" #include "random_walks/uniform_vaidya_walk.hpp" #include "random_walks/uniform_accelerated_billiard_walk.hpp" +#include "random_walks/exponential_hamiltonian_monte_carlo_walk.hpp" +#include "random_walks/exponential_hmc_leapfrog.hpp" #ifndef VOLESTIPY #include "random_walks/hamiltonian_monte_carlo_walk.hpp" #include "random_walks/langevin_walk.hpp" diff --git a/include/root_finders/quadratic_polynomial_solvers.hpp b/include/root_finders/quadratic_polynomial_solvers.hpp new file mode 100644 index 000000000..1b624faf8 --- /dev/null +++ b/include/root_finders/quadratic_polynomial_solvers.hpp @@ -0,0 +1,77 @@ +// VolEsti (volume computation and sampling library) +// Copyright (c) 2021 Vissarion Fisikopoulos +// Copyright (c) 2021 Apostolos Chalkis + + +// Licensed under GNU LGPL.3, see LICENCE file + + +#ifndef QUADRATIC_POLYNOMIAL_SOLVERS_H +#define QUADRATIC_POLYNOMIAL_SOLVERS_H + +// The function returns the sigh of the input number +template int sgn(T val) +{ + return (T(0) < val) - (val < T(0)); +} + +// The function compute the roots of a quadratic polynomial equation +template +void solve_quadratic_polynomial(NT a, NT b, NT c, NT &x1, NT &x2, bool &real) { + + NT Delta = b * b - 4.0 * a * c; + if (Delta < NT(0)) { + real = false; + return; + } + real = true; + if (b >= NT(0)){ + x1 = (- b - std::sqrt(Delta)) / (2.0 * a); + x2 = (2.0 * c) / (- b - std::sqrt(Delta)); + } else { + x1 = (2.0 * c) / (- b + std::sqrt(Delta)); + x2 = (- b + std::sqrt(Delta)) / (2.0 * a); + } +} + +// The function compute the roots of a quadratic polynomial equation +// using the algorithm in Revisiting the stability of computing the roots of a quadratic polynomial +// by Nicola Mastronardi and Paul Michel Van Dooren, Electronic transactions on numerical analysis, 2014 +template +void solve_qudratic_polynomial_stable(NT a, NT b, NT c, NT &x1, NT &x2, bool &real) +{ + real = true; + + if (c == NT(0)) { + x1 = NT(0); + x2 = -b/a; + return; + } + b = b / a; + c = c / a; + + NT alpha = sgn(b) * std::sqrt(std::abs(c)); + NT e = NT(sgn(c)); + NT beta = b / (NT(2)*alpha); + + if (e > NT(0) && beta < NT(1)) + { + real = false; + return; + } + + if (e < NT(0)) + { + x1 = beta + std::sqrt(beta*beta + NT(1)); + x2 = -NT(1)/x1; + } else + { + x1 = beta + std::sqrt((beta + NT(1)) * (beta - NT(1))); + x2 = NT(1) / x1; + } + x1 = -x1 * alpha; + x2 = -x2 * alpha; +} + +#endif + diff --git a/include/sampling/random_point_generators.hpp b/include/sampling/random_point_generators.hpp index 2bb711de6..a025bce15 100644 --- a/include/sampling/random_point_generators.hpp +++ b/include/sampling/random_point_generators.hpp @@ -204,5 +204,137 @@ struct LogconcaveRandomPointGenerator }; +template +< + typename Walk +> +struct ExponentialRandomPointGenerator +{ + template + < + typename Polytope, + typename Point, + typename NT, + typename PointList, + typename WalkPolicy, + typename RandomNumberGenerator + > + static void apply(Polytope const& P, + Point &p, // a point to start + Point const& c, // bias function + NT const& T, // temperature/variance + unsigned int const& rnum, + unsigned int const& walk_length, + PointList &randPoints, + WalkPolicy &policy, + RandomNumberGenerator &rng) + { + Walk walk(P, p, c, T, rng); + bool success; + for (unsigned int i=0; i + static void apply(Polytope const& P, + Point &p, // a point to start + Point const& c, // bias function + NT const& T, // temperature/variance + unsigned int const& rnum, + unsigned int const& walk_length, + PointList &randPoints, + WalkPolicy &policy, + RandomNumberGenerator &rng, + Parameters const& parameters) + { + Walk walk(P, p, c, T, rng, parameters); + bool success; + + for (unsigned int i=0; i + static void apply(Polytope const& P, + Point &p, // a point to start + Point const& c, // bias function + NT const& T, // temperature/variance + NT const& eta, + unsigned int const& rnum, + unsigned int const& walk_length, + PointList &randPoints, + WalkPolicy &policy, + RandomNumberGenerator &rng) + { + Walk walk(P, p, c, T, eta, rng); + for (unsigned int i=0; i + static void apply(Polytope const& P, + Point &p, // a point to start + Point const& c, // bias function + NT const& T, // temperature/variance + NT const& eta, + unsigned int const& rnum, + unsigned int const& walk_length, + PointList &randPoints, + WalkPolicy &policy, + RandomNumberGenerator &rng, + Parameters const& parameters) + { + Walk walk(P, p, c, T, eta, rng, parameters); + + for (unsigned int i=0; i RandomPointGenerator; - RandomPointGenerator::apply(P, p, nburns, walk_len, randPoints, - push_back_policy, rng); - randPoints.clear(); + if (nburns > 0) { + RandomPointGenerator::apply(P, p, nburns, walk_len, randPoints, + push_back_policy, rng); + randPoints.clear(); + } RandomPointGenerator::apply(P, p, rnum, walk_len, randPoints, push_back_policy, rng); @@ -85,10 +87,11 @@ void uniform_sampling(PointList &randPoints, typedef RandomPointGenerator RandomPointGenerator; Point p = starting_point; - - RandomPointGenerator::apply(P, p, nburns, walk_len, randPoints, - push_back_policy, rng, WalkType.param); - randPoints.clear(); + if (nburns > 0) { + RandomPointGenerator::apply(P, p, nburns, walk_len, randPoints, + push_back_policy, rng, WalkType.param); + randPoints.clear(); + } RandomPointGenerator::apply(P, p, rnum, walk_len, randPoints, push_back_policy, rng, WalkType.param); } @@ -122,10 +125,11 @@ void uniform_sampling_boundary(PointList &randPoints, Point p = starting_point; typedef BoundaryRandomPointGenerator BoundaryRandomPointGenerator; - BoundaryRandomPointGenerator::apply(P, p, nburns, walk_len, - randPoints, push_back_policy, rng); - - randPoints.clear(); + if (nburns > 0) { + BoundaryRandomPointGenerator::apply(P, p, nburns, walk_len, + randPoints, push_back_policy, rng); + randPoints.clear(); + } unsigned int n = rnum / 2; BoundaryRandomPointGenerator::apply(P, p, rnum / 2, walk_len, randPoints, push_back_policy, rng); @@ -164,9 +168,11 @@ void gaussian_sampling(PointList &randPoints, Point p = starting_point; typedef GaussianRandomPointGenerator RandomPointGenerator; - RandomPointGenerator::apply(P, p, a, nburns, walk_len, randPoints, - push_back_policy, rng); - randPoints.clear(); + if (nburns > 0) { + RandomPointGenerator::apply(P, p, a, nburns, walk_len, randPoints, + push_back_policy, rng); + randPoints.clear(); + } RandomPointGenerator::apply(P, p, a, rnum, walk_len, randPoints, push_back_policy, rng); @@ -205,13 +211,13 @@ void gaussian_sampling(PointList &randPoints, Point p = starting_point; typedef GaussianRandomPointGenerator RandomPointGenerator; - RandomPointGenerator::apply(P, p, a, nburns, walk_len, randPoints, - push_back_policy, rng, WalkType.param); - randPoints.clear(); + if (nburns > 0) { + RandomPointGenerator::apply(P, p, a, nburns, walk_len, randPoints, + push_back_policy, rng, WalkType.param); + randPoints.clear(); + } RandomPointGenerator::apply(P, p, a, rnum, walk_len, randPoints, push_back_policy, rng, WalkType.param); - - } template < @@ -264,15 +270,182 @@ void logconcave_sampling(PointList &randPoints, Point p = starting_point; typedef LogconcaveRandomPointGenerator RandomPointGenerator; - RandomPointGenerator::apply(P, p, nburns, walk_len, randPoints, - push_back_policy, rng, F, f, params); + if (nburns > 0) { + RandomPointGenerator::apply(P, p, nburns, walk_len, randPoints, + push_back_policy, rng, F, f, params); - randPoints.clear(); + randPoints.clear(); + } RandomPointGenerator::apply(P, p, rnum, walk_len, randPoints, push_back_policy, rng, F, f, params); } +template +< + typename WalkTypePolicy, + typename PointList, + typename Polytope, + typename RandomNumberGenerator, + typename NT, + typename Point +> +void exponential_sampling(PointList &randPoints, + Polytope &P, + RandomNumberGenerator &rng, + const unsigned int &walk_len, + const unsigned int &rnum, + const Point &c, + const NT &a, + const Point &starting_point, + unsigned int const& nburns) +{ + + typedef typename WalkTypePolicy::template Walk + < + Polytope, + RandomNumberGenerator + > walk; + + PushBackWalkPolicy push_back_policy; + + Point p = starting_point; + + typedef ExponentialRandomPointGenerator RandomPointGenerator; + if (nburns > 0) { + RandomPointGenerator::apply(P, p, c, a, nburns, walk_len, randPoints, + push_back_policy, rng); + randPoints.clear(); + } + RandomPointGenerator::apply(P, p, c, a, rnum, walk_len, randPoints, + push_back_policy, rng); +} + + +template < + typename PointList, + typename Polytope, + typename RandomNumberGenerator, + typename WalkTypePolicy, + typename NT, + typename Point + > +void exponential_sampling(PointList &randPoints, + Polytope &P, + RandomNumberGenerator &rng, + WalkTypePolicy &WalkType, + const unsigned int &walk_len, + const unsigned int &rnum, + const Point &c, + const NT &a, + const Point &starting_point, + unsigned int const& nburns) +{ + + typedef typename WalkTypePolicy::template Walk + < + Polytope, + RandomNumberGenerator + > walk; + + PushBackWalkPolicy push_back_policy; + + Point p = starting_point; + + typedef ExponentialRandomPointGenerator RandomPointGenerator; + if (nburns > 0) { + RandomPointGenerator::apply(P, p, c, a, nburns, walk_len, randPoints, + push_back_policy, rng, WalkType.param); + randPoints.clear(); + } + RandomPointGenerator::apply(P, p, c, a, rnum, walk_len, randPoints, + push_back_policy, rng, WalkType.param); +} + + +template +< + typename WalkTypePolicy, + typename PointList, + typename Polytope, + typename RandomNumberGenerator, + typename NT, + typename Point +> +void exponential_sampling(PointList &randPoints, + Polytope &P, + RandomNumberGenerator &rng, + const unsigned int &walk_len, + const unsigned int &rnum, + const Point &c, + const NT &a, + const NT &eta, + const Point &starting_point, + unsigned int const& nburns) +{ + + typedef typename WalkTypePolicy::template Walk + < + Polytope, + RandomNumberGenerator + > walk; + + PushBackWalkPolicy push_back_policy; + + Point p = starting_point; + + typedef ExponentialRandomPointGenerator RandomPointGenerator; + if (nburns > 0) { + RandomPointGenerator::apply(P, p, c, a, eta, nburns, walk_len, randPoints, + push_back_policy, rng); + randPoints.clear(); + } + RandomPointGenerator::apply(P, p, c, a, eta, rnum, walk_len, randPoints, + push_back_policy, rng); +} + + +template < + typename PointList, + typename Polytope, + typename RandomNumberGenerator, + typename WalkTypePolicy, + typename NT, + typename Point + > +void exponential_sampling(PointList &randPoints, + Polytope &P, + RandomNumberGenerator &rng, + WalkTypePolicy &WalkType, + const unsigned int &walk_len, + const unsigned int &rnum, + const Point &c, + const NT &a, + const NT &eta, + const Point &starting_point, + unsigned int const& nburns) +{ + + typedef typename WalkTypePolicy::template Walk + < + Polytope, + RandomNumberGenerator + > walk; + + PushBackWalkPolicy push_back_policy; + + Point p = starting_point; + + typedef ExponentialRandomPointGenerator RandomPointGenerator; + if (nburns > 0) { + RandomPointGenerator::apply(P, p, c, a, eta, nburns, walk_len, randPoints, + push_back_policy, rng, WalkType.param); + randPoints.clear(); + } + RandomPointGenerator::apply(P, p, c, a, eta, rnum, walk_len, randPoints, + push_back_policy, rng, WalkType.param); +} + #endif