diff --git a/man/examples/boundary_samplers.R b/man/examples/boundary_samplers.R new file mode 100644 index 00000000..47dba6e5 --- /dev/null +++ b/man/examples/boundary_samplers.R @@ -0,0 +1,103 @@ +# VolEsti (volume computation and sampling library) + +# Copyright (c) 2012-2025 Vissarion Fisikopoulos +# Copyright (c) 2018-2025 Apostolos Chalkis +# Copyright (c) 2025-2025 Iva Janković + +# Contributed and/or modified by Iva Janković, as part of Google Summer of Code 2025 program. + +# Licensed under GNU LGPL.3, see LICENCE file + +# Example script for comparing volesti boundary samplers + +# Import required libraries +library(ggplot2) +library(volesti) +library(R.matlab) + +# Generate EColi polytope + +root <- rprojroot::find_root_file(criterion = rprojroot::has_file("DESCRIPTION")) +metabolic_polytope_mat <- readMat(paste(root , '/man/examples/data/polytope_e_coli.mat', sep="")) +A <- as.matrix(metabolic_polytope_mat$polytope[[1]]) +b <- as.matrix(metabolic_polytope_mat$polytope[[2]]) +center <- as.matrix(metabolic_polytope_mat$polytope[[3]]) +b_new <- b - A %*% center +P <- Hpolytope(A = A, b = c(b_new)) + +# Sampling time (Boundary Random Direction Hit and Run) +start_time_brdhr <- Sys.time() +points_brdhr= sample_points(P, n = 100000, random_walk = list("walk" = "BRDHR", "nburns" = 10, "walk_length" = 1)) +end_time_brdhr <- Sys.time() + +# Calculate Effective Sample size (Boundary Random Direction Hit and Run) +brdhr_ess = ess(points_brdhr) +min_ess_brdhr <- min(brdhr_ess) + +# Calculate PSRF (Boundary Random Direction Hit and Run) +brdhr_psrfs = psrf_univariate(points_brdhr) +max_psrf_brdhr = max(brdhr_psrfs) +elapsed_time_brdhr <- end_time_brdhr - start_time_brdhr +time_per_ind_brdhr <- elapsed_time_brdhr / min_ess_brdhr + +# +# Sampling time (Boundary Coordinate Direction Hit and Run) +start_time_bcdhr <- Sys.time() +points_bcdhr= sample_points(P, n = 100000, random_walk = list("walk" = "BCDHR", "nburns" = 10, "walk_length" = 1)) +end_time_bcdhr <- Sys.time() + +# Calculate Effective Sample size (Boundary Coordinate Direction Hit and Run) +bcdhr_ess = ess(points_bcdhr) +min_ess_bcdhr <- min(bcdhr_ess) + +# Calculate PSRF (Boundary Coordinate Direction Hit and Run) +bcdhr_psrfs = psrf_univariate(points_bcdhr) +max_psrf_bcdhr = max(bcdhr_psrfs) +elapsed_time_bcdhr <- end_time_bcdhr - start_time_bcdhr +time_per_ind_bcdhr <- elapsed_time_bcdhr / min_ess_bcdhr + +# + +# Sampling time (Shake and Bake) +start_time_sb <- Sys.time() +points_sb = sample_points(P, n = 100000, random_walk = list("walk" = "SB", "nburns" = 10, "walk_length" = 1)) +end_time_sb <- Sys.time() + +# Calculate Effective Sample size (Shake and Bake) +sb_ess = ess(points_sb) +min_ess_sb <- min(sb_ess) + +# Calculate PSRF (Shake and Bake) +sb_psrfs = psrf_univariate(points_sb) +max_psrf_sb = max(sb_psrfs) +elapsed_time_sb <- end_time_sb - start_time_sb +time_per_ind_sb <- elapsed_time_sb / min_ess_sb + +# + +# Sampling time (Billiard Shake and Bake) +start_time_bsb <- Sys.time() +points_bsb = sample_points(P, n = 100000, random_walk = list("walk" = "BSB", "nburns" = 10, "walk_length" = 1)) +end_time_bsb <- Sys.time() + +# Calculate Effective Sample size (Billiard Shake and Bake) +bsb_ess = ess(points_bsb) +min_ess_bsb <- min(bsb_ess) + +# Calculate PSRF (Billiard Shake and Bake) +bsb_psrfs = psrf_univariate(points_bsb) +max_psrf_bsb = max(bsb_psrfs) +elapsed_time_bsb <- end_time_bsb - start_time_bsb +time_per_ind_bsb <- elapsed_time_bsb / min_ess_bsb + + +#Final results + +results <- data.frame( + Method = c("Boundary RDHR", "Boundary CDHR","Shake and Bake", "Billiard Shake and Bake"), + Min_ESS = c(min_ess_brdhr, min_ess_bcdhr,min_ess_sb, min_ess_bsb), + Max_PSRF = c(max_psrf_brdhr, max_psrf_bcdhr, max_psrf_sb, max_psrf_bsb), + Time_per_ind_sample = c(time_per_ind_brdhr, time_per_ind_bcdhr,time_per_ind_sb, time_per_ind_bsb) +) + +print(results) diff --git a/src/sample_points.cpp b/src/sample_points.cpp index ab61f49a..1949340f 100644 --- a/src/sample_points.cpp +++ b/src/sample_points.cpp @@ -9,6 +9,7 @@ // Contributed and/or modified by Alexandros Manochis, as part of Google Summer of Code 2020 program. +// Contributed and/or modified by Iva Janković, as part of Google Summer of Code 2025 program. #include #include @@ -24,6 +25,7 @@ #include "ode_solvers/ode_solvers.hpp" #include "oracle_functors_rcpp.h" #include "preprocess/crhmc/constraint_problem.h" +#include "preprocess/feasible_point.hpp" enum random_walks { ball_walk, @@ -31,6 +33,8 @@ enum random_walks { cdhr, billiard, accelarated_billiard, + shake_and_bake, + billiard_shake_bake, dikin_walk, vaidya_walk, john_walk, @@ -60,7 +64,8 @@ void sample_from_polytope(Polytope &P, int type, RNGType &rng, PointList &randPo Point const& StartingPoint, unsigned int const& nburns, bool const& set_L, random_walks walk, NegativeGradientFunctor *F=NULL, NegativeLogprobFunctor *f=NULL, - HessianFunctor *h=NULL, ode_solvers solver_type = no_solver) + HessianFunctor *h=NULL, ode_solvers solver_type = no_solver, + int facet_index=-1, int nreflections=-1) { switch (walk) { @@ -135,6 +140,14 @@ void sample_from_polytope(Polytope &P, int type, RNGType &rng, PointList &randPo StartingPoint, nburns); } break; + case shake_and_bake: + shakeandbake_sampling(randPoints, P, rng, walkL, numpoints, + StartingPoint, nburns, facet_index); + break; + case billiard_shake_bake: + billiard_shakeandbake_sampling(randPoints, P, rng, walkL,nreflections, numpoints, + StartingPoint, nburns, facet_index); + break; case gaussian_hmc: if(set_L) { GaussianHamiltonianMonteCarloExactWalk WalkType(L); @@ -232,18 +245,19 @@ bool is_walk(Rcpp::Nullable random_walk, std::string str) { //' x) \code{'NUTS'} for NUTS Hamiltonian Monte Carlo sampler (logconcave densities), xi) \code{'HMC'} for Hamiltonian Monte Carlo (logconcave densities), //' xii) CRHMC for Riemannian HMC with H-polytope constraints (uniform and general logconcave densities), //' xiii) \code{'ULD'} for Underdamped Langevin Dynamics using the Randomized Midpoint Method (logconcave densities), -//' xiii) \code{'ExactHMC'} for exact Hamiltonian Monte Carlo with reflections (spherical Gaussian or exponential distribution). +//' xiii) \code{'ExactHMC'} for exact Hamiltonian Monte Carlo with reflections (spherical Gaussian or exponential distribution) +//' xiv) \code{'SB'} for Running variant of Shake and Bake algorithm and \code{'BSB'} for hybrid algorithm = Billiard + Running Shake and Bake, both for uniform boundary sampling. //' The default walk is \code{'aBiW'} for the uniform distribution, \code{'CDHR'} for the Gaussian distribution and H-polytopes and //' \code{'BiW'} or \code{'RDHR'} for the same distributions and V-polytopes and zonotopes. \code{'NUTS'} is the default sampler for logconcave densities and \code{'CRHMC'} //' for logconcave densities with H-polytope and sparse constrainted problems.} //' \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{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()}.Exception is for (Billiard) Shake and Bake where starting point is on the boundary of the polytope!} //' \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{nreflections}}{Optionally chosen an upper bound for the number of reflections for the Billiard SB} //' @param distribution Optional. A list that declares the target density and some related parameters as follows: //' \describe{ //' \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.} @@ -275,6 +289,9 @@ bool is_walk(Rcpp::Nullable random_walk, std::string str) { //' @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.} //' +//' @references \cite{Boender CGE, Caron Richard J, McDonald J Fred, Kan AHG Rinnooy, Romeijn H Edwin, Smith Robert L, Telgen Jan, Vorst ACF +//' \dQuote{"hake-and-bake algorithms for generating uniform points on the boundary of bounded polyhedra,"} \emph{Operations Research,} 1991.}. +//' //' @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 @@ -304,7 +321,8 @@ Rcpp::NumericMatrix sample_points(Rcpp::Reference P, typedef double NT; typedef Cartesian Kernel; typedef BoostRandomNumberGenerator RNGType; - typedef typename Kernel::Point Point; + typedef typename Kernel::Point Point; + typedef typename Point::FT FT; typedef HPolytope Hpolytope; typedef VPolytope Vpolytope; typedef Zonotope zonotope; @@ -359,6 +377,8 @@ Rcpp::NumericMatrix sample_points(Rcpp::Reference P, NT radius = 1.0; NT L; NT eta = -1; + int facet_index; // Used only for Shake and Bake sampling + int nreflections; // Used only for Billiard Shake and Bake sampling bool set_mode = false; bool gaussian = false; @@ -523,7 +543,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Reference P, } else if (logconcave) { walk = (type == 5) ? crhmc : nuts; } else if (gaussian) { - walk = (type == 1) ? cdhr : rdhr; + walk = (type == 1) ? cdhr : rdhr; } else { walk = (type == 1) ? accelarated_billiard : billiard; } @@ -575,6 +595,22 @@ Rcpp::NumericMatrix sample_points(Rcpp::Reference P, set_L = true; if (L <= 0.0) throw Rcpp::exception("L must be a postitive number!"); } + } else if (is_walk(random_walk, std::string("SB"))) { + if (gaussian || exponential) throw Rcpp::exception("Gaussian/ exponential Shake and Bake sampling is not supported!"); + if (type !=1) { + throw Rcpp::exception("Shake and Bake sampling is supported only for H-polytopes !"); + } + walk = shake_and_bake; + } else if (is_walk(random_walk, std::string("BSB"))) { + if (gaussian || exponential) throw Rcpp::exception("Gaussian/ exponential Billiard Shake and Bake sampling is not supported!"); + if (type !=1) { + throw Rcpp::exception("Billiard Shake and Bake sampling is supported only for H-polytopes !"); + } + if (Rcpp::as(random_walk).containsElementNamed("nreflections")) { + nreflections = Rcpp::as(Rcpp::as(random_walk)["nreflections"]); + if (nreflections <= 0.0) throw Rcpp::exception("Number of reflections must be a postitive number!"); + } + walk = billiard_shake_bake; } else if (is_walk(random_walk, std::string("BRDHR"))) { if (gaussian || exponential) throw Rcpp::exception("Gaussian sampling from the boundary is not supported!"); walk = brdhr; @@ -639,26 +675,37 @@ Rcpp::NumericMatrix sample_points(Rcpp::Reference P, // Hpolytope Hpolytope HP(dim, Rcpp::as(P.slot("A")), Rcpp::as(P.slot("b"))); - InnerBall = HP.ComputeInnerBall(); - if (InnerBall.second < 0.0) throw Rcpp::exception("Unable to compute a feasible point."); - if (!set_starting_point || (!set_mode && gaussian)) { - if (!set_starting_point) StartingPoint = InnerBall.first; - if (!set_mode && gaussian) mode = InnerBall.first; + if ((walk == shake_and_bake) || (walk == billiard_shake_bake)) { + if (!set_starting_point){ + FT tolerance = static_cast(1e-10); + auto results = compute_boundary_point(HP, rng, tolerance); + VT boundary_point= results.first; + facet_index=results.second; + StartingPoint=boundary_point; + } } - if (HP.is_in(StartingPoint) == 0) { - throw Rcpp::exception("The given point is not in the interior of the polytope!"); - } - if (gaussian) { - StartingPoint = StartingPoint - mode; - HP.shift(mode.getCoefficients()); + else{ + InnerBall = HP.ComputeInnerBall(); + if (InnerBall.second < 0.0) throw Rcpp::exception("Unable to compute a feasible point."); + if (!set_starting_point || (!set_mode && gaussian)) { + if (!set_starting_point) StartingPoint = InnerBall.first; + if (!set_mode && gaussian) mode = InnerBall.first; + } + if (HP.is_in(StartingPoint) == 0) { + throw Rcpp::exception("The given point is not in the interior of the polytope!"); + } + if (gaussian) { + StartingPoint = StartingPoint - mode; + HP.shift(mode.getCoefficients()); + } } if (functor_defined) { sample_from_polytope(HP, type, rng, randPoints, walkL, numpoints, gaussian, a, L, c, - StartingPoint, nburns, set_L, walk, F, f, h, solver); + StartingPoint, nburns, set_L, walk, F, f, h, solver,facet_index,nreflections); } else { sample_from_polytope(HP, type, rng, randPoints, walkL, numpoints, gaussian, a, L, c, - StartingPoint, nburns, set_L, walk, G, g, hess_g, solver); + StartingPoint, nburns, set_L, walk, G, g, hess_g, solver,facet_index,nreflections); } break; }