diff --git a/R/HpolytopeSparseClass.R b/R/HpolytopeSparseClass.R new file mode 100644 index 00000000..394396f4 --- /dev/null +++ b/R/HpolytopeSparseClass.R @@ -0,0 +1,51 @@ +#' An R class to represent an H-polytope +#' +#' An H-polytope is a convex polytope defined by a set of linear inequalities or equivalently a \eqn{d}-dimensional H-polytope with \eqn{m} facets is defined by a \eqn{m\times d} matrix A and a \eqn{m}-dimensional vector b, s.t.: \eqn{Ax\leq b}. +#' +#' \describe{ +#' \item{A}{An \eqn{m\times d} numerical matrix.} +#' +#' \item{b}{An \eqn{m}-dimensional vector b.} +#' +#' \item{volume}{The volume of the polytope if it is known, \eqn{NaN} otherwise by default.} +#' +#' \item{type}{A character with default value 'Hpolytope', to declare the representation of the polytope.} +#' } +#' +#' @examples +#' A = matrix(c(-1,0,0,-1,1,1), ncol=2, nrow=3, byrow=TRUE) +#' b = c(0,0,1) +#' P = Hpolytope(A = A, b = b) +#' +#' @name HpolytopeSparse-class +#' @rdname HpolytopeSparse-class +#' @exportClass HpolytopeSparse +library(Matrix) +HpolytopeSparse <- setClass ( + # Class name + "HpolytopeSparse", + + # Defining slot type + representation ( + Aineq = "sparseMatrix", + bineq = "numeric", + Aeq = "sparseMatrix", + beq = "numeric", + lb = "numeric", + ub = "numeric", + dimension = "numeric", + type = "character" + ), + + # Initializing slots + prototype = list( + Aineq = as(0, "CsparseMatrix"), + bineq = as.numeric(NULL), + Aeq = as(0, "CsparseMatrix"), + beq = as.numeric(NULL), + lb = as.numeric(NULL), + ub = as.numeric(NULL), + dimension = as.numeric(NaN), + type = "HpolytopeSparse" + ) +) diff --git a/R/RcppExports.R b/R/RcppExports.R index f16e5f11..2ada6786 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -311,7 +311,17 @@ 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 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} }{ 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{'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). +#' 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()}.} @@ -360,7 +370,7 @@ rounding <- function(P, method = NULL, seed = NULL) { #' # gaussian distribution from the 2d unit simplex in H-representation with variance = 2 #' A = matrix(c(-1,0,0,-1,1,1), ncol=2, nrow=3, byrow=TRUE) #' b = c(0,0,1) -#' P = Hpolytope(A = A, b = b) +#' P = Hpolytope(A=A,b=b) #' points = sample_points(P, n = 100, distribution = list("density" = "gaussian", "variance" = 2)) #' #' # uniform points from the boundary of a 2-dimensional random H-polytope diff --git a/man/examples/nuts_rand_poly.R b/man/examples/nuts_rand_poly.R index 22187931..6c10914e 100644 --- a/man/examples/nuts_rand_poly.R +++ b/man/examples/nuts_rand_poly.R @@ -25,8 +25,8 @@ f <- function(x) (norm_vec(x)^2 + sum(x)) # Negative log-probability gradient oracle grad_f <- function(x) (2 * x + 1) -dimension <- 50 -facets <- 200 +dimension <- 5 +facets <- 20 # Create domain of truncation H <- gen_rand_hpoly(dimension, facets) @@ -39,11 +39,11 @@ P <- Hpolytope(A = Tr$Mat[1:nrow(Tr$Mat), 2:ncol(Tr$Mat)], b = Tr$Mat[,1]) x_min = matrix(0, dimension, 1) # Warm start point from truncated Gaussian -warm_start <- sample_points(P, n = 1, random_walk = list("nburns" = 5000), +warm_start <- sample_points(P, n = 1, random_walk = list("nburns" = 500), distribution = list("density" = "gaussian", "variance" = 1/2, "mode" = x_min)) # Sample points -n_samples <- 20000 +n_samples <- 1000 samples <- sample_points(P, n = n_samples, random_walk = list("walk" = "NUTS", "solver" = "leapfrog", "starting_point" = warm_start[,1]), distribution = list("density" = "logconcave", "negative_logprob" = f, "negative_logprob_gradient" = grad_f)) @@ -53,4 +53,6 @@ hist(samples[1,], probability=TRUE, breaks = 100) psrfs <- psrf_univariate(samples) n_ess <- ess(samples) - +cat("ESS=", n_ess, "\n") +cat("\n") +cat("PSRF=", psrfs, "\n") diff --git a/man/examples/simple_hmc_rand_poly.R b/man/examples/simple_hmc_rand_poly.R index 975062d4..bbba743b 100644 --- a/man/examples/simple_hmc_rand_poly.R +++ b/man/examples/simple_hmc_rand_poly.R @@ -25,8 +25,8 @@ f <- function(x) (norm_vec(x)^2 + sum(x)) # Negative log-probability gradient oracle grad_f <- function(x) (2 * x + 1) -dimension <- 50 -facets <- 200 +dimension <- 5 +facets <- 20 # Create domain of truncation H <- gen_rand_hpoly(dimension, facets) diff --git a/man/examples/sparse_crhmc.R b/man/examples/sparse_crhmc.R index 37337198..b3662435 100644 --- a/man/examples/sparse_crhmc.R +++ b/man/examples/sparse_crhmc.R @@ -31,7 +31,7 @@ A = matrix(c(1, -1), ncol=1, nrow=2, byrow=TRUE) b = c(2,1) # Create domain of truncation -P <- volesti::Hpolytope$new(A, b) +P <- volesti::Hpolytope(A=A, b=b) # Mode of logconcave density x_min <- c(-0.5) @@ -41,7 +41,7 @@ L <- 2 m <- 2 # Sample points -n_samples <- 80000 +n_samples <- 1000 n_burns <- n_samples / 2 cat("---Sampling without hessian\n") pts <- sample_points(P, n = n_samples, random_walk = list("walk" = "CRHMC", "step_size" = 0.3, "nburns" = n_burns, "walk_length" = 1, "solver" = "implicit_midpoint"), distribution = list("density" = "logconcave", "negative_logprob" = f, "negative_logprob_gradient" = grad_f, "L_" = L, "m" = m)) @@ -79,16 +79,16 @@ invisible(capture.output(dev.off())) walk="CRHMC" library(Matrix) -bineq=matrix(c(10,10,10,10,10), nrow=5, ncol=1, byrow=TRUE) +bineq=c(10,10,10,10,10) Aineq = matrix(c(1,0,-0.25,-1,2.5,1,0.4,-1,-0.9,0.5), nrow=5, ncol=2, byrow = TRUE) Aineq = as( Aineq, 'dgCMatrix' ) -beq=matrix(,nrow=0, ncol=1, byrow=TRUE) +beq=(0) Aeq = matrix(, nrow=0, ncol=2, byrow = TRUE) Aeq=as( Aeq, 'dgCMatrix' ) lb=-100000*c(1,1); ub=100000*c(1,1); cat("---Sampling the normal distribution in a pentagon\n") -P <- volesti::sparse_constraint_problem$new(Aineq, bineq,Aeq, beq, lb, ub) +P <- volesti::HpolytopeSparse(Aineq=Aineq, bineq=bineq, Aeq=Aeq, beq=beq, lb=lb, ub=ub) points <- sample_points(P, n = n_samples, random_walk = list("walk" = "CRHMC", "step_size" = 0.3, "nburns" = n_burns, "walk_length" = 1, "solver" = "implicit_midpoint"), distribution = list("density" = "logconcave", "variance" = 8)) jpeg("pentagon.jpg") plot(ggplot(data.frame( x=points[1,], y=points[2,] )) + diff --git a/src/sample_points.cpp b/src/sample_points.cpp index b349d1ba..b8e967bf 100644 --- a/src/sample_points.cpp +++ b/src/sample_points.cpp @@ -24,6 +24,7 @@ #include "ode_solvers/ode_solvers.hpp" #include "oracle_functors_rcpp.h" #include "preprocess/crhmc/constraint_problem.h" + enum random_walks { ball_walk, rdhr, @@ -57,7 +58,7 @@ void sample_from_polytope(Polytope &P, int type, RNGType &rng, PointList &randPo unsigned int const& walkL, unsigned int const& numpoints, bool const& gaussian, NT const& a, NT const& L, Point const& c, Point const& StartingPoint, unsigned int const& nburns, - bool const& set_L, NT const& eta, random_walks walk, + bool const& set_L, random_walks walk, NegativeGradientFunctor *F=NULL, NegativeLogprobFunctor *f=NULL, HessianFunctor *h=NULL, ode_solvers solver_type = no_solver) { @@ -217,7 +218,8 @@ void sample_from_polytope(Polytope &P, int type, RNGType &rng, PointList &randPo break; case crhmc:{ - execute_crhmc(P, rng, randPoints, walkL, numpoints, nburns, F, f, h); + execute_crhmc(P, rng, randPoints, walkL, numpoints, nburns, F, f, h); break; } case nuts: @@ -288,7 +290,7 @@ void sample_from_polytope(Polytope &P, int type, RNGType &rng, PointList &randPo //' \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{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{ @@ -330,7 +332,7 @@ void sample_from_polytope(Polytope &P, int type, RNGType &rng, PointList &randPo //' # gaussian distribution from the 2d unit simplex in H-representation with variance = 2 //' A = matrix(c(-1,0,0,-1,1,1), ncol=2, nrow=3, byrow=TRUE) //' b = c(0,0,1) -//' P = Hpolytope$new(A,b) +//' P = Hpolytope(A=A,b=b) //' points = sample_points(P, n = 100, distribution = list("density" = "gaussian", "variance" = 2)) //' //' # uniform points from the boundary of a 2-dimensional random H-polytope @@ -341,7 +343,7 @@ void sample_from_polytope(Polytope &P, int type, RNGType &rng, PointList &randPo //' //' @export // [[Rcpp::export]] -Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, +Rcpp::NumericMatrix sample_points(Rcpp::Reference P, Rcpp::Nullable n, Rcpp::Nullable random_walk = R_NilValue, Rcpp::Nullable distribution = R_NilValue, @@ -360,9 +362,6 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, typedef Eigen::SparseMatrix SpMat; typedef constraint_problem sparse_problem; - unsigned int type = Rcpp::as(P).field("type"), dim = Rcpp::as(P).field("dimension"), - walkL = 1, numpoints, nburns = 0; - RcppFunctor::GradientFunctor *F = NULL; RcppFunctor::FunctionFunctor *f = NULL; RcppFunctor::HessianFunctor *h = NULL; @@ -372,6 +371,34 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, GaussianFunctor::HessianFunctor *hess_g = NULL; bool functor_defined = true; + unsigned int dim; + unsigned int type; + + std::string type_str = Rcpp::as(P.slot("type")); + + if (type_str.compare(std::string("Hpolytope")) == 0) { + dim = Rcpp::as(P.slot("A")).cols(); + type = 1; + } else if (type_str.compare(std::string("Vpolytope")) == 0) { + dim = Rcpp::as(P.slot("V")).cols(); + type = 2; + } else if (type_str.compare(std::string("Zonotope")) == 0) { + dim = Rcpp::as(P.slot("G")).cols(); + type = 3; + } else if (type_str.compare(std::string("VpolytopeIntersection")) == 0) { + dim = Rcpp::as(P.slot("V1")).cols(); + type = 4; + } else if (type_str.compare(std::string("HpolytopeSparse")) == 0) { + dim = Rcpp::as(P.slot("Aineq")).cols(); + type = 5; + } else { + throw Rcpp::exception("Unknown polytope representation!"); + } + + unsigned int numpoints; + unsigned int nburns = 0; + unsigned int walkL = 1; + RNGType rng(dim); if (seed.isNotNull()) { @@ -668,8 +695,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, switch(type) { case 1: { // Hpolytope - Hpolytope HP(dim, Rcpp::as(Rcpp::as(P).field("A")), - Rcpp::as(Rcpp::as(P).field("b"))); + 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."); @@ -686,18 +712,17 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, } if (functor_defined) { sample_from_polytope(HP, type, rng, randPoints, walkL, numpoints, gaussian, a, L, c, - StartingPoint, nburns, set_L, eta, walk, F, f, h, solver); + StartingPoint, nburns, set_L, walk, F, f, h, solver); } else { sample_from_polytope(HP, type, rng, randPoints, walkL, numpoints, gaussian, a, L, c, - StartingPoint, nburns, set_L, eta, walk, G, g, hess_g, solver); + StartingPoint, nburns, set_L, walk, G, g, hess_g, solver); } break; } case 2: { // Vpolytope - Vpolytope VP(dim, Rcpp::as(Rcpp::as(P).field("V")), - VT::Ones(Rcpp::as(Rcpp::as(P).field("V")).rows())); + Vpolytope VP(dim, Rcpp::as(P.slot("V")), VT::Ones(Rcpp::as(P.slot("V")).rows())); InnerBall = VP.ComputeInnerBall(); if (InnerBall.second < 0.0) throw Rcpp::exception("Unable to compute a feasible point."); @@ -712,13 +737,12 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, VP.shift(mode.getCoefficients()); } sample_from_polytope(VP, type, rng, randPoints, walkL, numpoints, gaussian, a, L, c, - StartingPoint, nburns, set_L, eta, walk, F, f, h, solver); + StartingPoint, nburns, set_L, walk, F, f, h, solver); break; } case 3: { // Zonotope - zonotope ZP(dim, Rcpp::as(Rcpp::as(P).field("G")), - VT::Ones(Rcpp::as(Rcpp::as(P).field("G")).rows())); + zonotope ZP(dim, Rcpp::as(P.slot("G")), VT::Ones(Rcpp::as(P.slot("G")).rows())); InnerBall = ZP.ComputeInnerBall(); if (InnerBall.second < 0.0) throw Rcpp::exception("Unable to compute a feasible point."); @@ -733,15 +757,15 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, ZP.shift(mode.getCoefficients()); } sample_from_polytope(ZP, type, rng, randPoints, walkL, numpoints, gaussian, a, L, c, - StartingPoint, nburns, set_L, eta, walk, F, f, h, solver); + StartingPoint, nburns, set_L, walk, F, f, h, solver); break; } case 4: { // Intersection of two V-polytopes - Vpolytope VP1(dim, Rcpp::as(Rcpp::as(P).field("V1")), - VT::Ones(Rcpp::as(Rcpp::as(P).field("V1")).rows())); - Vpolytope VP2(dim, Rcpp::as(Rcpp::as(P).field("V2")), - VT::Ones(Rcpp::as(Rcpp::as(P).field("V2")).rows())); + Vpolytope VP1(dim, Rcpp::as(P.slot("V1")), + VT::Ones(Rcpp::as(P.slot("V1")).rows())); + Vpolytope VP2(dim, Rcpp::as(P.slot("V2")), + VT::Ones(Rcpp::as(P.slot("V2")).rows())); InterVP VPcVP(VP1, VP2); if (!VPcVP.is_feasible()) throw Rcpp::exception("Empty set!"); @@ -756,26 +780,30 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, VPcVP.shift(mode.getCoefficients()); } sample_from_polytope(VPcVP, type, rng, randPoints, walkL, numpoints, gaussian, a, L, c, - StartingPoint, nburns, set_L, eta, walk, F, f, h, solver); + StartingPoint, nburns, set_L, walk, F, f, h, solver); break; } case 5: { - // Sparse constraint_problem - SpMat Aeq = Rcpp::as(Rcpp::as(P).field("Aeq")); - VT beq= Rcpp::as(Rcpp::as(P).field("beq")); - SpMat Aineq = Rcpp::as(Rcpp::as(P).field("Aineq")); - VT bineq= Rcpp::as(Rcpp::as(P).field("bineq")); - VT lb= Rcpp::as(Rcpp::as(P).field("lb")); - VT ub= Rcpp::as(Rcpp::as(P).field("ub")); - sparse_problem problem(dim, Aeq, beq, Aineq, bineq, lb, ub); - if(walk!=crhmc){throw Rcpp::exception("Sparse problems are supported only by the CRHMC walk.");} - if (functor_defined) { - execute_crhmc, RcppFunctor::GradientFunctor,RcppFunctor::FunctionFunctor, RcppFunctor::HessianFunctor, CRHMCWalk, 8>(problem, rng, randPoints, walkL, numpoints, nburns, F, f, h); - } - else { - execute_crhmc, GaussianFunctor::GradientFunctor,GaussianFunctor::FunctionFunctor, GaussianFunctor::HessianFunctor, CRHMCWalk, 8>(problem, rng, randPoints, walkL, numpoints, nburns, G, g, hess_g); - } - break; + // Sparse constraint_problem + SpMat Aeq = Rcpp::as(P.slot("Aeq")); + VT beq= Rcpp::as(P.slot("beq")); + SpMat Aineq = Rcpp::as(P.slot("Aineq")); + VT bineq= Rcpp::as(P.slot("bineq")); + VT lb= Rcpp::as(P.slot("lb")); + VT ub= Rcpp::as(P.slot("ub")); + sparse_problem problem(dim, Aeq, beq, Aineq, bineq, lb, ub); + if(walk!=crhmc){throw Rcpp::exception("Sparse problems are supported only by the CRHMC walk.");} + if (functor_defined) { + execute_crhmc, RcppFunctor::GradientFunctor, + RcppFunctor::FunctionFunctor, RcppFunctor::HessianFunctor, CRHMCWalk, 1> + (problem, rng, randPoints, walkL, numpoints, nburns, F, f, h); + } + else { + execute_crhmc, GaussianFunctor::GradientFunctor, + GaussianFunctor::FunctionFunctor, GaussianFunctor::HessianFunctor, CRHMCWalk, 1> + (problem, rng, randPoints, walkL, numpoints, nburns, G, g, hess_g); + } + break; } }