diff --git a/R-proj/R/RcppExports.R b/R-proj/R/RcppExports.R index d15f28082..898e40401 100644 --- a/R-proj/R/RcppExports.R +++ b/R-proj/R/RcppExports.R @@ -46,6 +46,7 @@ copula <- function(r1 = NULL, r2 = NULL, sigma = NULL, m = NULL, n = NULL) { #' \item{\code{radius} }{ The radius of the \eqn{d}-dimensional hypersphere. The default value is \eqn{1}.} #' } #' @param n The number of points that the function is going to sample. +#' @param seed Optional. A fixed seed for the number generator. #' #' @references \cite{R.Y. Rubinstein and B. Melamed, #' \dQuote{Modern simulation and modeling} \emph{ Wiley Series in Probability and Statistics,} 1998.} @@ -57,8 +58,8 @@ copula <- function(r1 = NULL, r2 = NULL, sigma = NULL, m = NULL, n = NULL) { #' # 100 uniform points from the 2-d unit ball #' points = direct_sampling(n = 100, body = list("type" = "ball", "dimension" = 2)) #' @export -direct_sampling <- function(body = NULL, n = NULL) { - .Call(`_volesti_direct_sampling`, body, n) +direct_sampling <- function(body = NULL, n = NULL, seed = NULL) { + .Call(`_volesti_direct_sampling`, body, n, seed) } #' Compute the exact volume of (a) a zonotope (b) an arbitrary simplex in V-representation or (c) if the volume is known and declared by the input object. @@ -171,13 +172,14 @@ rotating <- function(P, T = NULL, seed = NULL) { #' Internal rcpp function for the rounding of a convex polytope #' #' @param P A convex polytope (H- or V-representation or zonotope). +#' @param seed Optional. A fixed seed for the number generator. #' #' @section warning: #' Do not use this function. #' #' @return A numerical matrix that describes the rounded polytope and contains the round value. -rounding <- function(P) { - .Call(`_volesti_rounding`, P) +rounding <- function(P, seed = NULL) { + .Call(`_volesti_rounding`, P, seed) } #' Sample uniformly or normally distributed points from a convex Polytope (H-polytope, V-polytope or a zonotope). @@ -201,6 +203,7 @@ rounding <- function(P) { #' \item{\code{variance} }{ The variance of the multidimensional spherical gaussian. 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 Chebychev ball.} #' } +#' @param seed Optional. A fixed seed for the number generator. #' #' @return A \eqn{d\times n} matrix that contains, column-wise, the sampled points from the convex polytope P. #' @examples @@ -219,8 +222,8 @@ rounding <- function(P) { #' points = sample_points(P, n = 5000, random_walk = list("walk" = "BRDHR")) #' #' @export -sample_points <- function(P = NULL, n = NULL, random_walk = NULL, distribution = NULL) { - .Call(`_volesti_sample_points`, P, n, random_walk, distribution) +sample_points <- function(P = NULL, n = NULL, random_walk = NULL, distribution = NULL, seed = NULL) { + .Call(`_volesti_sample_points`, P, n, random_walk, distribution, seed) } #' The main function for volume approximation of a convex Polytope (H-polytope, V-polytope or a zonotope) @@ -234,24 +237,11 @@ sample_points <- function(P = NULL, n = NULL, random_walk = NULL, distribution = #' \item{\code{error} }{ A numeric value to set the upper bound for the approximation error. The default value is \eqn{1} for \code{'SOB'} and \eqn{0.1} otherwise.} #' \item{\code{random_walk} }{ A string that declares the random walk method: a) \code{'CDHR'} for Coordinate Directions Hit-and-Run, b) \code{'RDHR'} for Random Directions Hit-and-Run, c) \code{'BaW'} for Ball Walk, or \code{'BiW'} for Billiard walk. The default walk is \code{'CDHR'} for H-polytopes and \code{'BiW'} for the other representations.} #' \item{\code{walk_length} }{ An integer to set the number of the steps for the random walk. The default value is \eqn{\lfloor 10 + d/10\rfloor} for \code{'SOB'} and \eqn{1} otherwise.} -#' \item{\code{inner_ball} }{ A \eqn{d+1} numeric vector that contains an inner ball. The first \eqn{d} coordinates corresponds to the center and the last one to the radius of the ball. If it is not given then for H-polytopes the Chebychev ball is computed, for V-polytopes \eqn{d+1} vertices are picked randomly and the Chebychev ball of the defined simplex is computed. For a zonotope that is defined by the Minkowski sum of \eqn{m} segments we compute the maximal \eqn{r} s.t.: \eqn{re_i\in Z} for all \eqn{i=1,\dots ,d}, then the ball centered at the origin with radius \eqn{r/\sqrt{d}} is an inscribed ball.} -#' \item{\code{len_win} }{ The length of the sliding window for CG algorithm. The default value is \eqn{500+4dimension^2}.} -#' \item{\code{C} }{ A constant for the lower bound of \eqn{variance/mean^2} in schedule annealing of CG algorithm. The default value is \eqn{2}.} -#' \item{\code{M} }{ The number of points we sample in each step of schedule annealing in CG algorithm. The default value is \eqn{500C + dimension^2 / 2}.} -#' \item{\code{ratio} }{ Parameter of schedule annealing of CG algorithm, larger ratio means larger steps in schedule annealing. The default value is \eqn{1 - 1/dimension}.} -#' \item{\code{frac} }{ The fraction of the total error to spend in the first gaussian in CG algorithm. The default value is \eqn{0.1}.} -#' \item{\code{BW_rad} }{ The radius for the ball walk. The default value is \eqn{4r/dimension}, where \eqn{r} is the radius of the inscribed ball of the polytope.} -#' \item{\code{ub} }{ The lower bound for the ratios in MMC in CB algorithm. The default value is \eqn{0.1}.} -#' \item{\code{lb} }{ The upper bound for the ratios in MMC in CB algorithm. The default value is \eqn{0.15}.} -#' \item{\code{N} }{ An integer that controls the number of points \eqn{\nu N} generated in each convex body in annealing schedule of algorithm CB.} -#' \item{\code{nu} }{ The degrees of freedom for the t-student distribution in t-tests in CB algorithm. The default value is \eqn{10}.} -#' \item{\code{alpha} }{ The significance level for the t-tests in CB algorithm. The default values is 0.2.} -#' \item{\code{prob} }{ The probability is used for the empirical confidence interval in ratio estimation of CB algorithm. The default value is \eqn{0.75}.} +#' \item{\code{win_len} }{ The length of the sliding window for CG algorithm. The default value is \eqn{500+4dimension^2}.} #' \item{\code{hpoly} }{ A boolean parameter to use H-polytopes in MMC of CB algorithm. The default value is \code{FALSE}.} -#' \item{\code{minmaxW} }{ A boolean parameter to use the sliding window with a minmax values as a stopping criterion.} -#' \item{\code{L} }{The maximum length of the billiard trajectory.} #' } #' @param rounding Optional. A boolean parameter for rounding. The default value is \code{TRUE} for V-polytopes and \code{FALSE} otherwise. +#' @param seed Optional. A fixed seed for the number generator. #' #' @references \cite{I.Z.Emiris and V. Fisikopoulos, #' \dQuote{Practical polytope volume approximation,} \emph{ACM Trans. Math. Soft.,} 2014.}, @@ -274,8 +264,8 @@ sample_points <- function(P = NULL, n = NULL, random_walk = NULL, distribution = #' Z = gen_rand_zonotope(2, 4) #' vol = volume(Z, settings = list("random_walk" = "RDHR", "walk_length" = 5)) #' @export -volume <- function(P, settings = NULL, rounding = NULL) { - .Call(`_volesti_volume`, P, settings, rounding) +volume <- function(P, settings = NULL, rounding = NULL, seed = NULL) { + .Call(`_volesti_volume`, P, settings, rounding, seed) } #' An internal Rccp function for the over-approximation of a zonotope @@ -283,12 +273,13 @@ volume <- function(P, settings = NULL, rounding = NULL) { #' @param Z A zonotope. #' @param fit_ratio Optional. A boolean parameter to request the computation of the ratio of fitness. #' @param settings Optional. A list that declares the values of the parameters of CB algorithm. +#' @param seed Optional. A fixed seed for the number generator. #' #' @section warning: #' Do not use this function. #' #' @return A List that contains a numerical matrix that describes the PCA approximation as a H-polytope and the ratio of fitness. -zono_approx <- function(Z, fit_ratio = NULL, settings = NULL) { - .Call(`_volesti_zono_approx`, Z, fit_ratio, settings) +zono_approx <- function(Z, fit_ratio = NULL, settings = NULL, seed = NULL) { + .Call(`_volesti_zono_approx`, Z, fit_ratio, settings, seed) } diff --git a/R-proj/R/round_polytope.R b/R-proj/R/round_polytope.R index 562f91a83..53a5f06e4 100644 --- a/R-proj/R/round_polytope.R +++ b/R-proj/R/round_polytope.R @@ -3,6 +3,7 @@ #' Given a convex H or V polytope or a zonotope as input this functionbrings the polytope in well rounded position based on minimum volume enclosing ellipsoid of a pointset. #' #' @param P A convex polytope. It is an object from class (a) Hpolytope or (b) Vpolytope or (c) Zonotope. +#' @param seed Optional. A fixed seed for the number generator. #' #' @return A list with 2 elements: (a) a polytope of the same class as the input polytope class and (b) the element "round_value" which is the determinant of the square matrix of the linear transformation that was applied on the polytope that is given as input. #' diff --git a/R-proj/R/zonotope_approximation.R b/R-proj/R/zonotope_approximation.R index 33931562e..a66f05391 100644 --- a/R-proj/R/zonotope_approximation.R +++ b/R-proj/R/zonotope_approximation.R @@ -7,22 +7,11 @@ #' @param settings Optional. A list that declares the values of the parameters of CB algorithm as follows: #' \itemize{ #' \item{\code{error} }{ A numeric value to set the upper bound for the approximation error. The default value is \eqn{1} for \code{'SOB'} and \eqn{0.1} otherwise.} -#' \item{\code{random_walk} }{ A string that declares the random walk method: a) \code{'CDHR'} for Coordinate Directions Hit-and-Run, b) \code{'RDHR'} for Random Directions Hit-and-Run, c) \code{'BaW'} for Ball Walk, or \code{'BiW'} for Billiard walk. The default walk is \code{'CDHR'} for H-polytopes and \code{'BiW'} for the other representations.} #' \item{\code{walk_length} }{ An integer to set the number of the steps for the random walk. The default value is \eqn{\lfloor 10 + d/10\rfloor} for \code{'SOB'} and \eqn{1} otherwise.} -#' \item{\code{rounding} }{ A boolean parameter for rounding. The default value is \code{FALSE}.} -#' \item{\code{inner_ball} }{ A \eqn{d+1} numeric vector that contains an inner ball. The first \eqn{d} coordinates corresponds to the center and the last one to the radius of the ball. If it is not given then for H-polytopes the Chebychev ball is computed, for V-polytopes \eqn{d+1} vertices are picked randomly and the Chebychev ball of the defined simplex is computed. For a zonotope that is defined by the Minkowski sum of \eqn{m} segments we compute the maximal \eqn{r} s.t.: \eqn{re_i\in Z} for all \eqn{i=1,\dots ,d}, then the ball centered at the origin with radius \eqn{r/\sqrt{d}} is an inscribed ball.} -#' \item{\code{len_win} }{ The length of the sliding window for CG algorithm. The default value is \eqn{500+4dimension^2}.} -#' \item{\code{BW_rad} }{ The radius for the ball walk. The default value is \eqn{4r/dimension}, where \eqn{r} is the radius of the inscribed ball of the polytope.} -#' \item{\code{ub} }{ The lower bound for the ratios in MMC in CB algorithm. The default value is \eqn{0.1}.} -#' \item{\code{lb} }{ The upper bound for the ratios in MMC in CB algorithm. The default value is \eqn{0.15}.} -#' \item{\code{N} }{ An integer that controls the number of points \eqn{\nu N} generated in each convex body in annealing schedule of algorithm CB.} -#' \item{\code{nu} }{ The degrees of freedom for the t-student distribution in t-tests in CB algorithm. The default value is \eqn{10}.} -#' \item{\code{alpha} }{ The significance level for the t-tests in CB algorithm. The default values is 0.2.} -#' \item{\code{prob} }{ The probability is used for the empirical confidence interval in ratio estimation of CB algorithm. The default value is \eqn{0.75}.} +#' \item{\code{win_len} }{ The length of the sliding window for CG algorithm. The default value is \eqn{500+4dimension^2}.} #' \item{\code{hpoly} }{ A boolean parameter to use H-polytopes in MMC of CB algorithm. The default value is \code{FALSE}.} -#' \item{\code{minmaxW} }{ A boolean parameter to use the sliding window with a minmax values as a stopping criterion.} -#' \item{\code{L} }{The maximum length of the billiard trajectory.} #' } +#' @param seed Optional. A fixed seed for the number generator. #' #' @return A list that contains the approximation body in H-representation and the ratio of fitness #' @@ -32,9 +21,9 @@ #' retList = zonotope_approximation(Z = Z, fit_ratio = TRUE) #' #' @export -zonotope_approximation <- function(Z, fit_ratio = NULL, settings = NULL){ +zonotope_approximation <- function(Z, fit_ratio = NULL, settings = NULL, seed = NULL){ - ret_list = zono_approx(Z, fit_ratio, settings) + ret_list = zono_approx(Z, fit_ratio, settings, seed) Mat = ret_list$Mat diff --git a/R-proj/man/volume.Rd b/R-proj/man/volume.Rd index 9104a62a9..2b51d4d11 100644 --- a/R-proj/man/volume.Rd +++ b/R-proj/man/volume.Rd @@ -15,22 +15,8 @@ volume(P, settings = NULL, rounding = NULL) \item{\code{error} }{ A numeric value to set the upper bound for the approximation error. The default value is \eqn{1} for \code{'SOB'} and \eqn{0.1} otherwise.} \item{\code{random_walk} }{ A string that declares the random walk method: a) \code{'CDHR'} for Coordinate Directions Hit-and-Run, b) \code{'RDHR'} for Random Directions Hit-and-Run, c) \code{'BaW'} for Ball Walk, or \code{'BiW'} for Billiard walk. The default walk is \code{'CDHR'} for H-polytopes and \code{'BiW'} for the other representations.} \item{\code{walk_length} }{ An integer to set the number of the steps for the random walk. The default value is \eqn{\lfloor 10 + d/10\rfloor} for \code{'SOB'} and \eqn{1} otherwise.} -\item{\code{inner_ball} }{ A \eqn{d+1} numeric vector that contains an inner ball. The first \eqn{d} coordinates corresponds to the center and the last one to the radius of the ball. If it is not given then for H-polytopes the Chebychev ball is computed, for V-polytopes \eqn{d+1} vertices are picked randomly and the Chebychev ball of the defined simplex is computed. For a zonotope that is defined by the Minkowski sum of \eqn{m} segments we compute the maximal \eqn{r} s.t.: \eqn{re_i\in Z} for all \eqn{i=1,\dots ,d}, then the ball centered at the origin with radius \eqn{r/\sqrt{d}} is an inscribed ball.} -\item{\code{len_win} }{ The length of the sliding window for CG algorithm. The default value is \eqn{500+4dimension^2}.} -\item{\code{C} }{ A constant for the lower bound of \eqn{variance/mean^2} in schedule annealing of CG algorithm. The default value is \eqn{2}.} -\item{\code{M} }{ The number of points we sample in each step of schedule annealing in CG algorithm. The default value is \eqn{500C + dimension^2 / 2}.} -\item{\code{ratio} }{ Parameter of schedule annealing of CG algorithm, larger ratio means larger steps in schedule annealing. The default value is \eqn{1 - 1/dimension}.} -\item{\code{frac} }{ The fraction of the total error to spend in the first gaussian in CG algorithm. The default value is \eqn{0.1}.} -\item{\code{BW_rad} }{ The radius for the ball walk. The default value is \eqn{4r/dimension}, where \eqn{r} is the radius of the inscribed ball of the polytope.} -\item{\code{ub} }{ The lower bound for the ratios in MMC in CB algorithm. The default value is \eqn{0.1}.} -\item{\code{lb} }{ The upper bound for the ratios in MMC in CB algorithm. The default value is \eqn{0.15}.} -\item{\code{N} }{ An integer that controls the number of points \eqn{\nu N} generated in each convex body in annealing schedule of algorithm CB.} -\item{\code{nu} }{ The degrees of freedom for the t-student distribution in t-tests in CB algorithm. The default value is \eqn{10}.} -\item{\code{alpha} }{ The significance level for the t-tests in CB algorithm. The default values is 0.2.} -\item{\code{prob} }{ The probability is used for the empirical confidence interval in ratio estimation of CB algorithm. The default value is \eqn{0.75}.} +\item{\code{win_len} }{ The length of the sliding window for CG algorithm. The default value is \eqn{500+4dimension^2}.} \item{\code{hpoly} }{ A boolean parameter to use H-polytopes in MMC of CB algorithm. The default value is \code{FALSE}.} -\item{\code{minmaxW} }{ A boolean parameter to use the sliding window with a minmax values as a stopping criterion.} -\item{\code{L} }{The maximum length of the billiard trajectory.} }} \item{rounding}{Optional. A boolean parameter for rounding. The default value is \code{TRUE} for V-polytopes and \code{FALSE} otherwise.} diff --git a/R-proj/src/direct_sampling.cpp b/R-proj/src/direct_sampling.cpp index fb41115c7..c0c7390ea 100644 --- a/R-proj/src/direct_sampling.cpp +++ b/R-proj/src/direct_sampling.cpp @@ -15,7 +15,7 @@ #include #include #include -#include "samplers.h" +#include "new_volume.hpp" #include "simplex_samplers.h" @@ -30,6 +30,7 @@ //' \item{\code{radius} }{ The radius of the \eqn{d}-dimensional hypersphere. The default value is \eqn{1}.} //' } //' @param n The number of points that the function is going to sample. +//' @param seed Optional. A fixed seed for the number generator. //' //' @references \cite{R.Y. Rubinstein and B. Melamed, //' \dQuote{Modern simulation and modeling} \emph{ Wiley Series in Probability and Statistics,} 1998.} @@ -43,19 +44,36 @@ //' @export // [[Rcpp::export]] Rcpp::NumericMatrix direct_sampling(Rcpp::Nullable body = R_NilValue, - Rcpp::Nullable n = R_NilValue) { + Rcpp::Nullable n = R_NilValue, + Rcpp::Nullable seed = R_NilValue) { typedef double NT; typedef Cartesian Kernel; + typedef boost::mt19937 RNGType2; typedef typename Kernel::Point Point; - typedef boost::mt19937 RNGType; typedef Eigen::Matrix VT; typedef Eigen::Matrix MT; + typedef BoostRandomNumberGenerator RNGType; int dim, numpoints; NT radius = 1.0; std::list randPoints; + if (!Rcpp::as(body).containsElementNamed("dimension")) { + throw Rcpp::exception("Dimension has to be given as input!"); + } + dim = Rcpp::as(Rcpp::as(body)["dimension"]); + if (dim <=1) throw Rcpp::exception("Dimension has to be larger than 1!"); + RNGType rng(dim); + + if (seed.isNotNull()) { + unsigned seed2 = std::chrono::system_clock::now().time_since_epoch().count(); + std::cout<<"seed = "<::signaling_NaN() : Rcpp::as(seed); + //RNGType rng2(5); + numpoints = (!n.isNotNull()) ? 100 : Rcpp::as(n); if (!n.isNotNull()) { @@ -66,13 +84,8 @@ Rcpp::NumericMatrix direct_sampling(Rcpp::Nullable body = R_NilValue } if (body.isNotNull()) { - if (!Rcpp::as(body).containsElementNamed("dimension")) { - throw Rcpp::exception("Dimension has to be given as input!"); - } - dim = Rcpp::as(Rcpp::as(body)["dimension"]); - if (dim <=1) throw Rcpp::exception("Dimension has to be larger than 1!"); if (Rcpp::as(body).containsElementNamed("radius")) { radius = Rcpp::as(Rcpp::as(body)["radius"]); @@ -84,25 +97,26 @@ Rcpp::NumericMatrix direct_sampling(Rcpp::Nullable body = R_NilValue throw Rcpp::exception("The kind of body has to be given as input!"); } + //RNGType rng(dim); if (Rcpp::as(Rcpp::as(body)["type"]).compare(std::string("hypersphere"))==0) { for (unsigned int k = 0; k < numpoints; ++k) { - randPoints.push_back(get_point_on_Dsphere(dim, radius)); + randPoints.push_back(GetPointOnDsphere::apply(dim, radius, rng)); } } else if (Rcpp::as(Rcpp::as(body)["type"]).compare(std::string("ball"))==0) { for (unsigned int k = 0; k < numpoints; ++k) { - randPoints.push_back(get_point_in_Dsphere(dim, radius)); + randPoints.push_back(GetPointInDsphere::apply(dim, radius, rng)); } } else if (Rcpp::as(Rcpp::as(body)["type"]).compare(std::string("unit_simplex"))==0) { - Sam_Unit(dim, numpoints, randPoints); + Sam_Unit(dim, numpoints, randPoints, seed3); } else if (Rcpp::as(Rcpp::as(body)["type"]).compare(std::string("canonical_simplex"))==0) { - Sam_Canon_Unit(dim, numpoints, randPoints); + Sam_Canon_Unit(dim, numpoints, randPoints, seed3); } else { diff --git a/R-proj/src/exact_vol.cpp b/R-proj/src/exact_vol.cpp index eb58a9840..ce18edd8f 100644 --- a/R-proj/src/exact_vol.cpp +++ b/R-proj/src/exact_vol.cpp @@ -54,7 +54,6 @@ double exact_vol(Rcpp::Nullable P) { typedef double NT; typedef Cartesian Kernel; typedef typename Kernel::Point Point; - typedef boost::mt19937 RNGType; typedef Eigen::Matrix VT; typedef Eigen::Matrix MT; @@ -67,8 +66,6 @@ double exact_vol(Rcpp::Nullable P) { if (type == 2) { - typedef VPolytope Vpolytope; - Vpolytope VP; dim = Rcpp::as(P).field("dimension"); if (Rcpp::as(Rcpp::as(P).field("V")).rows() == diff --git a/R-proj/src/inner_ball.cpp b/R-proj/src/inner_ball.cpp index b09a946f8..66c9d012b 100644 --- a/R-proj/src/inner_ball.cpp +++ b/R-proj/src/inner_ball.cpp @@ -11,10 +11,7 @@ #include #include #include -#include "hpolytope.h" -#include "vpolytope.h" -#include "zpolytope.h" -#include "vpolyintersectvpoly.h" +#include "new_volume.hpp" //' Compute an inscribed ball of a convex polytope //' @@ -41,19 +38,17 @@ Rcpp::NumericVector inner_ball(Rcpp::Reference P) { typedef double NT; typedef Cartesian Kernel; typedef typename Kernel::Point Point; - typedef boost::mt19937 RNGType; + typedef BoostRandomNumberGenerator RNGType; typedef HPolytope Hpolytope; - typedef VPolytope Vpolytope; + typedef VPolytope Vpolytope; typedef Zonotope zonotope; - typedef IntersectionOfVpoly InterVP; + typedef IntersectionOfVpoly InterVP; typedef Eigen::Matrix VT; typedef Eigen::Matrix MT; - unsigned int n = P.field("dimension"); + unsigned int n = P.field("dimension"), type = P.field("type"); std::pair InnerBall; - int type = P.field("type"); - switch (type) { case 1: { // Hpolytope diff --git a/R-proj/src/poly_gen.cpp b/R-proj/src/poly_gen.cpp index 0a0ff9150..f51763eb4 100644 --- a/R-proj/src/poly_gen.cpp +++ b/R-proj/src/poly_gen.cpp @@ -46,7 +46,7 @@ Rcpp::NumericMatrix poly_gen (int kind_gen, bool Vpoly_gen, bool Zono_gen, int d typedef typename Kernel::Point Point; typedef boost::mt19937 RNGType; typedef HPolytope Hpolytope; - typedef VPolytope Vpolytope; + typedef VPolytope Vpolytope; typedef Zonotope zonotope; double seed2 = (!seed.isNotNull()) ? std::numeric_limits::signaling_NaN() : Rcpp::as(seed); diff --git a/R-proj/src/rotating.cpp b/R-proj/src/rotating.cpp index 1420cfdfb..daeebdcfe 100644 --- a/R-proj/src/rotating.cpp +++ b/R-proj/src/rotating.cpp @@ -15,11 +15,7 @@ #include #include #include -#include "hpolytope.h" -#include "vpolytope.h" -#include "zpolytope.h" -#include "vpolyintersectvpoly.h" -#include "samplers.h" +#include "new_volume.hpp" #include "rotating.h" #include "extractMatPoly.h" @@ -40,18 +36,17 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable Kernel; typedef typename Kernel::Point Point; - typedef boost::mt19937 RNGType; + typedef BoostRandomNumberGenerator RNGType; typedef HPolytope Hpolytope; - typedef VPolytope Vpolytope; + typedef VPolytope Vpolytope; typedef Zonotope zonotope; - typedef IntersectionOfVpoly InterVP; + typedef IntersectionOfVpoly InterVP; typedef Eigen::Matrix VT; typedef Eigen::Matrix MT; MT TransorfMat; Rcpp::NumericMatrix Mat; - unsigned int n = P.field("dimension"); - int type = P.field("type"); + unsigned int n = P.field("dimension"), type = P.field("type"); int seed2 = (!seed.isNotNull()) ? std::chrono::system_clock::now().time_since_epoch().count() : Rcpp::as(seed); diff --git a/R-proj/src/rounding.cpp b/R-proj/src/rounding.cpp index 232894535..ae54c925e 100644 --- a/R-proj/src/rounding.cpp +++ b/R-proj/src/rounding.cpp @@ -15,63 +15,53 @@ #include #include #include -#include "vars.h" -#include "hpolytope.h" -#include "vpolytope.h" -#include "zpolytope.h" -#include "samplers.h" -#include "rounding.h" -#include "vpolyintersectvpoly.h" +#include "new_volume.hpp" +#include "new_gaussian_volume.hpp" #include "extractMatPoly.h" //' Internal rcpp function for the rounding of a convex polytope //' //' @param P A convex polytope (H- or V-representation or zonotope). +//' @param seed Optional. A fixed seed for the number generator. //' //' @section warning: //' Do not use this function. //' //' @return A numerical matrix that describes the rounded polytope and contains the round value. // [[Rcpp::export]] -Rcpp::List rounding (Rcpp::Reference P){ +Rcpp::List rounding (Rcpp::Reference P, Rcpp::Nullable seed = R_NilValue){ typedef double NT; typedef Cartesian Kernel; typedef typename Kernel::Point Point; - typedef boost::mt19937 RNGType; + typedef BoostRandomNumberGenerator RNGType; typedef HPolytope Hpolytope; - typedef VPolytope Vpolytope; + typedef VPolytope Vpolytope; typedef Zonotope zonotope; - typedef IntersectionOfVpoly InterVP; typedef Eigen::Matrix VT; typedef Eigen::Matrix MT; + bool cdhr = false; + unsigned int n = P.field("dimension"), walkL, type = P.field("type"); + + RNGType rng(n); + if (seed.isNotNull()) { + unsigned seed2 = Rcpp::as(seed); + rng.set_seed(seed2); + } + Hpolytope HP; Vpolytope VP; zonotope ZP; - InterVP VPcVP; - - bool rand_only=false, - NN=false, - birk=false, - verbose=false, - cdhr=false, rdhr = false, ball_walk = false, billiard = false; - NT delta = -1.0, diam = -1.0; - - unsigned int n = P.field("dimension"); - unsigned int rnum = std::pow(1.0,-2.0) * 400 * n * std::log(n); - unsigned int walkL = 10+n/10; std::pair InnerBall; Rcpp::NumericMatrix Mat; - int type = P.field("type"); if (type == 1) { - walkL = 10 + 10/n; + walkL = 10 + 10*n; cdhr = true; } else { - walkL = 5; - billiard = true; + walkL = 2; } switch (type) { @@ -79,20 +69,17 @@ Rcpp::List rounding (Rcpp::Reference P){ // Hpolytope HP.init(n, Rcpp::as(P.field("A")), Rcpp::as(P.field("b"))); InnerBall = HP.ComputeInnerBall(); - //if (billiard && diam < 0.0) HP.comp_diam(diam, InnerBall.second); break; } case 2: { VP.init(n, Rcpp::as(P.field("V")), VT::Ones(Rcpp::as(P.field("V")).rows())); InnerBall = VP.ComputeInnerBall(); - VP.comp_diam(diam, 0.0); break; } case 3: { // Zonotope ZP.init(n, Rcpp::as(P.field("G")), VT::Ones(Rcpp::as(P.field("G")).rows())); InnerBall = ZP.ComputeInnerBall(); - ZP.comp_diam(diam, 0.0); break; } case 4: { @@ -100,30 +87,38 @@ Rcpp::List rounding (Rcpp::Reference P){ } } - unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); - // the random engine with this seed - RNGType rng(seed); - boost::random::uniform_real_distribution<>(urdist); - boost::random::uniform_real_distribution<> urdist1(-1,1); - - // initialization - vars var(rnum,n,walkL,1,0.0,0.0,0,0.0,0,InnerBall.second,diam,rng,urdist,urdist1, - delta,verbose,rand_only,false,NN,birk,ball_walk,cdhr,rdhr,billiard); std::pair< std::pair, NT > round_res; - switch (type) { case 1: { - round_res = rounding_min_ellipsoid(HP, InnerBall, var); + if (cdhr) { + std::pair , NT> res = round_polytope(HP, InnerBall, walkL, + rng); + } else { + std::pair , NT> res = round_polytope(HP, InnerBall, walkL, + rng); + } Mat = extractMatPoly(HP); break; } case 2: { - round_res = rounding_min_ellipsoid(VP, InnerBall, var); + if (cdhr) { + std::pair , NT> res = round_polytope(VP, InnerBall, walkL, + rng); + } else { + std::pair , NT> res = round_polytope(VP, InnerBall, walkL, + rng); + } Mat = extractMatPoly(VP); break; } case 3: { - round_res = rounding_min_ellipsoid(ZP, InnerBall, var); + if (cdhr) { + std::pair , NT> res = round_polytope(ZP, InnerBall, walkL, + rng); + } else { + std::pair , NT> res = round_polytope(ZP, InnerBall, walkL, + rng); + } Mat = extractMatPoly(ZP); break; } diff --git a/R-proj/src/sample_points.cpp b/R-proj/src/sample_points.cpp index e15969051..407d51f1a 100644 --- a/R-proj/src/sample_points.cpp +++ b/R-proj/src/sample_points.cpp @@ -10,19 +10,74 @@ #include #include #include -#include "cartesian_geom/cartesian_kernel.h" #include #include #include #include -#include "vars.h" -#include "hpolytope.h" -#include "vpolytope.h" -#include "zpolytope.h" -#include "samplers.h" -#include "gaussian_samplers.h" +#include "new_volume.hpp" +#include "new_gaussian_volume.hpp" #include "sample_only.h" -#include "vpolyintersectvpoly.h" + +template +void sample_from_polytope(Polytope &P, RNGType &rng, PointList &randPoints, unsigned int const& walkL, unsigned int const& numpoints, + bool const& gaussian, NT const& a, NT const& L, bool const& boundary, Point const& StartingPoint, unsigned int const& nburns, + bool const& set_L, bool const& cdhr, bool const& rdhr, bool const& billiard, bool const& ball_walk) +{ + if (boundary) { + if (cdhr) { + sampling_only_boundary (randPoints, P, rng, walkL, numpoints, + StartingPoint, nburns); + } else { + sampling_only_boundary (randPoints, P, rng, walkL, numpoints, + StartingPoint, nburns); + } + } else if (cdhr) { + if (gaussian) { + sampling_only_gaussian(randPoints, P, rng, walkL, numpoints, + a, StartingPoint, nburns); + } else { + sampling_only(randPoints, P, rng, walkL, numpoints, + StartingPoint, nburns); + } + } else if(rdhr){ + if (gaussian) { + sampling_only_gaussian(randPoints, P, rng, walkL, numpoints, + a, StartingPoint, nburns); + } else { + sampling_only(randPoints, P, rng, walkL, numpoints, + StartingPoint, nburns); + } + } else if (billiard) { + if (set_L) { + BilliardWalk WalkType(L); + sampling_only(randPoints, P, rng, WalkType, walkL, numpoints, StartingPoint, nburns); + } else { + sampling_only(randPoints, P, rng, walkL, numpoints, + StartingPoint, nburns); + } + } else { + if (set_L) { + + if (gaussian) { + GaussianBallWalk WalkType(L); + sampling_only_gaussian(randPoints, P, rng, WalkType, walkL, numpoints, a, + StartingPoint, nburns); + } else { + BallWalk WalkType(L); + sampling_only(randPoints, P, rng, WalkType, walkL, numpoints, + StartingPoint, nburns); + } + } else { + if (gaussian) { + sampling_only_gaussian(randPoints, P, rng, walkL, numpoints, + a, StartingPoint, nburns); + } else { + sampling_only(randPoints, P, rng, walkL, numpoints, + StartingPoint, nburns); + } + } + } +} //' Sample uniformly or normally distributed points from a convex Polytope (H-polytope, V-polytope or a zonotope). @@ -46,6 +101,7 @@ //' \item{\code{variance} }{ The variance of the multidimensional spherical gaussian. 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 Chebychev ball.} //' } +//' @param seed Optional. A fixed seed for the number generator. //' //' @return A \eqn{d\times n} matrix that contains, column-wise, the sampled points from the convex polytope P. //' @examples @@ -68,28 +124,38 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue, Rcpp::Nullable n = R_NilValue, Rcpp::Nullable random_walk = R_NilValue, - Rcpp::Nullable distribution = R_NilValue){ + Rcpp::Nullable distribution = R_NilValue, + Rcpp::Nullable seed = R_NilValue){ typedef double NT; typedef Cartesian Kernel; + typedef BoostRandomNumberGenerator RNGType; typedef typename Kernel::Point Point; - typedef boost::mt19937 RNGType; typedef HPolytope Hpolytope; - typedef VPolytope Vpolytope; + typedef VPolytope Vpolytope; typedef Zonotope zonotope; - typedef IntersectionOfVpoly InterVP; + typedef IntersectionOfVpoly InterVP; typedef Eigen::Matrix VT; typedef Eigen::Matrix MT; + unsigned int type = Rcpp::as(P).field("type"),dim = Rcpp::as(P).field("dimension"), + walkL = 10 + dim / 10; + + RNGType rng(dim); + if (seed.isNotNull()) { + unsigned seed2 = Rcpp::as(seed); + rng.set_seed(seed2); + } + Hpolytope HP; Vpolytope VP; zonotope ZP; InterVP VPcVP; - int type, dim, numpoints, nburns = 0; - NT radius = 1.0, delta = -1.0, diam = -1.0; + unsigned int numpoints, nburns = 0; + NT radius = 1.0, L; bool set_mode = false, cdhr = false, rdhr = false, ball_walk = false, gaussian = false, - billiard = false, boundary = false, set_starting_point = false; + billiard = false, boundary = false, set_starting_point = false, set_L = false; std::list randPoints; std::pair InnerBall; Point mode(dim); @@ -107,9 +173,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue throw Rcpp::exception("No polytope is given as input!"); } - type = Rcpp::as(P).field("type"); - dim = Rcpp::as(P).field("dimension"); - int walkL = 10 + dim / 10; + if (!distribution.isNotNull() || !Rcpp::as(distribution).containsElementNamed("density")) { billiard = true; @@ -157,7 +221,9 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue rdhr = true; } else if (Rcpp::as(Rcpp::as(random_walk)["walk"]).compare(std::string("BaW")) == 0) { if (Rcpp::as(random_walk).containsElementNamed("BaW_rad")) { - delta = Rcpp::as(Rcpp::as(random_walk)["BaW_rad"]); + L = Rcpp::as(Rcpp::as(random_walk)["BaW_rad"]); + set_L = true; + if (L<=0.0) throw Rcpp::exception("BaW diameter must be a postitive number!"); } ball_walk = true; } else if (Rcpp::as(Rcpp::as(random_walk)["walk"]).compare(std::string("BiW")) == 0) { @@ -165,7 +231,9 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue billiard = true; walkL = 5; if (Rcpp::as(random_walk).containsElementNamed("L")) { - diam = Rcpp::as(Rcpp::as(random_walk)["L"]); + L = Rcpp::as(Rcpp::as(random_walk)["L"]); + set_L = true; + 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!"); @@ -205,24 +273,13 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue } } - bool rand_only=false, - NN=false, - birk=false, - verbose=false; - - unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); - // the random engine with this seed - RNGType rng(seed); - boost::random::uniform_real_distribution<>(urdist); - boost::random::uniform_real_distribution<> urdist1(-1,1); - switch(type) { case 1: { // Hpolytope HP.init(dim, Rcpp::as(Rcpp::as(P).field("A")), Rcpp::as(Rcpp::as(P).field("b"))); - if (!set_starting_point || (!set_mode && gaussian) || ball_walk || billiard) { + if (!set_starting_point || (!set_mode && gaussian)) { InnerBall = HP.ComputeInnerBall(); if (!set_starting_point) StartingPoint = InnerBall.first; if (!set_mode && gaussian) mode = InnerBall.first; @@ -230,12 +287,10 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue if (HP.is_in(StartingPoint) == 0) { throw Rcpp::exception("The given point is not in the interior of the polytope!"); } - if (billiard && diam < 0.0) HP.comp_diam(diam, InnerBall.second); HP.normalize(); if (gaussian) { StartingPoint = StartingPoint - mode; HP.shift(mode.getCoefficients()); - } break; } @@ -244,14 +299,13 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue VP.init(dim, Rcpp::as(Rcpp::as(P).field("V")), VT::Ones(Rcpp::as(Rcpp::as(P).field("V")).rows())); - if (!set_starting_point || (!set_mode && gaussian) || ball_walk) { + if (!set_starting_point || (!set_mode && gaussian)) { InnerBall = VP.ComputeInnerBall(); if (!set_starting_point) StartingPoint = InnerBall.first; if (!set_mode && gaussian) mode = InnerBall.first; } if (VP.is_in(StartingPoint) == 0) throw Rcpp::exception("The given point is not in the interior of the polytope!"); - if (billiard && diam < 0.0) VP.comp_diam(diam); if (gaussian) { StartingPoint = StartingPoint - mode; VP.shift(mode.getCoefficients()); @@ -263,14 +317,13 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue ZP.init(dim, Rcpp::as(Rcpp::as(P).field("G")), VT::Ones(Rcpp::as(Rcpp::as(P).field("G")).rows())); - if (!set_starting_point || (!set_mode && gaussian) || ball_walk) { + if (!set_starting_point || (!set_mode && gaussian)) { InnerBall = ZP.ComputeInnerBall(); if (!set_starting_point) StartingPoint = InnerBall.first; if (!set_mode && gaussian) mode = InnerBall.first; } if (ZP.is_in(StartingPoint) == 0) throw Rcpp::exception("The given point is not in the interior of the polytope!"); - if (billiard && diam < 0.0) ZP.comp_diam(diam); if (gaussian) { StartingPoint = StartingPoint - mode; ZP.shift(mode.getCoefficients()); @@ -293,9 +346,6 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue if (!set_mode && gaussian) mode = InnerBall.first; if (VPcVP.is_in(StartingPoint) == 0) throw Rcpp::exception("The given point is not in the interior of the polytope!"); - if (billiard && diam < 0.0) { - VPcVP.comp_diam(diam, InnerBall.second); - } if (gaussian) { StartingPoint = StartingPoint - mode; VPcVP.shift(mode.getCoefficients()); @@ -304,38 +354,25 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue } } - if (ball_walk && delta < 0.0) { - if (gaussian) { - delta = 4.0 * InnerBall.second / std::sqrt(std::max(NT(1.0), a) * NT(dim)); - } else { - delta = 4.0 * InnerBall.second / std::sqrt(NT(dim)); - } - } - - vars var1(1,dim,walkL,1,0.0,0.0,0,0.0,0,InnerBall.second,diam,rng,urdist,urdist1, - delta,verbose,rand_only,false,NN,birk,ball_walk,cdhr,rdhr, billiard); - vars_g var2(dim, walkL, 0, 0, 1, 0, InnerBall.second, rng, 0, 0, 0, delta, verbose, - rand_only, false, NN, birk, ball_walk, cdhr, rdhr); - switch (type) { case 1: { - sampling_only(randPoints, HP, walkL, numpoints, gaussian, - a, boundary, StartingPoint, nburns, var1, var2); + sample_from_polytope(HP, rng, randPoints, walkL, numpoints, gaussian, a, L, boundary, StartingPoint, nburns, + set_L, cdhr, rdhr, billiard, ball_walk); break; } case 2: { - sampling_only(randPoints, VP, walkL, numpoints, gaussian, - a, boundary, StartingPoint, nburns, var1, var2); + sample_from_polytope(VP, rng, randPoints, walkL, numpoints, gaussian, a, L, boundary, StartingPoint, nburns, + set_L, cdhr, rdhr, billiard, ball_walk); break; } case 3: { - sampling_only(randPoints, ZP, walkL, numpoints, gaussian, - a, boundary, StartingPoint, nburns, var1, var2); + sample_from_polytope(ZP, rng, randPoints, walkL, numpoints, gaussian, a, L, boundary, StartingPoint, nburns, + set_L, cdhr, rdhr, billiard, ball_walk); break; } case 4: { - sampling_only(randPoints, VPcVP, walkL, numpoints, gaussian, - a, boundary, StartingPoint, nburns, var1, var2); + sample_from_polytope(VPcVP, rng, randPoints, walkL, numpoints, gaussian, a, L, boundary, StartingPoint, nburns, + set_L, cdhr, rdhr, billiard, ball_walk); break; } } diff --git a/R-proj/src/volume.cpp b/R-proj/src/volume.cpp index b2050b5be..62820b505 100644 --- a/R-proj/src/volume.cpp +++ b/R-proj/src/volume.cpp @@ -13,105 +13,69 @@ #include #include #include -#include "volume.h" -#include "cooling_balls.h" -#include "cooling_hpoly.h" +#include "new_volume.hpp" +#include "new_gaussian_volume.hpp" +#include "new_cooling_balls.hpp" +#include "new_cooling_hpoly.hpp" -template -double generic_volume(Polytope& P, unsigned int walk_step, double e, - std::pair InnerBall, bool CG, bool CB, bool hpoly, unsigned int win_len, - unsigned int N, double C, double ratio, double frac, NT lb, NT ub, NT p, NT alpha, - unsigned int NN, unsigned int nu, bool win2, bool ball_walk, double delta, bool cdhr, - bool rdhr, bool billiard, double diam, bool rounding, int type) +template +double generic_volume(Polytope& P, RNGType &rng, unsigned int walk_length, NT e, + bool CG, bool CB, unsigned int win_len, + bool rounding, bool cdhr, bool rdhr, bool ball_walk, + bool billiard, int type) { - bool rand_only=false, - NNN=false, - birk=false, - verbose = false; - unsigned int n_threads=1; - NT round_val = 1.0, rmax = 0.0; + typedef typename Polytope::MT MT; + typedef typename Polytope::VT VT; + typedef typename Polytope::PointType Point; - //unsigned int m;//=A.nrow()-1; - unsigned int n = P.dimension();//=A.ncol()-1; - unsigned int rnum = std::pow(e,-2) * 400 * n * std::log(n); + typedef HPolytope Hpolytope; - unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); - // the random engine with this seed - typedef boost::mt19937 RNGType; - RNGType rng(seed); - boost::random::uniform_real_distribution<>(urdist); - boost::random::uniform_real_distribution<> urdist1(-1,1); + NT round_val = 1.0; + unsigned int n = P.dimension(); - std::pair InnerB; + if (rounding) { + std::pair InnerBall = P.ComputeInnerBall(); - if(InnerBall.second > 0.0) { //if it is given as an input - - InnerB = InnerBall; - - } else if(type == 2 && CB) { - if (rounding) { - - //std::cout<<"hello2"< var2(1, n, 1, n_threads, 0.0, e, 0, 0.0, 0, InnerB.second, 2 * P.get_max_vert_norm(), - rng, urdist, urdist1, -1, verbose, rand_only, rounding, NNN, birk, ball_walk, - cdhr, rdhr, billiard); - std::pair res_round = rounding_min_ellipsoid(P, InnerB, var2); - round_val = res_round.first; - - rounding = false; - InnerB.second = 0.0; - InnerB.first = Point(n); - get_vpoly_center(P); - rmax = P.get_max_vert_norm(); + if (type == 1) { + round_val = round_polytope(P, InnerBall, 10 + 10 * n, rng).second; } else { - InnerB.second = 0.0; - InnerB.first = Point(n); - get_vpoly_center(P); - rmax = P.get_max_vert_norm(); + round_val = round_polytope(P, InnerBall, 2, rng).second; } - } else { - // no internal ball or point is given as input - InnerB = P.ComputeInnerBall(); } - //set parameters for billiard and ball walk - if (billiard && diam < 0.0) { - P.comp_diam(diam, InnerB.second); - } else if (ball_walk && delta < 0.0) { - delta = 4.0 * InnerB.second / NT(n); - } - - // initialization - vars var(rnum,n,walk_step,n_threads,0.0,e,0,0.0,0, InnerB.second, diam, rng,urdist,urdist1, - delta,verbose,rand_only,rounding,NNN,birk,ball_walk,cdhr,rdhr, billiard); NT vol; if (CG) { - vars var2(rnum, n, 10 + n / 10, n_threads, 0.0, e, 0, 0.0, 0, InnerB.second, diam, rng, - urdist, urdist1, delta, verbose, rand_only, rounding, NNN, birk, ball_walk, cdhr, - rdhr, billiard); - vars_g var1(n, walk_step, N, win_len, 1, e, InnerB.second, rng, C, frac, ratio, delta, verbose, - rand_only, rounding, NN, birk, ball_walk, cdhr, rdhr); - vol = volume_gaussian_annealing(P, var1, var2, InnerB); + if (cdhr) { + vol = volume_gaussian_annealing(P, rng, e, walk_length); + } else if (rdhr) { + vol = volume_gaussian_annealing(P, rng, e, walk_length); + } else { + vol = volume_gaussian_annealing(P, rng, e, walk_length); + } } else if (CB) { - vars_ban var_ban(lb, ub, p, rmax, alpha, win_len, NN, nu, win2); - if (!hpoly) { - vol = vol_cooling_balls(P, var, var_ban, InnerB); + if (cdhr) { + vol = volume_cooling_balls(P, rng, e, walk_length, win_len); + } else if (rdhr) { + vol = volume_cooling_balls(P, rng, e, walk_length, win_len); + } else if (ball_walk) { + vol = volume_cooling_balls(P, rng, e, walk_length, win_len); } else { - vars_g varg(n, 1, N, 6 * n * n + 500, 1, e, InnerB.second, rng, C, frac, ratio, delta, - verbose, rand_only, false, false, birk, false, true, false); - vol = vol_cooling_hpoly < HPolytope < Point > > (P, var, var_ban, varg, InnerB); + vol = volume_cooling_balls(P, rng, e, walk_length, win_len); } - if (vol < 0.0) { - throw Rcpp::exception("Simulated annealing failed! Try to increase the walk length."); + } else { + if (cdhr) { + vol = volume_sequence_of_balls(P, rng, e, walk_length); + } else if (rdhr) { + vol = volume_sequence_of_balls(P, rng, e, walk_length); + } else if (ball_walk) { + vol = volume_sequence_of_balls(P, rng, e, walk_length); + } else { + vol = volume_sequence_of_balls(P, rng, e, walk_length); } - }else { - vol = volume(P, var, InnerB); } - - return vol * round_val; + vol *= round_val; + return vol; } //' The main function for volume approximation of a convex Polytope (H-polytope, V-polytope or a zonotope) @@ -125,24 +89,11 @@ double generic_volume(Polytope& P, unsigned int walk_step, double e, //' \item{\code{error} }{ A numeric value to set the upper bound for the approximation error. The default value is \eqn{1} for \code{'SOB'} and \eqn{0.1} otherwise.} //' \item{\code{random_walk} }{ A string that declares the random walk method: a) \code{'CDHR'} for Coordinate Directions Hit-and-Run, b) \code{'RDHR'} for Random Directions Hit-and-Run, c) \code{'BaW'} for Ball Walk, or \code{'BiW'} for Billiard walk. The default walk is \code{'CDHR'} for H-polytopes and \code{'BiW'} for the other representations.} //' \item{\code{walk_length} }{ An integer to set the number of the steps for the random walk. The default value is \eqn{\lfloor 10 + d/10\rfloor} for \code{'SOB'} and \eqn{1} otherwise.} -//' \item{\code{inner_ball} }{ A \eqn{d+1} numeric vector that contains an inner ball. The first \eqn{d} coordinates corresponds to the center and the last one to the radius of the ball. If it is not given then for H-polytopes the Chebychev ball is computed, for V-polytopes \eqn{d+1} vertices are picked randomly and the Chebychev ball of the defined simplex is computed. For a zonotope that is defined by the Minkowski sum of \eqn{m} segments we compute the maximal \eqn{r} s.t.: \eqn{re_i\in Z} for all \eqn{i=1,\dots ,d}, then the ball centered at the origin with radius \eqn{r/\sqrt{d}} is an inscribed ball.} -//' \item{\code{len_win} }{ The length of the sliding window for CG algorithm. The default value is \eqn{500+4dimension^2}.} -//' \item{\code{C} }{ A constant for the lower bound of \eqn{variance/mean^2} in schedule annealing of CG algorithm. The default value is \eqn{2}.} -//' \item{\code{M} }{ The number of points we sample in each step of schedule annealing in CG algorithm. The default value is \eqn{500C + dimension^2 / 2}.} -//' \item{\code{ratio} }{ Parameter of schedule annealing of CG algorithm, larger ratio means larger steps in schedule annealing. The default value is \eqn{1 - 1/dimension}.} -//' \item{\code{frac} }{ The fraction of the total error to spend in the first gaussian in CG algorithm. The default value is \eqn{0.1}.} -//' \item{\code{BW_rad} }{ The radius for the ball walk. The default value is \eqn{4r/dimension}, where \eqn{r} is the radius of the inscribed ball of the polytope.} -//' \item{\code{ub} }{ The lower bound for the ratios in MMC in CB algorithm. The default value is \eqn{0.1}.} -//' \item{\code{lb} }{ The upper bound for the ratios in MMC in CB algorithm. The default value is \eqn{0.15}.} -//' \item{\code{N} }{ An integer that controls the number of points \eqn{\nu N} generated in each convex body in annealing schedule of algorithm CB.} -//' \item{\code{nu} }{ The degrees of freedom for the t-student distribution in t-tests in CB algorithm. The default value is \eqn{10}.} -//' \item{\code{alpha} }{ The significance level for the t-tests in CB algorithm. The default values is 0.2.} -//' \item{\code{prob} }{ The probability is used for the empirical confidence interval in ratio estimation of CB algorithm. The default value is \eqn{0.75}.} +//' \item{\code{win_len} }{ The length of the sliding window for CG algorithm. The default value is \eqn{500+4dimension^2}.} //' \item{\code{hpoly} }{ A boolean parameter to use H-polytopes in MMC of CB algorithm. The default value is \code{FALSE}.} -//' \item{\code{minmaxW} }{ A boolean parameter to use the sliding window with a minmax values as a stopping criterion.} -//' \item{\code{L} }{The maximum length of the billiard trajectory.} //' } //' @param rounding Optional. A boolean parameter for rounding. The default value is \code{TRUE} for V-polytopes and \code{FALSE} otherwise. +//' @param seed Optional. A fixed seed for the number generator. //' //' @references \cite{I.Z.Emiris and V. Fisikopoulos, //' \dQuote{Practical polytope volume approximation,} \emph{ACM Trans. Math. Soft.,} 2014.}, @@ -168,27 +119,32 @@ double generic_volume(Polytope& P, unsigned int walk_step, double e, // [[Rcpp::export]] double volume (Rcpp::Reference P, Rcpp::Nullable settings = R_NilValue, - Rcpp::Nullable rounding = R_NilValue) { + Rcpp::Nullable rounding = R_NilValue, + Rcpp::Nullable seed = R_NilValue) { typedef double NT; typedef Cartesian Kernel; typedef typename Kernel::Point Point; - typedef boost::mt19937 RNGType; + typedef BoostRandomNumberGenerator RNGType; typedef HPolytope Hpolytope; - typedef VPolytope Vpolytope; + typedef VPolytope Vpolytope; typedef Zonotope zonotope; - typedef IntersectionOfVpoly InterVP; + typedef IntersectionOfVpoly InterVP; typedef Eigen::Matrix VT; typedef Eigen::Matrix MT; - unsigned int n = P.field("dimension"), walkL; - int type = P.field("type"); + unsigned int n = P.field("dimension"), walkL, type = P.field("type"); - bool CG = false, CB = false, cdhr = false, rdhr = false, ball_walk = false, round = false, win2 = false, - hpoly = false, billiard = false, set_mean_point = false; - unsigned int win_len = 4*n*n+500, N = 500 * 2 + n * n / 2, NN = 120 + (n*n)/10, nu = 10, cg_params = 0, cb_params = 0; + RNGType rng(n); + if (seed.isNotNull()) { + unsigned seed2 = Rcpp::as(seed); + rng.set_seed(seed2); + } - NT C = 2.0, ratio = 1.0-1.0/(NT(n)), frac = 0.1, e, delta = -1.0, lb = 0.1, ub = 0.15, p = 0.75, rmax = 0.0, - alpha = 0.2, diam = -1.0; + bool CG = false, CB = false, cdhr = false, rdhr = false, ball_walk = false, round = false, + hpoly = false, billiard = false; + unsigned int win_len = 4*n*n+500; + + NT e; if (!Rcpp::as(settings).containsElementNamed("algorithm")) { if (type == 2 || type == 3) { @@ -237,8 +193,7 @@ double volume (Rcpp::Reference P, } else { if (CB) { billiard = true; - win_len = 170; - NN = 125; + win_len = 200; } else { rdhr = true; } @@ -262,8 +217,7 @@ double volume (Rcpp::Reference P, } } else { billiard = true; - win_len = 170; - NN = 125; + win_len = 200; } }else { throw Rcpp::exception("Unknown walk type!"); @@ -279,84 +233,9 @@ double volume (Rcpp::Reference P, round = (!rounding.isNotNull()) ? false : Rcpp::as(rounding); } - if (Rcpp::as(settings).containsElementNamed("BW_rad")) { - delta = Rcpp::as(Rcpp::as(settings)["BW_rad"]); - } - if (Rcpp::as(settings).containsElementNamed("C")) { - C = Rcpp::as(Rcpp::as(settings)["C"]); - N = 500 * ((int) C) + n * n / 2; - cg_params++; - } - if (Rcpp::as(settings).containsElementNamed("M")) { - N = Rcpp::as(Rcpp::as(settings)["M"]); - cg_params++; - } - if (Rcpp::as(settings).containsElementNamed("len_win")) { - win_len = Rcpp::as(Rcpp::as(settings)["len_win"]); - if (!CB && !CG) Rf_warning("flag 'len_win' can be used only for CG or CB algorithms."); - } - if (Rcpp::as(settings).containsElementNamed("frac")) { - frac = Rcpp::as(Rcpp::as(settings)["frac"]); - cg_params++; - } - if (Rcpp::as(settings).containsElementNamed("ratio")) { - ratio = Rcpp::as(Rcpp::as(settings)["ratio"]); - cg_params++; - } - if (Rcpp::as(settings).containsElementNamed("hpoly")) { - hpoly = Rcpp::as(Rcpp::as(settings)["hpoly"]); - cb_params++; - if ((hpoly && !CB) || (type != 3 && CB && hpoly)) - Rf_warning("flag 'hpoly' can be used to only in MMC of CB algorithm for zonotopes."); - } - if (Rcpp::as(settings).containsElementNamed("lb")) { - lb = Rcpp::as(Rcpp::as(settings)["lb"]); - cb_params++; - } - if (Rcpp::as(settings).containsElementNamed("ub")) { - ub = Rcpp::as(Rcpp::as(settings)["ub"]); - cb_params++; - } - if (Rcpp::as(settings).containsElementNamed("nu")) { - nu = Rcpp::as(Rcpp::as(settings)["nu"]); - cb_params++; - } - if (Rcpp::as(settings).containsElementNamed("N")) { - NN = Rcpp::as(Rcpp::as(settings)["N"]); - cb_params++; - } - if (Rcpp::as(settings).containsElementNamed("minmaxW")) { - win2 = Rcpp::as(Rcpp::as(settings)["minmaxW"]); - cb_params++; - } - if (Rcpp::as(settings).containsElementNamed("prob")) { - p = Rcpp::as(Rcpp::as(settings)["prob"]); - cb_params++; - } - if (Rcpp::as(settings).containsElementNamed("alpha")) { - alpha = Rcpp::as(Rcpp::as(settings)["alpha"]); - cb_params++; - } - if (Rcpp::as(settings).containsElementNamed("L")) { - diam = Rcpp::as(Rcpp::as(settings)["L"]); - } - - if ((CB && cg_params > 0) || (CG && cb_params > 0) || (!CB & !CG && (cg_params > 0 || cb_params > 0))){ - Rf_warning("Maybe some algorithm's input parameters are not going to be used."); - } - - std::pair inner_ball; - inner_ball.second = -1.0; - if (Rcpp::as(settings).containsElementNamed("inner_ball")) { - if (Rcpp::as(Rcpp::as(settings)["inner_ball"]).size() != n+1) { - throw Rcpp::exception("Inscribed ball has to lie in the same dimension as the polytope P"); - } else { - set_mean_point = true; - VT temp = Rcpp::as(Rcpp::as(settings)["inner_ball"]); - inner_ball.first = Point(n, std::vector(&temp[0], temp.data()+temp.cols()*temp.rows()-1)); - inner_ball.second = temp(n); - if (inner_ball.second <= 0.0) throw Rcpp::exception("The radius of the given inscribed ball has to be a positive number."); - } + if (Rcpp::as(settings).containsElementNamed("win_len")) { + win_len = Rcpp::as(Rcpp::as(settings)["win_len"]); + if (!CB && !CG) Rf_warning("flag 'win_len' can be used only for CG or CB algorithms."); } switch(type) { @@ -364,31 +243,42 @@ double volume (Rcpp::Reference P, // Hpolytope Hpolytope HP; HP.init(n, Rcpp::as(P.field("A")), Rcpp::as(P.field("b"))); - if (set_mean_point){ - if (HP.is_in(inner_ball.first) == 0) throw Rcpp::exception("The center of the given inscribed ball is not in the interior of the polytope!"); - } - return generic_volume(HP, walkL, e, inner_ball, CG, CB, hpoly, win_len, N, C, ratio, frac, lb, ub, p, - alpha, NN, nu, win2, ball_walk, delta, cdhr, rdhr, billiard, diam, round, type); + return generic_volume(HP, rng, walkL, e, CG, CB, win_len, round, + cdhr, rdhr, ball_walk, billiard, type); } case 2: { // Vpolytope Vpolytope VP; VP.init(n, Rcpp::as(P.field("V")), VT::Ones(Rcpp::as(P.field("V")).rows())); - if (set_mean_point){ - if (VP.is_in(inner_ball.first) == 0) throw Rcpp::exception("The center of the given inscribed ball is not in the interior of the polytope!"); - } - return generic_volume(VP, walkL, e, inner_ball, CG, CB, hpoly, win_len, N, C, ratio, frac, lb, ub, p, - alpha, NN, nu, win2, ball_walk, delta, cdhr, rdhr, billiard, diam, round, type); + return generic_volume(VP, rng, walkL, e, CG, CB, win_len, round, + cdhr, rdhr, ball_walk, billiard, type); } case 3: { // Zonotope zonotope ZP; ZP.init(n, Rcpp::as(P.field("G")), VT::Ones(Rcpp::as(P.field("G")).rows())); - if (set_mean_point){ - if (ZP.is_in(inner_ball.first) == 0) throw Rcpp::exception("The center of the given inscribed ball is not in the interior of the polytope!"); + if (Rcpp::as(settings).containsElementNamed("hpoly")) { + hpoly = Rcpp::as(Rcpp::as(settings)["hpoly"]); + if (hpoly && !CB) + Rf_warning("flag 'hpoly' can be used to only in MMC of CB algorithm for zonotopes."); + } else if (ZP.num_of_generators() / ZP.dimension() < 5 ) { + hpoly = true; + } else { + hpoly = false; } - return generic_volume(ZP, walkL, e, inner_ball, CG, CB, hpoly, win_len, N, C, ratio, frac, lb, ub, p, - alpha, NN, nu, win2, ball_walk, delta, cdhr, rdhr, billiard, diam, round, type); + if (hpoly) { + if (cdhr) { + return volume_cooling_hpoly(ZP, rng, e, walkL, win_len); + } else if (rdhr) { + return volume_cooling_hpoly(ZP, rng, e, walkL, win_len); + } else if (ball_walk) { + return volume_cooling_hpoly(ZP, rng, e, walkL, win_len); + } else { + return volume_cooling_hpoly(ZP, rng, e, walkL, win_len); + } + } + return generic_volume(ZP, rng, walkL, e, CG, CB, win_len, round, + cdhr, rdhr, ball_walk, billiard, type); } case 4: { // Intersection of two V-polytopes @@ -399,11 +289,8 @@ double volume (Rcpp::Reference P, VP2.init(n, Rcpp::as(P.field("V2")), VT::Ones(Rcpp::as(P.field("V2")).rows())); VPcVP.init(VP1, VP2); if (!VPcVP.is_feasible()) throw Rcpp::exception("Empty set!"); - if (set_mean_point){ - if (VPcVP.is_in(inner_ball.first) == 0) throw Rcpp::exception("The center of the given inscribed ball is not in the interior of the polytope!"); - } - return generic_volume(VPcVP, walkL, e, inner_ball, CG, CB, hpoly, win_len, N, C, ratio, frac, lb, ub, p, - alpha, NN, nu, win2, ball_walk, delta, cdhr, rdhr, billiard, diam, round, type); + return generic_volume(VPcVP, rng, walkL, e, CG, CB, win_len, round, + cdhr, rdhr, ball_walk, billiard, type); } } diff --git a/R-proj/src/zonotope_approximation.cpp b/R-proj/src/zonotope_approximation.cpp index 24f6217db..896cb10d6 100644 --- a/R-proj/src/zonotope_approximation.cpp +++ b/R-proj/src/zonotope_approximation.cpp @@ -11,15 +11,17 @@ #include #include #include -#include "volume.h" -#include "cooling_balls.h" -#include "cooling_hpoly.h" +#include "new_volume.hpp" +#include "new_gaussian_volume.hpp" +#include "new_cooling_balls.hpp" +#include "new_cooling_hpoly.hpp" //' An internal Rccp function for the over-approximation of a zonotope //' //' @param Z A zonotope. //' @param fit_ratio Optional. A boolean parameter to request the computation of the ratio of fitness. //' @param settings Optional. A list that declares the values of the parameters of CB algorithm. +//' @param seed Optional. A fixed seed for the number generator. //' //' @section warning: //' Do not use this function. @@ -28,25 +30,27 @@ // [[Rcpp::export]] Rcpp::List zono_approx (Rcpp::Reference Z, Rcpp::Nullable fit_ratio = R_NilValue, - Rcpp::Nullable settings = R_NilValue) { + Rcpp::Nullable settings = R_NilValue, + Rcpp::Nullable seed = R_NilValue) { typedef double NT; typedef Cartesian Kernel; + typedef BoostRandomNumberGenerator RNGType; typedef typename Kernel::Point Point; - typedef boost::mt19937 RNGType; typedef HPolytope Hpolytope; typedef Zonotope zonotope; typedef Eigen::Matrix VT; typedef Eigen::Matrix MT; + int n = Rcpp::as(Z.field("dimension")), k = Rcpp::as(Z.field("G")).rows(), win_len = 200, walkL = 1; - int n = Rcpp::as(Z.field("dimension")), k = Rcpp::as(Z.field("G")).rows(); - double e = 0.1, delta = -1.0, lb = 0.1, ub = 0.15, p = 0.75, rmax = 0.0, alpha = 0.2, diam = -1.0; - int win_len = 3 * n * n + 400, NN = 220 + (n * n) / 10, nu = 10, walkL = 1; - bool ball_walk = false, verbose = false, cdhr = false, rdhr = false, billiard = false, round = false, win2 = false, - hpoly = false, set_mean_point = false; - + RNGType rng(n); + if (seed.isNotNull()) { + unsigned seed2 = Rcpp::as(seed); + rng.set_seed(seed2); + } - NT ratio = std::numeric_limits::signaling_NaN(); + NT e = 0.1, ratio = std::numeric_limits::signaling_NaN();; + bool hpoly = false; MT X(n, 2 * k); X << Rcpp::as(Z.field("G")).transpose(), -Rcpp::as(Z.field("G")).transpose(); @@ -57,7 +61,6 @@ Rcpp::List zono_approx (Rcpp::Reference Z, MT A(n, 2 * n); A << -MT::Identity(n, n), MT::Identity(n, n); MT Mat(2 * n, n + 1); - Mat << Gred_ii, A.transpose() * svd.matrixU().transpose(); Hpolytope HP; @@ -73,109 +76,25 @@ Rcpp::List zono_approx (Rcpp::Reference Z, Rcpp::as(settings)["walk_length"]); e = (!Rcpp::as(settings).containsElementNamed("error")) ? 0.1 : Rcpp::as( Rcpp::as(settings)["error"]); - round = (!Rcpp::as(settings).containsElementNamed("rounding")) ? false : Rcpp::as( - Rcpp::as(settings)["rounding"]); - //if (rounding.isNotNull()) round = Rcpp::as(rounding); - if (!Rcpp::as(settings).containsElementNamed("random_walk")) { - billiard = true; - win_len = 170; - NN = 125; - } else if (Rcpp::as(Rcpp::as(settings)["random_walk"]).compare(std::string("BiW")) == 0) { - billiard = true; - win_len = 170; - NN = 125; - } else if (Rcpp::as(Rcpp::as(settings)["random_walk"]).compare(std::string("RDHR")) == 0) { - rdhr = true; - } else if (Rcpp::as(Rcpp::as(settings)["random_walk"]).compare(std::string("CDHR")) == 0) { - cdhr = true; - } else if (Rcpp::as(Rcpp::as(settings)["random_walk"]).compare(std::string("Baw")) == 0) { - ball_walk = true; - } else { - throw Rcpp::exception("Unknown walk type!"); - } - - - if (Rcpp::as(settings).containsElementNamed("BW_rad")) { - delta = Rcpp::as(Rcpp::as(settings)["BW_rad"]); - } - if (Rcpp::as(settings).containsElementNamed("len_win")) { - win_len = Rcpp::as(Rcpp::as(settings)["len_win"]); - } - if (Rcpp::as(settings).containsElementNamed("lb")) { - lb = Rcpp::as(Rcpp::as(settings)["lb"]); - } - if (Rcpp::as(settings).containsElementNamed("ub")) { - ub = Rcpp::as(Rcpp::as(settings)["ub"]); - } - if (Rcpp::as(settings).containsElementNamed("nu")) { - nu = Rcpp::as(Rcpp::as(settings)["nu"]); - } - if (Rcpp::as(settings).containsElementNamed("N")) { - NN = Rcpp::as(Rcpp::as(settings)["N"]); - } - if (Rcpp::as(settings).containsElementNamed("minmaxW")) { - win2 = Rcpp::as(Rcpp::as(settings)["minmaxW"]); - } - if (Rcpp::as(settings).containsElementNamed("prob")) { - p = Rcpp::as(Rcpp::as(settings)["prob"]); - } - if (Rcpp::as(settings).containsElementNamed("alpha")) { - alpha = Rcpp::as(Rcpp::as(settings)["alpha"]); - } - if (Rcpp::as(settings).containsElementNamed("hpoly")) { - hpoly = Rcpp::as(Rcpp::as(settings)["hpoly"]); - } - if (Rcpp::as(settings).containsElementNamed("L")) { - diam = Rcpp::as(Rcpp::as(settings)["L"]); - } + win_len = (!Rcpp::as(settings).containsElementNamed("win_len")) ? 200 : Rcpp::as( + Rcpp::as(settings)["win_len"]); zonotope ZP; ZP.init(n, Rcpp::as(Z.field("G")), VT::Ones(Rcpp::as(Z.field("G")).rows())); - unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); - // the random engine with this seed - typedef boost::mt19937 RNGType; - RNGType rng(seed); - boost::random::uniform_real_distribution<>(urdist); - boost::random::uniform_real_distribution<> urdist1(-1, 1); - - - std::pair InnerB; - InnerB.second = -1.0; - if (Rcpp::as(settings).containsElementNamed("inner_ball")) { - if (Rcpp::as(Rcpp::as(settings)["inner_ball"]).size() != n + 1) { - throw Rcpp::exception("Inscribed ball has to lie in the same dimension as the polytope P"); - } else { - set_mean_point = true; - VT temp = Rcpp::as(Rcpp::as(settings)["inner_ball"]); - InnerB.first = Point(n, std::vector(&temp[0], temp.data() + temp.cols() * temp.rows() - 1)); - InnerB.second = temp(n); - if (InnerB.second <= 0.0) - throw Rcpp::exception("The radius of the given inscribed ball has to be a positive number."); - } + if (Rcpp::as(settings).containsElementNamed("hpoly")) { + hpoly = Rcpp::as(Rcpp::as(settings)["hpoly"]); + } else if (ZP.num_of_generators() / ZP.dimension() < 5 ) { + hpoly = true; } else { - InnerB = ZP.ComputeInnerBall(); + hpoly = false; } - if (billiard && diam < 0.0) { - ZP.comp_diam(diam, 0.0); - } else if (ball_walk && delta < 0.0) { - delta = 4.0 * InnerB.second / NT(n); - } - - vars var(1, n, walkL, 1, 0.0, e, 0, 0.0, 0, InnerB.second, diam, rng, - urdist, urdist1, delta, false, false, round, false, false, ball_walk, cdhr, rdhr, - billiard); - vars_ban var_ban(lb, ub, p, 0.0, alpha, win_len, NN, nu, win2); - NT vol; if (!hpoly) { - vol = vol_cooling_balls(ZP, var, var_ban, InnerB); + vol = volume_cooling_balls(ZP, rng, e, walkL, win_len); } else { - vars_g varg(n, 1, 1000 + n * n / 2, 6 * n * n + 500, 1, e, InnerB.second, rng, 2.0, 0.1, - 1.0 - 1.0 / (NT(n)), delta, false, false, false, false, false, false, true, - false); - vol = vol_cooling_hpoly(ZP, var, var_ban, varg, InnerB); + vol = volume_cooling_hpoly(ZP, rng, e, walkL, win_len); } ratio = std::pow(vol_red / vol, 1.0 / NT(n)); } diff --git a/include/convex_bodies/ballintersectconvex.h b/include/convex_bodies/ballintersectconvex.h index f0863b43c..11cc45d30 100644 --- a/include/convex_bodies/ballintersectconvex.h +++ b/include/convex_bodies/ballintersectconvex.h @@ -90,7 +90,7 @@ class BallIntersectPolytope { comp_diam(diam); } - NT radius() { + NT radius() const { return B.radius(); } diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index 1bfd53d12..1e3e1f847 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -87,6 +87,10 @@ class HPolytope{ return A.rows(); } + int num_of_generators() const { + return 0; + } + // return the matrix A MT get_mat() const { diff --git a/include/convex_bodies/vpolyintersectvpoly.h b/include/convex_bodies/vpolyintersectvpoly.h index a06e50a89..5f3f17cea 100644 --- a/include/convex_bodies/vpolyintersectvpoly.h +++ b/include/convex_bodies/vpolyintersectvpoly.h @@ -83,6 +83,10 @@ class IntersectionOfVpoly { return NT(0); } + int num_of_generators() const { + return 0; + } + //std::vector get_vertices() const { // return vecV; //} @@ -310,7 +314,7 @@ class IntersectionOfVpoly { P2.linear_transformIt(T); } - std::vector get_dists(const NT &radius) { + std::vector get_dists(const NT &radius) const { std::vector res(upper_bound_of_hyperplanes(), radius); return res; } diff --git a/include/convex_bodies/vpolytope.h b/include/convex_bodies/vpolytope.h index 6d709721d..a8fafbeaf 100644 --- a/include/convex_bodies/vpolytope.h +++ b/include/convex_bodies/vpolytope.h @@ -75,6 +75,10 @@ class VPolytope{ return 0; } + int num_of_generators() const { + return 0; + } + // compute the number of facets of the cyclic polytope in dimension _d with the same number of vertices // this is an upper bound for the number of the facets from McMullen's Upper Bound Theorem unsigned int upper_bound_of_hyperplanes() const { diff --git a/include/generators/h_polytopes_gen.h b/include/generators/h_polytopes_gen.h index c40f43ed4..837e77c02 100644 --- a/include/generators/h_polytopes_gen.h +++ b/include/generators/h_polytopes_gen.h @@ -17,7 +17,7 @@ Polytope random_hpoly(unsigned int dim, unsigned int m, double seed = std::numer typedef typename Polytope::MT MT; typedef typename Polytope::VT VT; typedef typename Polytope::NT NT; - typedef typename Polytope::PolytopePoint Point; + typedef typename Polytope::PointType Point; unsigned rng_seed = std::chrono::system_clock::now().time_since_epoch().count(); RNGType rng(rng_seed); diff --git a/include/samplers/sample_only.h b/include/samplers/sample_only.h index 23034fc1d..da9190e48 100644 --- a/include/samplers/sample_only.h +++ b/include/samplers/sample_only.h @@ -22,25 +22,163 @@ #ifndef SAMPLE_ONLY_H #define SAMPLE_ONLY_H -template -void sampling_only(PointList &randPoints, Polytope &P, const unsigned int walk_len, - const unsigned int rnum, bool gaussian, const NT &a, const bool boundary, - const Point &starting_point, int nburns, UParameters const& var1, GParameters const& var2) { +template +void sampling_only(PointList &randPoints, Polytope &P, RandomNumberGenerator &rng, const unsigned int &walk_len, + const unsigned int &rnum, const Point &starting_point, unsigned int const& nburns) { + + typedef typename WalkTypePolicy::template Walk + < + Polytope, + RandomNumberGenerator + > walk; + + //RandomNumberGenerator rng(P.dimension()); + PushBackWalkPolicy push_back_policy; - typedef typename UParameters::RNGType RNGType; Point p = starting_point; - if (nburns > 0 ) rand_point_generator(P, p, 1, nburns, randPoints, var1); + typedef RandomPointGenerator RandomPointGenerator; + 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); + + +} +template < + typename RandomNumberGenerator, + typename PointList, + typename Polytope, + typename WalkTypePolicy, + typename Point +> +void sampling_only(PointList &randPoints, Polytope &P, RandomNumberGenerator &rng, WalkTypePolicy &WalkType, const unsigned int &walk_len, + const unsigned int &rnum, + const Point &starting_point, unsigned int const& nburns) +{ + typedef typename WalkTypePolicy::template Walk + < + Polytope, + RandomNumberGenerator + > walk; + + //RandomNumberGenerator rng(P.dimension()); + PushBackWalkPolicy push_back_policy; + typedef RandomPointGenerator RandomPointGenerator; + + Point p = starting_point; + + RandomPointGenerator::apply(P, p, nburns, walk_len, randPoints, + push_back_policy, rng, WalkType.param); randPoints.clear(); - if (boundary) { - boundary_rand_point_generator(P, p, rnum/2, walk_len, randPoints, var1); - } else if (!gaussian){ - rand_point_generator(P, p, rnum, walk_len, randPoints, var1); - } else { - rand_gaussian_point_generator(P, p, rnum, walk_len, randPoints, a, var2); - } + RandomPointGenerator::apply(P, p, rnum, walk_len, randPoints, + push_back_policy, rng, WalkType.param); +} + +//---------------------------------------------------------------------------------- + +template +void sampling_only_boundary(PointList &randPoints, Polytope &P, RandomNumberGenerator &rng, const unsigned int &walk_len, + const unsigned int &rnum, const Point &starting_point, unsigned int const& nburns) { + typedef typename WalkTypePolicy::template Walk + < + Polytope, + RandomNumberGenerator + > walk; + + //RandomNumberGenerator rng(P.dimension()); + PushBackWalkPolicy push_back_policy; + + Point p = starting_point; + + typedef BoundaryRandomPointGenerator BoundaryRandomPointGenerator; + BoundaryRandomPointGenerator::apply(P, p, nburns, walk_len, + randPoints, push_back_policy, rng); + randPoints.clear(); + BoundaryRandomPointGenerator::apply(P, p, rnum / 2, walk_len, + randPoints, push_back_policy, rng); } + +template < + typename WalkTypePolicy, + typename RandomNumberGenerator, + typename PointList, + typename Polytope, + typename NT, + typename Point + > +void sampling_only_gaussian(PointList &randPoints, Polytope &P, RandomNumberGenerator &rng, const unsigned int &walk_len, + const unsigned int &rnum, const NT &a, + const Point &starting_point, unsigned int const& nburns) { + + typedef typename WalkTypePolicy::template Walk + < + Polytope, + RandomNumberGenerator + > walk; + + //RandomNumberGenerator rng(P.dimension()); + PushBackWalkPolicy push_back_policy; + + Point p = starting_point; + + typedef GaussianRandomPointGenerator RandomPointGenerator; + 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); + + +} + + +template < + typename RandomNumberGenerator, + typename PointList, + typename Polytope, + typename WalkTypePolicy, + typename NT, + typename Point + > +void sampling_only_gaussian(PointList &randPoints, Polytope &P, RandomNumberGenerator &rng, WalkTypePolicy &WalkType, const unsigned int &walk_len, + const unsigned int &rnum, const NT &a, + const Point &starting_point, unsigned int const& nburns) { + + typedef typename WalkTypePolicy::template Walk + < + Polytope, + RandomNumberGenerator + > walk; + + //RandomNumberGenerator rng(P.dimension()); + PushBackWalkPolicy push_back_policy; + + Point p = starting_point; + + typedef GaussianRandomPointGenerator RandomPointGenerator; + 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); + + +} + + #endif diff --git a/include/samplers/simplex_samplers.h b/include/samplers/simplex_samplers.h index 28774d15c..e5f953967 100644 --- a/include/samplers/simplex_samplers.h +++ b/include/samplers/simplex_samplers.h @@ -22,15 +22,23 @@ #define SIMPLEX_SAMPLERS_H template -void Sam_Unit(unsigned int dim, unsigned int num, std::list &points){ +void Sam_Unit(unsigned int dim, unsigned int num, std::list &points, + double seed = std::numeric_limits::signaling_NaN()){ unsigned int j,i,x_rand,M=2147483647,pr,divisors,pointer; // M is the largest possible integer std::vector x_vec; std::vector y; boost::random::uniform_int_distribution<> uidist(1,M); - unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); - RNGType rng(seed); + unsigned rng_seed = std::chrono::system_clock::now().time_since_epoch().count(); + RNGType rng(rng_seed); + if (!std::isnan(seed)) { + std::cout<<"seed = "< &points){ } template -void Sam_Canon_Unit(unsigned int dim, unsigned int num, std::list &points){ +void Sam_Canon_Unit(unsigned int dim, unsigned int num, std::list &points, + double seed = std::numeric_limits::signaling_NaN()){ unsigned int j,i,x_rand,M=2147483647,pointer; // M is the largest possible integer //std::vector x_vec; @@ -162,8 +171,15 @@ void Sam_Canon_Unit(unsigned int dim, unsigned int num, std::list &points dim--; boost::random::uniform_int_distribution<> uidist(1,M); - unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); - RNGType rng(seed); + unsigned rng_seed = std::chrono::system_clock::now().time_since_epoch().count(); + RNGType rng(rng_seed); + if (!std::isnan(seed)) { + std::cout<<"seed = "< x_vec2; NT Ti,sum; diff --git a/include/volume/new_basic_sampling_features.hpp b/include/volume/new_basic_sampling_features.hpp index 76f96fe33..5556bc475 100644 --- a/include/volume/new_basic_sampling_features.hpp +++ b/include/volume/new_basic_sampling_features.hpp @@ -46,6 +46,10 @@ struct BoostRandomNumberGenerator return _ndist(_rng); } + void set_seed(unsigned rng_seed){ + _rng.seed(rng_seed); + } + private : RNGType _rng; boost::random::uniform_real_distribution _urdist; @@ -79,6 +83,10 @@ struct BoostRandomNumberGenerator return _ndist(_rng); } + void set_seed(unsigned rng_seed){ + _rng.seed(rng_seed); + } + private : RNGType _rng; boost::random::uniform_real_distribution _urdist; @@ -178,6 +186,20 @@ struct GetPointInDsphere } }; +template +struct GetPointOnDsphere +{ + template + inline static Point apply(unsigned int const& dim, + NT const& radius, + RandomNumberGenerator &rng) + { + Point p = GetDirection::apply(dim, rng); + if (radius != 0) p *= radius; + return p; + } +}; + ////////////////////////////// Random Point Generators /// @@ -205,11 +227,35 @@ struct RandomPointGenerator WalkPolicy &policy, RandomNumberGenerator &rng, Parameters const& parameters) + { + Walk walk(P, p, rng, parameters); + for (unsigned int i=0; i + static void apply(Polytope& P, + Point &p, // a point to start + unsigned int const& rnum, + unsigned int const& walk_length, + PointList &randPoints, + WalkPolicy &policy, + RandomNumberGenerator &rng) { Walk walk(P, p, rng); for (unsigned int i=0; i struct cooling_ball_parameters { - cooling_ball_parameters() + cooling_ball_parameters(unsigned int win_len) : lb(0.1) , ub(0.15) , p(0.75) , rmax(0) , alpha(0.2) - , win_len(1000) - , N(150) + , win_len(win_len) + , N(125) , nu(10) , window2(false) {} @@ -368,7 +368,6 @@ bool get_sequence_of_polytopeballs(Polytope& P, template bool is_max_error(NT const& a, NT const& b, NT const& error) { - std::cout<<"(b-a)/a = "<<(b-a)/a<<", error = "< -NT estimate_ratio(PolyBall1& Pb1, +NT estimate_ratio(PolyBall1 const& Pb1, PolyBall2 const& Pb2, NT const& ratio, NT const& error, @@ -627,7 +626,7 @@ template typename NT, typename RNG > -NT estimate_ratio_interval(PolyBall1& Pb1, +NT estimate_ratio_interval(PolyBall1 const& Pb1, PolyBall2 const& Pb2, NT const& ratio, NT const& error, @@ -664,13 +663,15 @@ NT estimate_ratio_interval(PolyBall1& Pb1, template < - typename WalkTypePolicy = BilliardWalk, - typename RandomNumberGenerator = BoostRandomNumberGenerator, + typename WalkTypePolicy, + typename RandomNumberGenerator, typename Polytope > double volume_cooling_balls(Polytope const& Pin, + RandomNumberGenerator &rng, double const& error = 1.0, - unsigned int const& walk_length = 1) + unsigned int const& walk_length = 1, + unsigned int const& win_len = 250) { typedef typename Polytope::PointType Point; typedef typename Point::FT NT; @@ -687,8 +688,8 @@ double volume_cooling_balls(Polytope const& Pin, typedef RandomPointGenerator RandomPointGenerator; auto P(Pin); - RandomNumberGenerator rng(P.dimension()); - cooling_ball_parameters parameters; + //RandomNumberGenerator rng(P.dimension()); + cooling_ball_parameters parameters(win_len); int n = P.dimension(); NT prob = parameters.p; @@ -777,4 +778,21 @@ double volume_cooling_balls(Polytope const& Pin, return vol; } + + +template + < + typename WalkTypePolicy, + typename RandomNumberGenerator, + typename Polytope + > +double volume_cooling_balls(Polytope const& Pin, + double const& error = 0.1, + unsigned int const& walk_length = 1) +{ + RandomNumberGenerator rng(Pin.dimension()); + return volume_cooling_balls(Pin, rng, error, walk_length); +} + + #endif diff --git a/include/volume/new_cooling_hpoly.hpp b/include/volume/new_cooling_hpoly.hpp index 4b38363ee..8300ab014 100644 --- a/include/volume/new_cooling_hpoly.hpp +++ b/include/volume/new_cooling_hpoly.hpp @@ -221,15 +221,16 @@ void compute_hpoly_for_mmc(Zonotope &P, HPolytope &HP) { template < - typename WalkTypePolicy = BilliardWalk, - typename RandomNumberGenerator = BoostRandomNumberGenerator, + typename WalkTypePolicy, typename HPolytope, - typename Zonotope - + typename Zonotope, + typename RandomNumberGenerator > double volume_cooling_hpoly (Zonotope const& Pin, + RandomNumberGenerator &rng, double const& error = 1.0, - unsigned int const& walk_length = 1) + unsigned int const& walk_length = 1, + unsigned int const& win_len = 200) { typedef typename Zonotope::PointType Point; @@ -254,9 +255,9 @@ double volume_cooling_hpoly (Zonotope const& Pin, typedef RandomPointGenerator CdhrRandomPointGenerator; auto P(Pin); - RandomNumberGenerator rng(P.dimension()); - RandomNumberGenerator rng_diam(P.num_of_generators()); - cooling_ball_parameters parameters; + //RandomNumberGenerator rng(P.dimension()); + //RandomNumberGenerator rng_diam(P.num_of_generators()); + cooling_ball_parameters parameters(win_len); int n = P.dimension(); NT prob = parameters.p, ratio; @@ -289,7 +290,7 @@ double volume_cooling_hpoly (Zonotope const& Pin, HPolytope HP2(HP); std::pair InnerBall = HP2.ComputeInnerBall(); - std::pair< std::pair, NT > res = round_polytope(HP2, InnerBall, + std::pair< std::pair, NT > res = round_polytope(HP2, InnerBall, 10 + 10 * n, rng); //TODO: rounding to HP2 //NT vol = res.second * volume_cooling_balls(HP2, Her, 1); @@ -351,4 +352,20 @@ double volume_cooling_hpoly (Zonotope const& Pin, } + +template + < + typename WalkTypePolicy, + typename RandomNumberGenerator, + typename HPolytope, + typename Polytope + > +double volume_cooling_hpoly(Polytope const& Pin, + double const& error = 0.1, + unsigned int const& walk_length = 1) +{ + RandomNumberGenerator rng(Pin.dimension()); + return volume_cooling_hpoly(Pin, rng, error, walk_length); +} + #endif diff --git a/include/volume/new_gaussian_volume.hpp b/include/volume/new_gaussian_volume.hpp index c23e27e01..20bf4d674 100644 --- a/include/volume/new_gaussian_volume.hpp +++ b/include/volume/new_gaussian_volume.hpp @@ -129,6 +129,25 @@ NT chord_random_point_generator_exp_coord(const NT &l, struct GaussianBallWalk { + GaussianBallWalk(double L) + : param(L, true) + {} + + GaussianBallWalk() + : param(0, false) + {} + + struct parameters + { + parameters(double L, bool set) + : m_L(L), set_delta(set) + {} + double m_L; + bool set_delta; + }; + + parameters param; + template < typename Polytope, @@ -139,11 +158,19 @@ struct Walk typedef typename Polytope::PointType Point; typedef typename Point::FT NT; - Walk (Polytope const& P, Point&, NT const&, RandomNumberGenerator&) + Walk (Polytope const& P, Point &p, NT const& a, RandomNumberGenerator &rng) { _delta = ((P.InnerBall()).second * NT(4)) / NT(P.dimension()); } + Walk (Polytope const& P, Point &p, NT const& a, RandomNumberGenerator &rng, parameters const& params) { + if (params.set_delta) { + _delta = params.m_L; + } else { + _delta = ((P.InnerBall()).second * NT(4)) / NT(P.dimension()); + } + } + template inline void apply(BallPolytope const& P, Point &p, // a point to start @@ -184,6 +211,9 @@ private : struct GaussianRDHRWalk { + struct parameters {}; + parameters param; + template < typename Polytope, @@ -197,6 +227,9 @@ struct Walk Walk(Polytope const&, Point &, NT const&, RandomNumberGenerator &) {} + Walk(Polytope const&, Point &, NT const&, RandomNumberGenerator &, parameters &) + {} + template < typename BallPolytope @@ -228,6 +261,9 @@ struct Walk struct GaussianCDHRWalk { + struct parameters {}; + parameters param; + template < typename Polytope, @@ -245,6 +281,11 @@ struct Walk initialize(P, p, a_i, rng); } + Walk(Polytope const& P, Point & p, NT const& a_i, RandomNumberGenerator &rng, parameters &) + { + initialize(P, p, a_i, rng); + } + template < typename BallPolytope @@ -344,6 +385,17 @@ struct update_delta> } }; +template +struct update_delta> +{ + template + static void apply(BilliardWalk::Walk walk, + NT L) + { + walk.update_delta(L); + } +}; + ////////////////////////////// Random Point Generators /// @@ -373,9 +425,41 @@ struct GaussianRandomPointGenerator RandomNumberGenerator &rng) { Walk walk(P, p, a_i, rng); - update_delta - ::apply(walk, 4.0 * P.InnerBall().second - / std::sqrt(std::max(NT(1.0), a_i) * NT(P.dimension()))); + //update_delta + // ::apply(walk, 4.0 * P.InnerBall().second + // / std::sqrt(std::max(NT(1.0), a_i) * NT(P.dimension()))); + + for (unsigned int i=0; i + static void apply(Polytope const& P, + Point &p, // a point to start + NT const& a_i, + unsigned int const& rnum, + unsigned int const& walk_length, + PointList &randPoints, + WalkPolicy &policy, + RandomNumberGenerator &rng, + Parameters const& parameters) + { + Walk walk(P, p, a_i, rng, parameters); + //update_delta + //::apply(walk, 4.0 * P.InnerBall().second + // / std::sqrt(std::max(NT(1.0), a_i) * NT(P.dimension()))); for (unsigned int i=0; i, + typename WalkTypePolicy, + typename RandomNumberGenerator, typename Polytope > double volume_gaussian_annealing(Polytope const& Pin, + RandomNumberGenerator &rng, double const& error = 0.1, unsigned int const& walk_length = 1) { @@ -669,7 +754,7 @@ double volume_gaussian_annealing(Polytope const& Pin, unsigned int n = P.dimension(); unsigned int m = P.num_of_hyperplanes(); gaussian_annealing_parameters parameters(P.dimension()); - RandomNumberGenerator rng(n); + //RandomNumberGenerator rng(n); // Consider Chebychev center as an internal point auto InnerBall = P.ComputeInnerBall(); @@ -812,4 +897,19 @@ double volume_gaussian_annealing(Polytope const& Pin, } +template + < + typename WalkTypePolicy, + typename RandomNumberGenerator, + typename Polytope + > +double volume_gaussian_annealing(Polytope const& Pin, + double const& error = 0.1, + unsigned int const& walk_length = 1) +{ + RandomNumberGenerator rng(Pin.dimension()); + return volume_gaussian_annealing(Pin, rng, error, walk_length); +} + + #endif diff --git a/include/volume/new_rounding.hpp b/include/volume/new_rounding.hpp index eab14ede6..f548da544 100644 --- a/include/volume/new_rounding.hpp +++ b/include/volume/new_rounding.hpp @@ -13,16 +13,15 @@ template < typename WalkTypePolicy, - typename RandomNumberGenerator = BoostRandomNumberGenerator, typename MT, typename VT, typename Polytope, typename Point, typename NT, - typename RNG + typename RandomNumberGenerator > std::pair< std::pair, NT > round_polytope(Polytope &P, std::pair &InnerBall, - const unsigned int &walk_length, RNG &rng) + const unsigned int &walk_length, RandomNumberGenerator &rng) { typedef typename WalkTypePolicy::template Walk < diff --git a/include/volume/new_volume.hpp b/include/volume/new_volume.hpp index a1a09d677..d57ad17c6 100644 --- a/include/volume/new_volume.hpp +++ b/include/volume/new_volume.hpp @@ -17,10 +17,10 @@ #include #include -#include "random.hpp" -#include "random/uniform_int.hpp" -#include "random/normal_distribution.hpp" -#include "random/uniform_real_distribution.hpp" +//#include "random.hpp" +//#include "random/uniform_int.hpp" +//#include "random/normal_distribution.hpp" +//#include "random/uniform_real_distribution.hpp" #include "cartesian_geom/cartesian_kernel.h" //#include "vars.h" @@ -33,29 +33,35 @@ #include "zonoIntersecthpoly.h" #include "vpolyintersectvpoly.h" #include "new_rounding.hpp" -#include "samplers.h" -#include "rounding.h" -#include "gaussian_samplers.h" -#include "gaussian_annealing.h" +//#include "samplers.h" +//#include "rounding.h" +//#include "gaussian_samplers.h" +//#include "gaussian_annealing.h" #include "khach.h" + /////////////////// Random Walks // ball walk with uniform target distribution struct BallWalk { BallWalk(double L) - : param(L) + : param(L, true) + {} + + BallWalk() + : param(0, false) {} struct parameters { - parameters(double L) - : m_L(L) + parameters(double L, bool set) + : m_L(L), set_delta(set) {} double m_L; + bool set_delta; }; parameters param; @@ -75,30 +81,70 @@ struct BallWalk typedef Zonotope zonotope; typedef ZonoIntersectHPoly ZonoHPoly; - Walk (Polytope const&, Point&, RandomNumberGenerator&) {} - Walk (BallPolytope const&, Point &, RandomNumberGenerator &) {} + Walk (Polytope const& P, Point&, RandomNumberGenerator&) { + _delta = ((P.InnerBall()).second * NT(4)) / NT(P.dimension()); + } + + Walk (BallPolytope const& P, Point &, RandomNumberGenerator &) { + _delta = ((P.InnerBall()).second * NT(4)) / NT(P.dimension()); + } + Walk (BallType const&, Point &, RandomNumberGenerator &) {} - Walk(ZonoHPoly const& P, Point & p, RandomNumberGenerator &) {} - template + Walk(ZonoHPoly const& P, Point & p, RandomNumberGenerator &) { + _delta = ((P.InnerBall()).second * NT(4)) / NT(P.dimension()); + } + + Walk (Polytope const& P, Point&, RandomNumberGenerator&, parameters const& params) { + if (params.set_delta) { + _delta = params.m_L; + } else { + _delta = ((P.InnerBall()).second * NT(4)) / NT(P.dimension()); + } + } + Walk (BallPolytope const& P, Point &, RandomNumberGenerator &, parameters const& params) { + if (params.set_delta) { + _delta = params.m_L; + } else { + _delta = ((P.InnerBall()).second * NT(4)) / NT(P.dimension()); + } + } + Walk (BallType const&, Point &, RandomNumberGenerator &, parameters &) {} + + Walk(ZonoHPoly const& P, Point & p, RandomNumberGenerator &, parameters const& params) { + if (params.set_delta) { + _delta = params.m_L; + } else { + _delta = ((P.InnerBall()).second * NT(4)) / NT(P.dimension()); + } + } + + template inline void apply(BallPolytope const& P, Point &p, // a point to start unsigned int const& walk_length, - RandomNumberGenerator &rng, - Parameters const& parameters) + RandomNumberGenerator &rng) { - const NT delta = ((P.InnerBall()).second * NT(4)) / NT(P.dimension()); + //const NT delta = ((P.InnerBall()).second * NT(4)) / NT(P.dimension()); for (auto j=0u; j::apply(P.dimension(), - delta, + _delta, rng); y += p; if (P.is_in(y) == -1) p = y; } //std::cout << "use" << parameters.m_L << std::endl; } + + inline void update_delta(NT delta) + { + _delta = delta; + } + + private: + double _delta; }; }; @@ -107,6 +153,8 @@ struct BallWalk // random directions hit-and-run walk with uniform target distribution struct RDHRWalk { + struct parameters {}; + parameters param; template < @@ -123,23 +171,40 @@ struct Walk typedef Zonotope zonotope; typedef ZonoIntersectHPoly ZonoHPoly; - Walk(Polytope const& P, Point & p, RandomNumberGenerator &rng) + Walk(Polytope const& P, Point &p, RandomNumberGenerator &rng) { initialize(P, p, rng); } - Walk(BallPolytope const& P, Point & p, RandomNumberGenerator &rng) + Walk(BallPolytope const& P, Point &p, RandomNumberGenerator &rng) { initialize(P, p, rng); } - Walk(ZonoHPoly const& P, Point & p, RandomNumberGenerator &rng) + Walk(ZonoHPoly const& P, Point &p, RandomNumberGenerator &rng) { initialize(P, p, rng); } Walk (BallType const&, Point &, RandomNumberGenerator &) {} + Walk(Polytope const& P, Point &p, RandomNumberGenerator &rng, parameters &) + { + initialize(P, p, rng); + } + + Walk(BallPolytope const& P, Point &p, RandomNumberGenerator &rng, parameters &) + { + initialize(P, p, rng); + } + + Walk(ZonoHPoly const& P, Point &p, RandomNumberGenerator &rng, parameters &) + { + initialize(P, p, rng); + } + + Walk (BallType const&, Point &, RandomNumberGenerator &, parameters &) {} + template < typename BallPolytope @@ -188,6 +253,8 @@ private : // random directions hit-and-run walk with uniform target distribution struct CDHRWalk { + struct parameters {}; + parameters param; template < @@ -204,23 +271,40 @@ struct Walk typedef BallIntersectPolytope BallPolytope; typedef ZonoIntersectHPoly ZonoHPoly; - Walk(Polytope const& P, Point & p, RandomNumberGenerator &rng) + Walk(Polytope const& P, Point &p, RandomNumberGenerator &rng) { initialize(P, p, rng); } - Walk(BallPolytope const& P, Point & p, RandomNumberGenerator &rng) + Walk(BallPolytope const& P, Point &p, RandomNumberGenerator &rng) { initialize(P, p, rng); } - Walk(ZonoHPoly const& P, Point & p, RandomNumberGenerator &rng) + Walk(ZonoHPoly const& P, Point &p, RandomNumberGenerator &rng) { initialize(P, p, rng); } Walk (BallType const&, Point &, RandomNumberGenerator &) {} + Walk(Polytope const& P, Point & p, RandomNumberGenerator &rng, parameters const& params) + { + initialize(P, p, rng); + } + + Walk(BallPolytope const& P, Point & p, RandomNumberGenerator &rng, parameters const& params) + { + initialize(P, p, rng); + } + + Walk(ZonoHPoly const& P, Point & p, RandomNumberGenerator &rng, parameters const& params) + { + initialize(P, p, rng); + } + + Walk (BallType const&, Point &, RandomNumberGenerator &, parameters &) {} + template < typename BallPolytope @@ -447,10 +531,10 @@ template struct compute_diameter> { template -static NT compute(HPolytope &P) +static NT compute(HPolytope const& P) { NT diameter = NT(4) * std::sqrt(NT(P.dimension())) * P.InnerBall().second; - P.set_diameter(diameter); + //P.set_diameter(diameter); return diameter; } }; @@ -459,7 +543,7 @@ template struct compute_diameter> { template -static NT compute(VPolytope &P) +static NT compute(VPolytope const& P) { typedef typename VPolytope::MT MT; NT diameter = NT(0), diam_iter; @@ -472,7 +556,7 @@ static NT compute(VPolytope &P) } } } - P.set_diameter(diameter); + //P.set_diameter(diameter); return diameter; } }; @@ -481,7 +565,7 @@ template struct compute_diameter> { template -static NT compute(Zonotope &P) +static NT compute(Zonotope const& P) { typedef typename Zonotope::MT MT; typedef typename Zonotope::VT VT; @@ -507,7 +591,7 @@ static NT compute(Zonotope &P) for (int j = 0; j < k; ++j) x0(j) = (obj_fun(j) < 0.0) ? -1.0 : 1.0; NT diameter = NT(2) * (V.transpose() * x0).norm(); - P.set_diameter(diameter); + //P.set_diameter(diameter); return diameter; } }; @@ -516,10 +600,10 @@ template struct compute_diameter, RandomNumberGenerator>> { template -static NT compute(IntersectionOfVpoly, RandomNumberGenerator> &P) +static NT compute(IntersectionOfVpoly, RandomNumberGenerator> const& P) { NT diameter = NT(2) * NT(P.dimension()) * P.InnerBall().second; - P.set_diameter(diameter); + //P.set_diameter(diameter); return diameter; } }; @@ -528,10 +612,10 @@ template struct compute_diameter>> { template -static NT compute(BallIntersectPolytope> &P) +static NT compute(BallIntersectPolytope> const& P) { NT diameter = NT(2) * P.radius(); - P.set_diameter(diameter); + //P.set_diameter(diameter); return diameter; } }; @@ -540,7 +624,7 @@ template struct compute_diameter, HPolytope>> { template -static NT compute(ZonoIntersectHPoly, HPolytope> &P) +static NT compute(ZonoIntersectHPoly, HPolytope> const& P) { typedef typename ZonoIntersectHPoly, HPolytope>::VT VT; typedef typename ZonoIntersectHPoly, HPolytope>::MT MT; @@ -585,7 +669,7 @@ static NT compute(ZonoIntersectHPoly, HPolytope> &P) if (iter_norm > max_norm) max_norm = iter_norm; } NT diameter = NT(2) * max_norm; - P.set_diameter(diameter); + //P.set_diameter(diameter); return diameter; } }; @@ -594,6 +678,25 @@ static NT compute(ZonoIntersectHPoly, HPolytope> &P) // billiard walk for uniform distribution struct BilliardWalk { + BilliardWalk(double L) + : param(L, true) + {} + + BilliardWalk() + : 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 < @@ -610,23 +713,68 @@ struct Walk typedef Ball BallType; typedef BallIntersectPolytope BallPolytope; - Walk(Polytope& P, Point & p, RandomNumberGenerator &rng) + Walk(Polytope const& P, Point &p, RandomNumberGenerator &rng) { + _L = compute_diameter::template compute(P); + initialize(P, p, rng); + } + + Walk(BallPolytope const& P, Point &p, RandomNumberGenerator &rng) + { + _L = compute_diameter::template compute(P); + initialize(P, p, rng); + } + + Walk(ZonoHPoly const& P, Point & p, RandomNumberGenerator &rng) + { + _L = compute_diameter::template compute(P); + initialize(P, p, rng); + } + + Walk(Polytope const& P, Point & p, RandomNumberGenerator &rng, parameters const& params) + { + if(params.set_L) + { + _L = params.m_L; + } + else + { + _L = compute_diameter::template compute(P); + } initialize(P, p, rng); } - Walk(BallPolytope& P, Point & p, RandomNumberGenerator &rng) + Walk(BallPolytope const& P, Point & p, RandomNumberGenerator &rng, parameters const& params) { + if(params.set_L) + { + _L = params.m_L; + } + else + { + _L = compute_diameter::template compute(P); + } initialize(P, p, rng); } - Walk(ZonoHPoly& P, Point & p, RandomNumberGenerator &rng) + Walk(ZonoHPoly const& P, Point & p, RandomNumberGenerator &rng, parameters const& params) { + if(params.set_L) + { + _L = params.m_L; + } + else + { + _L = compute_diameter::template compute(P); + } initialize(P, p, rng); } + Walk (BallType const&, Point &, RandomNumberGenerator &, parameters &) {} + Walk (BallType const&, Point &, RandomNumberGenerator &) {} + template < typename GenericPolytope @@ -637,13 +785,13 @@ struct Walk RandomNumberGenerator &rng) { unsigned int n = P.dimension(); - NT diameter = P.get_diameter(); - NT T = rng.sample_urdist() * diameter; + //NT diameter = P.get_diameter(); + NT T = rng.sample_urdist() * _L; const NT dl = 0.995; for (auto j=0u; j::apply(n, rng); Point p0 = _p; int it = 0; @@ -662,31 +810,34 @@ struct Walk P.compute_reflection(_v, _p, pbpair.second); it++; } - if (it == 10*n) _p = p0; + if (it == 30*n) _p = p0; } p = _p; } + inline void update_delta(NT L) + { + _L = L; + } + private : template < typename GenericPolytope > - inline void initialize(GenericPolytope& P, + inline void initialize(GenericPolytope const& P, Point &p, RandomNumberGenerator &rng) { unsigned int n = P.dimension(); const NT dl = 0.995; - NT diameter = compute_diameter::template compute(P); - _lambdas.setZero(P.num_of_hyperplanes()); _Av.setZero(P.num_of_hyperplanes()); _p = p; _v = GetDirection::apply(n, rng); - NT T = rng.sample_urdist() * diameter; + NT T = rng.sample_urdist() * _L; Point p0 = _p; int it = 0; @@ -702,7 +853,7 @@ private : T -= _lambda_prev; P.compute_reflection(_v, _p, pbpair.second); - while (it < 10*n) + while (it < 30*n) { std::pair pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev); @@ -720,6 +871,7 @@ private : if (it == 10*n) _p = p0; } + double _L; Point _p; Point _v; NT _lambda_prev; @@ -737,13 +889,12 @@ private : template < - //typename WalkTypePolicy = CDHRWalk, - typename RandomNumberGenerator = BoostRandomNumberGenerator, - typename Polytope, - typename WalkType + typename WalkTypePolicy, + typename RandomNumberGenerator, + typename Polytope > double volume_sequence_of_balls(Polytope const& Pin, - WalkType walk, + RandomNumberGenerator &rng, double const& error = 1.0, unsigned int const& walk_length = 1, unsigned int const& n_threads = 1) @@ -753,18 +904,18 @@ double volume_sequence_of_balls(Polytope const& Pin, typedef Ball Ball; typedef BallIntersectPolytope BallPoly; - typedef typename WalkType::template Walk + typedef typename WalkTypePolicy::template Walk < Polytope, RandomNumberGenerator - > WalkTypeInner; + > walk; - typedef RandomPointGenerator RandomPointGenerator; + typedef RandomPointGenerator RandomPointGenerator; auto P(Pin); //copy and work with P because we are going to shift unsigned int n = P.dimension(); unsigned int rnum = std::pow(error, -2) * 400 * n * std::log(n); - RandomNumberGenerator rng(P.dimension()); + //RandomNumberGenerator rng(P.dimension()); //Compute the Chebychev ball (largest inscribed ball) with center and radius auto InnerBall = P.ComputeInnerBall(); @@ -790,14 +941,14 @@ double volume_sequence_of_balls(Polytope const& Pin, std::list randPoints; //ds for storing rand points PushBackWalkPolicy push_back_policy; - RandomPointGenerator::apply(P, p, 1, 50*n, randPoints, push_back_policy, rng, walk.param); + RandomPointGenerator::apply(P, p, 1, 50*n, randPoints, push_back_policy, rng); #ifdef VOLESTI_DEBUG double tstart2 = (double)clock()/(double)CLOCKS_PER_SEC; std::cout<<"\nCompute "< counting_policy(nump_PBSmall, PBSmall); RandomPointGenerator::apply(PBLarge, p_gen, rnum-nump_PBLarge, walk_length, randPoints, - counting_policy, rng, walk.param); + counting_policy, rng); nump_PBSmall = counting_policy.get_nump_PBSmall(); @@ -920,4 +1071,21 @@ double volume_sequence_of_balls(Polytope const& Pin, return vol; } + +template + < + typename WalkTypePolicy, + typename RandomNumberGenerator, + typename Polytope + > +double volume_sequence_of_balls(Polytope const& Pin, + double const& error = 1.0, + unsigned int const& walk_length = 1, + unsigned int const& n_threads = 1) +{ + RandomNumberGenerator rng(Pin.dimension()); + return volume_sequence_of_balls(Pin, rng, error, walk_length, n_threads); +} + + #endif diff --git a/test/benchmarks_cb.cpp b/test/benchmarks_cb.cpp index 492015a6a..85bd39e51 100644 --- a/test/benchmarks_cb.cpp +++ b/test/benchmarks_cb.cpp @@ -77,23 +77,23 @@ int main() tstart = (double)clock()/(double)CLOCKS_PER_SEC; std::cout << "BallWalk (cube) = " - << volume_cooling_balls(HP, e, walk_len) << " , "; + << volume_cooling_balls(HP, e, walk_len) << " , "; std::cout << (double)clock()/(double)CLOCKS_PER_SEC - tstart << std::endl; tstart = (double)clock()/(double)CLOCKS_PER_SEC; tstart = (double)clock()/(double)CLOCKS_PER_SEC; std::cout << "RDHRWalk (cube) = " - << volume_cooling_balls(HP, e, walk_len) << " , "; + << volume_cooling_balls(HP, e, walk_len) << " , "; std::cout << (double)clock()/(double)CLOCKS_PER_SEC - tstart << std::endl; tstart = (double)clock()/(double)CLOCKS_PER_SEC; std::cout << "CDHRWalk (cube) = " - << volume_cooling_balls(HP, e, walk_len) << " , "; + << volume_cooling_balls(HP, e, walk_len) << " , "; std::cout << (double)clock()/(double)CLOCKS_PER_SEC - tstart << std::endl; tstart = (double)clock()/(double)CLOCKS_PER_SEC; std::cout << "BirdWalk (cube) = " - << volume_cooling_balls(HP, e, walk_len) << " , "; + << volume_cooling_balls(HP, e, walk_len) << " , "; std::cout << (double)clock()/(double)CLOCKS_PER_SEC - tstart << std::endl; // OLD Implementation @@ -170,25 +170,25 @@ int main() VP.init(VP.dimension(), VP.get_mat(), VP.get_vec()); tstart = (double)clock()/(double)CLOCKS_PER_SEC; std::cout << "Ball (cross) = " - << volume_cooling_balls(VP) << " , "; + << volume_cooling_balls(VP) << " , "; std::cout << (double)clock()/(double)CLOCKS_PER_SEC - tstart << std::endl; VP.init(VP.dimension(), VP.get_mat(), VP.get_vec()); tstart = (double)clock()/(double)CLOCKS_PER_SEC; std::cout << "RDHR (cross) = " - << volume_cooling_balls(VP) << " , "; + << volume_cooling_balls(VP) << " , "; std::cout << (double)clock()/(double)CLOCKS_PER_SEC - tstart << std::endl; VP.init(VP.dimension(), VP.get_mat(), VP.get_vec()); tstart = (double)clock()/(double)CLOCKS_PER_SEC; std::cout << "CDHR (cross) = " - << volume_cooling_balls(VP) << " , "; + << volume_cooling_balls(VP) << " , "; std::cout << (double)clock()/(double)CLOCKS_PER_SEC - tstart << std::endl; VP.init(VP.dimension(), VP.get_mat(), VP.get_vec()); tstart = (double)clock()/(double)CLOCKS_PER_SEC; std::cout << "Blrd (cross) = " - << volume_cooling_balls(VP) << " , "; + << volume_cooling_balls(VP) << " , "; std::cout << (double)clock()/(double)CLOCKS_PER_SEC - tstart << std::endl; // OLD Implementation diff --git a/test/benchmarks_cg.cpp b/test/benchmarks_cg.cpp index 5c2d51ad3..f73e14eaa 100644 --- a/test/benchmarks_cg.cpp +++ b/test/benchmarks_cg.cpp @@ -77,18 +77,18 @@ int main() tstart = (double)clock()/(double)CLOCKS_PER_SEC; std::cout << "BallWalk (cube) = " - << volume_gaussian_annealing(HP, e, walk_len) << " , "; + << volume_gaussian_annealing(HP, e, walk_len) << " , "; std::cout << (double)clock()/(double)CLOCKS_PER_SEC - tstart << std::endl; tstart = (double)clock()/(double)CLOCKS_PER_SEC; tstart = (double)clock()/(double)CLOCKS_PER_SEC; std::cout << "CDHRWalk (cube) = " - << volume_gaussian_annealing(HP, e, walk_len) << " , "; + << volume_gaussian_annealing(HP, e, walk_len) << " , "; std::cout << (double)clock()/(double)CLOCKS_PER_SEC - tstart << std::endl; tstart = (double)clock()/(double)CLOCKS_PER_SEC; std::cout << "RDHRWalk (cube) = " - << volume_gaussian_annealing(HP, e, walk_len) << " , "; + << volume_gaussian_annealing(HP, e, walk_len) << " , "; std::cout << (double)clock()/(double)CLOCKS_PER_SEC - tstart << std::endl; // OLD Implementation diff --git a/test/benchmarks_sob.cpp b/test/benchmarks_sob.cpp index 9520c72a9..6fc65bf25 100644 --- a/test/benchmarks_sob.cpp +++ b/test/benchmarks_sob.cpp @@ -78,17 +78,17 @@ int main() tstart = (double)clock()/(double)CLOCKS_PER_SEC; std::cout << "CDHRWalk (cube) = " - << volume_sequence_of_balls(HP, e, walk_len) << " , "; + << volume_sequence_of_balls(HP, e, walk_len) << " , "; std::cout << (double)clock()/(double)CLOCKS_PER_SEC - tstart << std::endl; tstart = (double)clock()/(double)CLOCKS_PER_SEC; std::cout << "RDHRWalk (cube) = " - << volume_sequence_of_balls(HP, e, walk_len) << " , "; + << volume_sequence_of_balls(HP, e, walk_len) << " , "; std::cout << (double)clock()/(double)CLOCKS_PER_SEC - tstart << std::endl; tstart = (double)clock()/(double)CLOCKS_PER_SEC; std::cout << "BilliardWalk (cube) = " - << volume_sequence_of_balls(HP, e, walk_len) << " , "; + << volume_sequence_of_balls(HP, e, walk_len) << " , "; std::cout << (double)clock()/(double)CLOCKS_PER_SEC - tstart << std::endl; std::cout << std::endl; @@ -141,25 +141,25 @@ int main() VP.init(VP.dimension(), VP.get_mat(), VP.get_vec()); tstart = (double)clock()/(double)CLOCKS_PER_SEC; std::cout << "Ball (cross) = " - << volume_sequence_of_balls(VP) << " , "; + << volume_sequence_of_balls(VP) << " , "; std::cout << (double)clock()/(double)CLOCKS_PER_SEC - tstart << std::endl; VP.init(VP.dimension(), VP.get_mat(), VP.get_vec()); tstart = (double)clock()/(double)CLOCKS_PER_SEC; std::cout << "RDHR (cross) = " - << volume_sequence_of_balls(VP) << " , "; + << volume_sequence_of_balls(VP) << " , "; std::cout << (double)clock()/(double)CLOCKS_PER_SEC - tstart << std::endl; VP.init(VP.dimension(), VP.get_mat(), VP.get_vec()); tstart = (double)clock()/(double)CLOCKS_PER_SEC; std::cout << "CDHR (cross) = " - << volume_sequence_of_balls(VP) << " , "; + << volume_sequence_of_balls(VP) << " , "; std::cout << (double)clock()/(double)CLOCKS_PER_SEC - tstart << std::endl; VP.init(VP.dimension(), VP.get_mat(), VP.get_vec()); tstart = (double)clock()/(double)CLOCKS_PER_SEC; std::cout << "Blrd (cross) = " - << volume_sequence_of_balls(VP) << " , "; + << volume_sequence_of_balls(VP) << " , "; std::cout << (double)clock()/(double)CLOCKS_PER_SEC - tstart << std::endl; // OLD Implementation diff --git a/test/new_rounding_test.cpp b/test/new_rounding_test.cpp index 0a7c8c517..95d84fd79 100644 --- a/test/new_rounding_test.cpp +++ b/test/new_rounding_test.cpp @@ -9,6 +9,10 @@ #include #include #include "misc.h" +#include "random.hpp" +#include "random/uniform_int.hpp" +#include "random/normal_distribution.hpp" +#include "random/uniform_real_distribution.hpp" #include "new_volume.hpp" #include "new_gaussian_volume.hpp" #include "new_cooling_balls.hpp" @@ -51,7 +55,7 @@ void rounding_test(Polytope &HP, RNGType rng(d); std::pair InnerBall = HP.ComputeInnerBall(); - std::pair< std::pair, NT > res = round_polytope(HP, InnerBall, 10 + 10 * d, rng); + std::pair< std::pair, NT > res = round_polytope(HP, InnerBall, 10 + 10 * d, rng); // Setup the parameters int walk_len = 1; diff --git a/test/new_volume_example.cpp b/test/new_volume_example.cpp index 72a6922da..642df1da7 100644 --- a/test/new_volume_example.cpp +++ b/test/new_volume_example.cpp @@ -36,13 +36,19 @@ int main() typedef HPolytope Hpolytope; typedef VPolytope Vpolytope; + + typedef BoostRandomNumberGenerator RNG; Hpolytope HPoly = gen_cube(10, false); BallWalk BW(3); // Estimate the volume + + //VPolytope VP2 = VP; + + //VP.init(VP.dimension(), VP2.get_mat(), VP2.get_vec()); double tstart = (double)clock()/(double)CLOCKS_PER_SEC; - std::cout << "Ball = " - << volume_sequence_of_balls<>(HPoly, BW) << " , "; + std::cout << "Ball (cross) = " + << volume_cooling_balls(HPoly) << " , "; std::cout << (double)clock()/(double)CLOCKS_PER_SEC - tstart << std::endl; return 0; diff --git a/test/volume_cb_hpolytope.cpp b/test/volume_cb_hpolytope.cpp index a87f89a68..e62c6198f 100644 --- a/test/volume_cb_hpolytope.cpp +++ b/test/volume_cb_hpolytope.cpp @@ -9,6 +9,10 @@ #include #include #include "misc.h" +#include "random.hpp" +#include "random/uniform_int.hpp" +#include "random/normal_distribution.hpp" +#include "random/uniform_real_distribution.hpp" #include "new_volume.hpp" #include "new_gaussian_volume.hpp" #include "new_cooling_balls.hpp" diff --git a/test/volume_cb_vpoly_intersection_vpoly.cpp b/test/volume_cb_vpoly_intersection_vpoly.cpp index 195a419a4..cc88bd046 100644 --- a/test/volume_cb_vpoly_intersection_vpoly.cpp +++ b/test/volume_cb_vpoly_intersection_vpoly.cpp @@ -9,6 +9,10 @@ #include #include #include "misc.h" +#include "random.hpp" +#include "random/uniform_int.hpp" +#include "random/normal_distribution.hpp" +#include "random/uniform_real_distribution.hpp" #include "new_volume.hpp" #include "new_gaussian_volume.hpp" #include "new_cooling_balls.hpp" diff --git a/test/volume_cb_vpolytope.cpp b/test/volume_cb_vpolytope.cpp index ccff704d8..535c42e1c 100644 --- a/test/volume_cb_vpolytope.cpp +++ b/test/volume_cb_vpolytope.cpp @@ -9,6 +9,10 @@ #include #include #include "misc.h" +#include "random.hpp" +#include "random/uniform_int.hpp" +#include "random/normal_distribution.hpp" +#include "random/uniform_real_distribution.hpp" #include "new_volume.hpp" #include "new_gaussian_volume.hpp" #include "new_cooling_balls.hpp" diff --git a/test/volume_cb_zonotopes.cpp b/test/volume_cb_zonotopes.cpp index 08178076f..73f0a5646 100644 --- a/test/volume_cb_zonotopes.cpp +++ b/test/volume_cb_zonotopes.cpp @@ -9,6 +9,10 @@ #include #include #include "misc.h" +#include "random.hpp" +#include "random/uniform_int.hpp" +#include "random/normal_distribution.hpp" +#include "random/uniform_real_distribution.hpp" #include "new_volume.hpp" #include "new_gaussian_volume.hpp" #include "new_cooling_balls.hpp" diff --git a/test/volume_cg_hpolytope.cpp b/test/volume_cg_hpolytope.cpp index b64dccef5..6b9b89d6a 100644 --- a/test/volume_cg_hpolytope.cpp +++ b/test/volume_cg_hpolytope.cpp @@ -9,6 +9,10 @@ #include #include #include "misc.h" +#include "random.hpp" +#include "random/uniform_int.hpp" +#include "random/normal_distribution.hpp" +#include "random/uniform_real_distribution.hpp" #include "new_volume.hpp" #include "new_gaussian_volume.hpp" #include "known_polytope_generators.h" diff --git a/test/volume_cg_vpolytope.cpp b/test/volume_cg_vpolytope.cpp index 74896da3d..be04ccacc 100644 --- a/test/volume_cg_vpolytope.cpp +++ b/test/volume_cg_vpolytope.cpp @@ -9,6 +9,10 @@ #include #include #include "misc.h" +#include "random.hpp" +#include "random/uniform_int.hpp" +#include "random/normal_distribution.hpp" +#include "random/uniform_real_distribution.hpp" #include "new_volume.hpp" #include "new_gaussian_volume.hpp" #include "new_cooling_balls.hpp" diff --git a/test/volume_sob_hpolytope.cpp b/test/volume_sob_hpolytope.cpp index 7ee35dc3d..354e50081 100644 --- a/test/volume_sob_hpolytope.cpp +++ b/test/volume_sob_hpolytope.cpp @@ -9,6 +9,10 @@ #include #include #include "misc.h" +#include "random.hpp" +#include "random/uniform_int.hpp" +#include "random/normal_distribution.hpp" +#include "random/uniform_real_distribution.hpp" #include "new_volume.hpp" #include "known_polytope_generators.h" diff --git a/test/volume_sob_vpolytope.cpp b/test/volume_sob_vpolytope.cpp index 18d168ab3..a829bb574 100644 --- a/test/volume_sob_vpolytope.cpp +++ b/test/volume_sob_vpolytope.cpp @@ -9,6 +9,10 @@ #include #include #include "misc.h" +#include "random.hpp" +#include "random/uniform_int.hpp" +#include "random/normal_distribution.hpp" +#include "random/uniform_real_distribution.hpp" #include "new_volume.hpp" #include "known_polytope_generators.h"